// elastic_bar.edp
//
//  Discussion:
//
//    This code models a case of linear elasticity in 2D.
//
//  Location:
//
//    https://people.sc.fsu.edu/~jburkardt/freefem_src/elastic_bar/elastic_bar.edp  
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    27 February 2021
//
//  Author:
//
//    Roberto Font, Francisco Periago
//
//  Reference:
//
//    Roberto Font, Francisco Periago,
//    The Finite Element Method with FreeFem++ for Beginners,
//    Electronic Journal of Mathematics and Technology,
//    Volume 7, Number 4, June 2013, pages 289-307.
//
cout << endl;
cout << "elastic_bar:" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Model the deformation of an elastic bar under stress." << endl;
//
//  Domain:
//  Lower side, right side, upper side, left side, interior circle.
//
border a1(t=0,20){x=t;y=-1;label=1;};
border a2(t=-1, 1){x=20;y=t;label=2;};
border a3(t=0,20){x=20-t;y=1;label=3;}; 
border a4(t=-1,1){x=0;y=-t;label=4;};
border C(t=0, 2*pi){x=10+0.5*cos(t); y=0.5*sin(t);};
//
//  Create and plot the mesh.
//
mesh Th=buildmesh(a1(100)+a2(10)+a3(100)+a4(10)+C(-10));
plot ( Th, wait = true, ps = "elastic_bar_mesh.ps" );
//
//  Define a piecewise quadratic finite element space on the mesh.
//
fespace Vh(Th,P2);
Vh u,v,uu,vv;
//
//  Define some variables:
//    E: Young's modulus.
//    f: body force
//    lambda: Lame coefficient
//    mu: Lame coefficient
//    nu: Poisson's ratio
//    sqrt2: square root of 2
//
real sqrt2=sqrt(2.);
real E = 21e5;
real nu = 0.28;
real lambda = E*nu/((1+nu)*(1-2*nu));
real mu= E/(2*(1+nu));
real f = -1;
//
//  We can define macros to make programming easier
//
macro epsilon(u,v) [dx(u),dy(v),(dy(u)+dx(v))/sqrt2] // EOM
macro div(u,v) (dx(u)+dy(v)) // EOM
//
//  Solving the problem in variational form
//
solve lame ( [u,v],[uu,vv] ) = 
  int2d ( Th )
  (
    lambda*div(u,v)*div(uu,vv)+2.0*mu*( epsilon(u,v)'*epsilon(uu,vv) ) 
  )
  - int2d ( Th ) ( f * vv )
  + on ( 4, u = 0, v = 0 );
//
//  Plotting the result, the displacement of the bar.
//
real coef = 100; //A coefficient of amplification
mesh Thd = movemesh ( Th, [x+u*coef, y+v*coef] );
plot ( Thd, wait = true,
  ps = "elastic_bar_deformed.ps" );

real dxmin = u[].min; //Minimum of displacement in x direction
real dymin = v[].min; //Minimum of displacement in y direction
//
//  The von Mises stress
//
//  Define a second finite element space.
//
fespace Wh(Th,P1);
Wh Sigmavm;
//
//  Stress tensor (since it is symmetric, it is enough with 3 elements)
//
macro Sigma(u,v) [2*mu*dx(u)+lambda*(dx(u)+dy(v)),
2*mu*dy(v)+lambda*(dx(u)+dy(v)),mu*(dy(u)+dx(v))]//EOM
//
//  Von Mises stress
//
Sigmavm = sqrt(Sigma(u,v)[0]*Sigma(u,v)[0]-Sigma(u,v)[0]*Sigma(u,v)[1]
+Sigma(u,v)[1]*Sigma(u,v)[1]+3*Sigma(u,v)[2]*Sigma(u,v)[2]);
//
//  Plot the stress.
//
plot ( Sigmavm, wait = true, fill = true, 
  ps = "elastic_bar_stress.ps" );

real Sigmavmmax = Sigmavm[].max;
cout << "  maximum von Mises stress = " << Sigmavmmax << endl;
//
//  Terminate.
//
cout << "\n";
cout << "elastic_bar:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );