////////////////////////////////////////////////
// Visualization of the optimal Grip          //
// Copyright O.Pantz, G. Allaire (2005)       //
////////////////////////////////////////////////
string backupname="grip";	//Name under which the results are stored
// Material parameters
	real E=15;				//Young Modulus (always positif)
	real nu=0.4;				//Poisson ratio (between -1 and 1/2)
	real lambda=E*nu/((1.+nu)*(1.-2.*nu));	//Lame Moduli
	real mu=E/(2.*(1.+nu));

//Applied Loads
	func g1=1.; func g2=0.;

////////////////////////////
// Definition of the mesh //
////////////////////////////
	//Boundaries with (non null)-Neumann and Dirichlet conditions are not optimized
		int neumann=1; int dirichlet=2; int teeth=4;
	//Boundaries with free conditions are optimized
		int free=3;
	mesh Sh; Sh=readmesh(backupname+".msh");
	plot(Sh,wait=1);
//Finite elements spaces on the finer mesh
	fespace WSh(Sh,[P2,P2]);
//Finite Elements
	WSh [ut1,ut2],[v1,v2];

/////////////////////////
// Elasticity Problem  //
/////////////////////////
problem elasticity2([ut1,ut2],[v1,v2]) =
	int2d(Sh)(
		2.*mu*(dx(ut1)*dx(v1)+dy(ut2)*dy(v2)+((dx(ut2)+dy(ut1))*(dx(v2)+dy(v1)))/2.)
		+lambda*(dx(ut1)+dy(ut2))*(dx(v1)+dy(v2)))
	-int1d(Sh,neumann)(g1*v1+g2*v2)	           
	+on(dirichlet,ut1=0,ut2=0);

elasticity2;
plot(Sh,wait=1);

real examax=0.4;
mesh Shdep;
real exa=0.;
for (int i=0;i<5;i++){
	exa=exa+examax/10.;
	Shdep=movemesh(Sh,[x+exa*ut1,y+exa*ut2]);
	plot(Shdep,wait=1);}

plot(Shdep,ps="deformed-grip.eps",wait=1);
plot(Sh,ps="reference-grip.eps",wait=1);