////////////////////////////////////////
//                                    //
// Contre-exemple en optimisation de  //
// l'épaisseur d'une plaque élastique //
//                                    //
// Ecole Polytechnique, MAP 562       //
// Copyright G. Allaire, 2004, 2008   //
//                                    //
////////////////////////////////////////
ofstream object("contrex.obj") ;
ofstream grad("contrex.grad") ;
string sauve="contrex";		// Nom du fichier de sauvegarde
int niter=100;			// Nombre d'itérations
int n;				// Taille du maillage
real pas=1. ;			// Pas de descente
real lagrange=0.01;		// Multiplicateur de Lagrange pour le volume
real lagmin, lagmax ;		// Encadrement du multiplicateur de Lagrange
int inddico ;
real objectif;			// Fonction cout
real objectifold ;
real volume0;			// Volume initial
real volume,volume1;		// Volume de la forme courante
real exa=1;			// Coefficient d'exagération de la déformation
string legende;			// Légende pour les sorties graphiques
real E=100;			// Module de Young (toujours positif)
real nu=0.3;			// Coefficient de Poisson (toujours compris entre -1 et 1/2)
func g1=0;  			// Forces appliquées
func g2=100;

real[int] vviso(21);
for (int i=0;i<21;i++)
vviso[i]=i*0.05 ;
/////////////////////////////////////
// Calcul des coefficients de Lamé //
/////////////////////////////////////
real lambda=E*nu/((1.+nu)*(1.-2.*nu));
real mu=E/(2.*(1.+nu));

/////////////////////////////////////////
// Epaisseur mini et maxi de la plaque //
/////////////////////////////////////////
real hmin=0.1;
real hmax=1.0;
real hmoy=0.5;
func h0=hmoy ;
//func b1=(y>8./9.);
//func b2=(8./9.>y)*(y>7./9.);
//func b3=(7./9.>y)*(y>6./9.);
//func b4=(6./9.>y)*(y>5./9.);
//func b5=(5./9.>y)*(y>4./9.);
//func b6=(4./9.>y)*(y>3./9.);
//func b7=(3./9.>y)*(y>2./9.);
//func b8=(2./9.>y)*(y>1./9.);
//func b9=(1./9.>y);
//func h0= 0.1*(b1+b3+b5+b7+b9) + hmax*(b2+b4+b6+b8) ;

/////////////////////////////
// Définition du domaine   // 
/////////////////////////////
mesh Th;

border bd(t=0,1)       { x=+1; y=t;label=2; };   // bord droit de la forme
border bg(t=1,0)       { x=-1; y=t;label=1; };   // bord gauche de la forme
border bs(t=1,-1)      { x=t;  y=1;label=3; };   // bord supérieur de la forme
border bi(t=-1,1)      { x=t;  y=0;label=4; };   // bord inférieur de la forme

//////////////////////////////
// Construction du maillage //
//////////////////////////////
n = 2 ;
Th= buildmesh (bd(10*n)+bs(20*n)+bg(10*n)+bi(20*n));
plot(Th,wait=0);
savemesh(Th,sauve+".msh");

fespace Vh0(Th,P0);				//Définition de l'espace P0
fespace Vh2(Th,[P2,P2]);			//Définition de l'espace P2*P2

Vh2 [u,v],[w,s],[p,q],[uold,vold];
Vh0 h,hold,gradient;

h = h0 ;

///////////////////////
// Déplacement cible //
///////////////////////
real u0x=-1.;
real u0y=100.;
real cx=1.;
real cy=0.1;
func u0=u0x;
func v0=u0y;
///////////////////////////
// Définition du système //
///////////////////////////
problem elasticite([u,v],[w,s]) =
    int2d(Th)(
               2.*h*mu*(dx(u)*dx(w)+dy(v)*dy(s)+((dx(v)+dy(u))*(dx(s)+dy(w)))/2.)
              +h*lambda*(dx(u)+dy(v))*(dx(w)+dy(s))
	)
    -int1d(Th,2)(g2*w+g1*s)
    -int1d(Th,3)(g1*w+0.1*g2*s)
    -int1d(Th,4)(g1*w-0.1*g2*s)
    +on(1,u=0,v=0)
;
/////////////////////////////
// Définition de l'adjoint //
/////////////////////////////
problem adjoint([p,q],[w,s]) =
    int2d(Th)(
               2.*h*mu*(dx(p)*dx(w)+dy(q)*dy(s)+((dx(q)+dy(p))*(dx(s)+dy(w)))/2.)
              +h*lambda*(dx(p)+dy(q))*(dx(w)+dy(s))
	)
    +int1d(Th,2)( 2*cx*(u-u0)*w )
    +int1d(Th,3)( 2*cy*(v-v0)*s )
    +int1d(Th,4)( 2*cy*(v+v0)*s )
    +on(1,p=0,q=0)
;
//////////////////////////////
// Calcul du volume initial //
//////////////////////////////
volume0=int2d(Th)(h);

////////////////////////////
// Fonction cout initiale //
////////////////////////////
elasticite;
objectif=int1d(Th,2)( cx*(u-u0)*(u-u0) )
        +int1d(Th,3)( cy*(v-v0)*(v-v0) )
        +int1d(Th,4)( cy*(v+v0)*(v+v0) ) ;
object << objectif   << endl ;
cout<<"Initialisation. Cout: "<<objectif<<" Volume: "<<volume0<<endl;

