// Resolution d'un probleme d'elasticite

///////////////////////////

real E=15.;			//Module de Young (toujours positif)
real nu=0.35;			//Coefficient de Poisson (toujours compris entre -1 et 1/2)
real lambda=E*nu/((1.+nu)*(1.-2.*nu));//coef de Lame
real mu=E/(2.*(1.+nu));

func g1=0;  			//Forces appliquées 
func g2=-1;

func f1=0;
func f2=-1;
////////////////////////////

real exa=0.05;

int Nbnoeuds=10;

mesh Th=square(4*Nbnoeuds,Nbnoeuds,[4.*x,y]);
mesh Sh;

fespace Vh(Th,P1);
fespace Wh(Th,[P1,P1]);	//definition de l'espace 


Wh [uh1,uh2],[vh1,vh2];	//definition de l'inconnu et des fonctions tests

problem elasticite([uh1,uh2],[vh1,vh2],solver=LU) =
    int2d(Th)(
               2.*mu*(dx(uh1)*dx(vh1)+dy(uh2)*dy(vh2)+((dx(uh2)+dy(uh1))*(dx(vh2)+dy(vh1)))/2.)
              +lambda*(dx(uh1)+dy(uh2))*(dx(vh1)+dy(vh2))
	)
    -int1d(Th,2)(g1*vh1+g2*vh2)	           
	-int2d(Th)(f1*vh1+f2*vh2)
    +on(4,uh1=0,uh2=0)
;


elasticite;

Vh sigma;

sigma=((2*mu+lambda)*dx(uh1)+lambda*dy(uh2))^2
+((2*mu+lambda)*dy(uh2)+lambda*dx(uh1))^2
+2*mu^2*(dx(uh2)+dy(uh1))^2;

Sh=movemesh(Th,[x+exa*uh1,y+exa*uh2]);

plot(sigma,value=1,wait=1);
plot(Sh,wait=1);

real erreur=0.001;
Th=adaptmesh(Th,uh1,uh2,err=erreur);

elasticite;


sigma=((2*mu+lambda)*dx(uh1)+lambda*dy(uh2))^2
+((2*mu+lambda)*dy(uh2)+lambda*dx(uh1))^2
+2*mu^2*(dx(uh2)+dy(uh1))^2;

Sh=movemesh(Th,[x+exa*uh1,y+exa*uh2]);

plot(sigma,value=1,wait=1);
plot(Sh,wait=1);