////////////////////////////////////
//                                //
// Elasticite : contre-exemple    //
// au principe du maximum         //
// (visualisation par xd3d)       //
//                                //
// Ecole Polytechnique, MAP 431   //
// Copyright G. Allaire, 2003     //
//                                //
////////////////////////////////////
// Nom du fichier de sauvegarde
// des figures postscript
string sauve="elasticite";
/////////////////////////////////////
real energie ;
string legende ;
string fichier, fichier2 ;
//Module de Young (toujours positif)
real E=1;
//Coefficient de Poisson (toujours compris entre -1 et 1/2)
real nu=0.3;
// Coefficients de Lamé 
real lambda=E*nu/((1.+nu)*(1.-2.*nu));
real mu=E/(2.*(1.+nu));
/////////////////////////////
// Définition du domaine   //
/////////////////////////////
//  bord superieur 
border a1(t=2,0) { x=t; y=1;label=1; };
// bord gauche 
border a2(t=1,0)  { x=0; y=t;label=2; };
// bord inferieur
border a3(t=0,2) { x=t; y=0;label=3; };
// bord droit
border a4(t=0,1) { x=2; y=t;label=4; };
//////////////////////////////
// construction du maillage //
//////////////////////////////
int n=5 ;
mesh Sh;
Sh= buildmesh (a1(2*n)+a2(n)+a3(2*n)+a4(n));
//
plot(Sh,wait=1,bw=1,ps=sauve+".maillage.eps"); 
fichier=sauve+".msh";
savemesh(Sh,fichier);
////////////////////////////////
// Définition de l'espace P1  //
////////////////////////////////
fespace Vh1(Sh,[P1,P1]);
Vh1 [u,v] , [u1,v1] ;	
//Définition de l'espace P0
fespace Vh0(Sh,P0);
Vh0 gradient;
//////////////////////////////////////////////////
// Définition de la formulation variationnelle  //
//////////////////////////////////////////////////
problem elasticite([u,v],[u1,v1]) =
        int2d(Sh)( 2*mu*( dx(u)*dx(u1)+dy(v)*dy(v1) 
                 + 0.5*((dx(v)+dy(u))*(dx(v1)+dy(u1))) )
                 + lambda*(dx(u)+dy(v))*(dx(u1)+dy(v1)) )
       +on(2,u=0,v=0)+on(4,u=1)
;
///////////////////////////
// Résolution du système //
///////////////////////////
elasticite ;
//
gradient=2*mu*(dx(u)^2+dy(v)^2+((dx(v)+dy(u))^2)/2.)+lambda*(dx(u)+dy(v))^2;
//
plot([u,v],wait=1,value=true); 
/////////////////////////////
// Déformation du maillage //
/////////////////////////////
real pas = 2. ;
Sh = movemesh (Sh,[x+pas*u,y+pas*v]);
//
plot(Sh,wait=1); 
[u,v] = [u,v] ;
//
plot([u,v],wait=1,ps=sauve+".eps"); 
fichier=sauve+".def.msh";
savemesh(Sh,fichier);
fichier2=sauve+".bb";
int taille=(u[].n)/2 ;
{ ofstream file(sauve+".bb");
	file << taille << " \n";
	int j;
	for (j=0;j<taille ; j++)  
	file << u[][2*j] << " " << u[][2*j+1] << endl;  }  	
exec("xd3d -bord=3 -hidden -fich="+fichier2+" -fleche "+fichier);