//////////////////////////////////
// Adjoint et gradient initial  //
//////////////////////////////////

adjoint ;

gradient = 2*mu*(dx(u)*dx(p)+dy(v)*dy(q)+((dx(v)+dy(u))*(dx(q)+dy(p)))/2.)+lambda*(dx(u)+dy(v))*(dx(p)+dy(q)) ;

real norgrad=int2d(Th)(gradient*gradient);
norgrad=sqrt(norgrad) ;
grad << norgrad   << endl ;
pas=pas/norgrad ;
cout << "premier pas "<< pas << endl;


///////////////////////////////////
//     Boucle d'optimisation     //
///////////////////////////////////
 
int iter;

for (iter=1;iter< niter;iter=iter+1)  
{
cout <<"Iteration " <<iter <<" ----------------------------------------" <<endl;

hold = h ;
objectifold = objectif ;
[uold,vold] = [u,v] ;

h = hold - pas*gradient - lagrange ;
h = min(h,hmax) ;
h = max(h,hmin) ;

///////////////////////////////////////////////////////////
// Calcul d'un encadrement du multiplicateur de Lagrange //
///////////////////////////////////////////////////////////
volume=int2d(Th)(h);
volume1 = volume ;
//
if (volume1 < volume0)
{
   while (volume < volume0)
{ 
      lagmin = lagrange ;
      lagrange = lagrange - 0.01 ;
      h = hold - pas*gradient - lagrange ;
      h = min(h,hmax) ;
      h = max(h,hmin) ;
      volume=int2d(Th)(h) ; 
};
   lagmax = lagrange ;
};
//
if (volume1 > volume0)
{
   while (volume > volume0)
{
      lagmax = lagrange ;
      lagrange = lagrange + 0.01 ;
      h = hold - pas*gradient - lagrange ;
      h = min(h,hmax) ;
      h = max(h,hmin) ;
      volume=int2d(Th)(h) ;
};
   lagmin = lagrange ;
};
//
if (volume1 == volume0) 
{
   lagmin = lagrange ;
   lagmax = lagrange ;
};

//////////////////////////////////////////////////
// Dichotomie sur le multiplicateur de lagrange //
//////////////////////////////////////////////////

inddico=0;
while ((abs(1.-volume/volume0)) > 0.0000001)
{
   lagrange = (lagmax+lagmin)*0.5 ;
   h = hold - pas*gradient - lagrange ;
   h = min(h,hmax) ;
   h = max(h,hmin) ;
   volume=int2d(Th)(h) ;
   inddico=inddico+1;
   if (volume < volume0) 
      {lagmin = lagrange ;} ;
   if (volume > volume0)
      {lagmax = lagrange ;} ;
};

//cout<<"nbre iterations dichotomie: "<<inddico<<endl;

// Résolution du système de l'élasticité linéaire
elasticite;

//Calcul de la fonction cout
objectif=int1d(Th,2)( cx*(u-u0)*(u-u0) )
        +int1d(Th,3)( cy*(v-v0)*(v-v0) )
        +int1d(Th,4)( cy*(v+v0)*(v+v0) ) ;

if (objectif > objectifold*1.00001 )
{
// on refuse le pas
pas = pas / 2. ;
cout << "refus, nouveau pas " << pas << endl ;
h=hold ;
[u,v] = [uold,vold] ;
objectif = objectifold ;
}
else
{
// on accepte le pas
pas = pas * 1.1 ;
cout << "nouveau pas " << pas << endl ;

///////////////////////////////////////////////
// Calcul du gradient grâce à l'adjoint      //
///////////////////////////////////////////////
adjoint ;
gradient = 2*mu*(dx(u)*dx(p)+dy(v)*dy(q)+((dx(v)+dy(u))*(dx(q)+dy(p)))/2.)+lambda*(dx(u)+dy(v))*(dx(p)+dy(q)) ;
//legende="Iteration "+iter+", gradient";
//plot(Th,gradient,fill=1,value=true,viso=vviso,cmm=legende,wait=1); 

};

//////////////////////////////////////////////////
// Fin de la boucle, affichage et impressions   //
//////////////////////////////////////////////////

object << objectif   << endl ;
cout<<"Cout: "<<objectif<<" Volume: "<<volume<<" Lagrange: "<<lagrange<<endl;
legende="Iteration "+iter+", Cout "+objectif+", Volume="+volume;
plot(Th,h,fill=1,value=0,viso=vviso,cmm=legende,wait=0); 

};



//On sauve la forme finale
legende="Forme finale, Iteration "+iter+", Cout "+objectif+", Volume="+volume;
plot(Th,h,fill=1,value=0,viso=vviso,cmm=legende,ps=sauve+".eps");
{ ofstream file(sauve+".bb");
	file << h[].n << " \n";
	int j;
	for (j=0;j<h[].n ; j++)  
	file << h[][j] << endl;  }  	
//exec("xd3d -bord=2 -hidden -fich="+sauve+".bb -iso=11 -nbcol=9 -table=8 "+sauve);

/////////////////////////////
// Déformation du maillage //
/////////////////////////////
//real coef = 0.1 ;
//Th = movemesh (Th,[x+coef*u,y+coef*v]);
//
//plot(Th,wait=1); 
//[u,v] = [u,v] ;
//
//plot([u,v],wait=1,ps=sauve+"def.eps");