/////////////////////////////////////////////////////
// Méthode de convexification pour la minimisation //
//        de la compliance d'une console           //
//                                                 //
//         Ecole Polytechnique, MAP 562            //
//         Copyright G. Allaire, 2004              //
/////////////////////////////////////////////////////

ofstream fichobj("console.obj") ;
string sauve="console";		// nom du fichier de sauvegarde
int niternp=30;			// nombre d'itérations non pénalisées
int niterpp=20;			// nombre d'itérations pénalisées
int niter=niternp+niterpp;	// nombre d'itérations total
int n=4;			// Taille du maillage
real compliance;		// Compliance
string legende;			// Légende pour les sorties graphiques
real E=1;			// Module de Young
real nu=0.3;			// Coefficient de Poisson
real lagrange=1;		// Multiplicateur de Lagrange pour le volume
real lagmin, lagmax;		// Encadrement du multiplicateur de Lagrange
real masse0;			// Masse initiale
real masse,masse1;		// Masse de la forme courante
int inddico ;			// Indice dans la recherche du multiplicateur de Lagrange
real thetamin=0.01;		// Densité minimale de matière
real thetamoy=0.5;		// Densité moyenne de matière
real pi=4*atan(1);		// Valeur de pi=3.14159...
func f1=0;  			// Forces appliquées
func f2=-1;

// Définition des isovaleurs de tracé
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));

///////////////////////////// 			1:Condition de Dirichlet
// Définition du domaine   // 			2:Condition Libre
///////////////////////////// 			3:Condition de Neuman non nulle
mesh Th;
border bd1(t=0,0.45)		{ x=+1; y=t;label=2; };   // bord droit de la forme
border bd2(t=0.45,0.55)		{ x=1;  y=t;label=3; };
border bd3(t=0.55,1)		{ x=1;  y=t;label=2; };
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=2; };   // bord supérieur de la forme
border bi(t=-1,1)		{ x=t;  y=0;label=2; };   // bord inférieur de la forme

//////////////////////////////
// Construction du maillage //
//////////////////////////////
int m=int(4.5*n) ;
Th= buildmesh (bg(10*n)+bs(20*n)+bi(20*n)+bd1(m)+bd2(n)+bd3(m));
plot(Th,wait=0,ps=sauve+".maillage.eps");
savemesh(Th,sauve+".msh");

fespace Vh0(Th,P0);
fespace Vh2(Th,P2);
//fespace Vh2(Th,[P1,P1]);
Vh2 u,v,w,s;
//Vh2 [u,v],[w,s];
Vh0 coef,eu,ev,euv,theta,intermed;
// theta = densité de matière
// coef  = theta lors des itérations non pénalisées
// coef  = 0.5*(1.-cos(pi*theta)) lors des itérations pénalisées

func theta0=thetamoy;
theta=theta0;
coef=theta;
// Masse initiale
real volume=int2d(Th)(1.);
masse0 = int2d(Th)(theta);
masse0 = masse0/volume ;

///////////////////////////
// Définition du système //
///////////////////////////
problem elasticite([u,v],[w,s],solver=Crout) =
    int2d(Th)(
               2.*coef*mu*(dx(u)*dx(w)+dy(v)*dy(s)+((dx(v)+dy(u))*(dx(s)+dy(w)))/2.)
              +coef*lambda*(dx(u)+dy(v))*(dx(w)+dy(s))
	)
    -int1d(Th,3)(f1*w+f2*s)
	+on(1,u=0,v=0)   
;

////////////////////////////
//  Compliance initiale   //
////////////////////////////
elasticite;
compliance = int1d(Th,3)(f1*u+f2*v);
cout<<"Initialisation. Compliance: "<<compliance<<"  Masse: "<<masse0<<endl;
fichobj << compliance << endl ;

//////////////////////////////////
//     Boucle d'optimisation    //
//////////////////////////////////

int iter;

for (iter=1;iter< niter;iter=iter+1)

{

cout <<"Iteration " <<iter <<" ----------------------------------------" <<endl;

// Calcul des composantes du tenseur des deformations
eu = dx(u);
ev = dy(v);
euv = (dx(v)+dy(u))/2;				

intermed = coef*sqrt( eu*(lambda*(eu+ev)+2*mu*eu)
			+ev*(lambda*(eu+ev)+2*mu*ev)
			+2*euv*2*mu*euv);

theta = 1/sqrt(lagrange) * intermed;
      theta = min(1,theta);
      theta = max(theta,thetamin);

/////////////////////////////////////////////////
// Détermination du multiplicateur de Lagrange //
/////////////////////////////////////////////////
masse = int2d(Th)(theta)/volume;
masse1 = masse;

// Encadrement du multiplicateur de Lagrange
if (masse1 < masse0)
{
   lagmin = lagrange ;
   while (masse < masse0)
{
      lagrange = lagrange/2 ;
      theta = 1/sqrt(lagrange) * intermed;
      theta = min(1,theta);
      theta = max(theta,thetamin);
      masse = int2d(Th)(theta)/volume ;
};
   lagmax = lagrange ;
};
//
if (masse1 > masse0)
{
   lagmax = lagrange ;
   while (masse > masse0)
{
      lagrange = 2*lagrange ;
      theta = 1/sqrt(lagrange) * intermed;
      theta = min(1,theta);
      theta = max(theta,thetamin);
      masse = int2d(Th)(theta)/volume ;
};
   lagmin = lagrange ;
};
//
if (masse1 == masse0) 
{
   lagmin = lagrange ;
   lagmax = lagrange ;
};

// Dichotomie sur le multiplicateur de Lagrange
inddico=0;
while ((abs(1.-masse/masse0)) > 0.000001)
{
   lagrange = (lagmax+lagmin)*0.5 ;
      theta = 1/sqrt(lagrange) * intermed;
      theta = min(1,theta);
      theta = max(theta,thetamin);
   masse = int2d(Th)(theta)/volume ;
   inddico=inddico+1;
   if (masse < masse0) 
      {lagmin = lagrange ;} ;
   if (masse > masse0)
      {lagmax = lagrange ;} ;
};

cout<<"Nombre d'iterations de Lagrange "<<inddico<<endl;

////////////////
// Résolution //
////////////////
if (iter <= niternp)
   {coef=theta ;} ;
if (iter > niternp)
   {coef=0.5*(1.-cos(pi*theta)) ;} ;

elasticite;

// Calcul de la compliance
compliance = int1d(Th,3)(f1*u+f2*v);
fichobj << compliance << endl ;


////////////////////////////////////////
// Affichage de la densité de matière //
////////////////////////////////////////

legende="Iteration "+iter+", Compliance "+compliance+", Masse="+masse;
plot(Th,theta,fill=1,value=true,viso=vviso,cmm=legende,wait=0);
//plot(Th,theta,fill=1,value=true,viso=vviso,cmm=legende,wait=0,ps=sauve+iter+".eps");
if (iter == niternp)
//On sauve la forme finale non penalisée
{ ofstream file(sauve+".np.bb");
	file << theta[].n << " \n";
	int j;
	for (j=0;j<theta[].n ; j++)
	file << theta[][j] << endl;  }
;
//
};
/////////////////////////////////////
// Fin de la boucle d'optimisation //
/////////////////////////////////////

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