// Solve linear elasticity problem (34 lines, comments and blank lines included)

real L=4,H=1,dn=10;	// Length/Height of the domain + nb of dof per unit Length
real lambda=1,mu=1;	// Lame Moduli

func g1=0;		// surface loads
func g2=-1;
macro g [g1,g2]		// end of surface loads defintion

macro Neumann 2		// Label(s) of parts boundary where Loads g are applied
macro Dirichlet 4	// Label(s) of compled part of the boundary

//---------------------------------------------------------------------------//
mesh Th=square(L*dn,H*dn,[L*x,H*y]);

macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/sqrt(2.),dy(u[1])] //
macro div(u) (dx(u[0])+dy(u[1]))//

fespace Vh(Th,[P1,P1]);

macro u [u1,u2]//
macro v [v1,v2]//

Vh u,v;
// Problem definition
problem elasticity(u,v)=int2d(Th)(2.*mu*(e(u)'*e(v))+lambda*div(u)*div(v))-int1d(Th,Neumann)(g'*v)+on(Dirichlet,u1=0,u2=0);

// Solving the problem
elasticity;

// Output
real exa=0.01;
mesh Sh=movemesh(Th,[x+exa*u1,y+exa*u2]);
plot(Sh,cmm="Deformed structure");