// file Lame.edp
//
// SOLUTION OF THE LAME SYSTEM OF ELASTICITY:
//
// 

//
// THE DATA:
//

real E = 21.e5, nu = 0.28;                              // Young modulus, Poisson ratio
real mu= E/(2*(1+nu)), lambda = E*nu/((1+nu)*(1-2*nu)); // Lame coefficients

func f = -1.;                                           // External forces
func g = .1*x;                                          // Applied surface forces

//
// AUXILIAR VARIABLES:
//

real sqrt2=sqrt(2.), usqrt2=1./sqrt2, mu2=2.*mu, scale=100.;

//
// THE MESH:
//

mesh Th=square(20,10,[20*x,2*y-1]);    // The domain: (0.,20) X (-1.,1.)
// mesh Th=square(100,50,[20*x,2*y-1]);    // The domain: (0.,20) X (-1.,1.)

//
// THE FINITE ELEMENT SPACE:
//

fespace Vh(Th,P1);
// fespace Vh(Th,P2);
Vh u,v,uu,vv;

//
// SOME MACRO DEFINITIONS:
//

macro epsilon(u1,u2) [dx(u1),dy(u2),usqrt2*(dy(u1)+dx(u2))] // EOM
		
macro div(u,v) (dx(u)+dy(v)) // EOM

//
// THE PROBLEM:
//

problem lame([u,v],[uu,vv]) = int2d(Th)(lambda*div(u,v)*div(uu,vv)
        + mu2*(epsilon(u,v)'*epsilon(uu,vv)))
        - int2d(Th)(f*vv) - int1d(Th,1)(g*vv)
        + on(4,u=0.,v=0.);

//
// THE COMPUTATIONS:
//

lame;

//
// FIGURES:
//

plot([u,v],wait=1,coef=scale,cmm="THE DISPLACEMENTS");
//plot([u,v],wait=1,ps="lamevect.eps",coef=scale,cmm="THE DISPLACEMENTS");
plot(u,wait=1,fill=1,value=1,cmm="THE x-DISPLACEMENT");
plot(v,wait=1,fill=1,value=1,cmm="THE y-DISPLACEMENT");
mesh Thd = movemesh(Th,[x+scale*u, y+scale*v]);
//plot(Thd,wait=1,ps="lamedeform.eps",cmm="THE DEFORMED MESH");
plot(Th,Thd,wait=1,cmm="THE DEFORMED MESH");

real dxmin  = u[].min; 
real dymin  = v[].min; 

cout << " min displacements:    x = "<< dxmin<< " y=" << dymin << endl; 
cout << " displacement at (20,0)  = "<< u(20,0) << " " << v(20,0) << endl;