////////////////////////////////////
//                                //
//  Optimisation à la torsion de  //
//  la section d'une barre        //
//                                //
// Ecole Polytechnique, MAP 562   //
// Copyright G. Allaire, 2004     //
//                                //
////////////////////////////////////
ofstream compli("torsion.obj") ;
string sauve="torsion";		// Nom du fichier de sauvegarde
int niter=30;			// Nombre d'itérations
int n=3;			// Taille du maillage
real lagrange=1.;		// Multiplicateur de Lagrange pour le volume
real lagmin, lagmax ;		// Encadrement du multiplicateur de Lagrange
int inddico ;
real compliance;		// Compliance
real volume0;			// Volume initial
real volume,volume1;		// Volume de la forme courante
string legende;			// Legende pour les sorties graphiques
func g=1;  			// Force appliquée (couple de torsion)

real[int] vviso(21);
for (int i=0;i<21;i++)
vviso[i]=i*0.05 ;

///////////////////////////////////////////////////////////
// Inverses des modules de cisaillement des deux phases  //
///////////////////////////////////////////////////////////
real alpha=0.1;
real beta=1.0;
///////////////////////////////////////////////////////
// Initialisation de la proportion de la phase alpha //
///////////////////////////////////////////////////////
real thetamoy=0.5;
func theta0=thetamoy ;

/////////////////////////////////////////
// Définition du domaine               // 
// Conditions aux limites de Dirichlet //
/////////////////////////////////////////
mesh Th ;

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

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

fespace Vh0(Th,P0);			//Définition de l'espace P0
//fespace Vh1(Th,P2);			//Définition de l'espace P2
fespace Vh1(Th,P1);			//Définition de l'espace P1

Vh1 u,v;
Vh0 theta,densite;

///////////////////////////////////
// Proportion de la phase alpha  //
///////////////////////////////////
theta = theta0 ;
///////////////////////////
// Definition du système //
///////////////////////////
problem conduction(u,v) =
    int2d(Th)(
               alpha*beta*(dx(u)*dx(v)+dy(u)*dy(v))/(beta*theta+alpha*(1-theta))
	)
    -int2d(Th)(g*v)	           
    +on(1,u=0)
;

//////////////////////////////
// Calcul du volume initial //
//////////////////////////////
volume0=int2d(Th)(theta);

/////////////////////////
// Compliance initiale //
/////////////////////////
conduction;
compliance=int2d(Th)(g*u);
compli << compliance  << endl ;
cout<<"Initialisation. Compliance: "<<compliance<<" Volume: "<<volume0<<endl;


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

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


/////////////////////////////////////////
// Calcul de la densite d'energie      //
// Formule pour la proportion optimale //
/////////////////////////////////////////

densite = (dx(u)^2+dy(u)^2);
theta = ( sqrt(densite*alpha*beta*(beta-alpha)/lagrange) - alpha )/(beta-alpha);
theta = min(theta,1) ;
theta = max(theta,0) ;

///////////////////////////////////////////////////////////
// Calcul d'un encadrement du multiplicateur de Lagrange //
///////////////////////////////////////////////////////////
volume=int2d(Th)(theta);
volume1 = volume ;
//
if (volume1 < volume0)
{
   lagmin = lagrange ;
   while (volume < volume0)
{ 
      lagrange = lagrange/2 ;
      theta = ( sqrt(densite*alpha*beta*(beta-alpha)/lagrange) - alpha )/(beta-alpha);
      theta = min(theta,1) ;
      theta = max(theta,0) ;
      volume=int2d(Th)(theta) ; 
};
   lagmax = lagrange ;
};
//
if (volume1 > volume0)
{
   lagmax = lagrange ;
   while (volume > volume0)
{
      lagrange = 2*lagrange ;
      theta = ( sqrt(densite*alpha*beta*(beta-alpha)/lagrange) - alpha )/(beta-alpha);
      theta = min(theta,1) ;
      theta = max(theta,0) ;
      volume=int2d(Th)(theta) ;
};
   lagmin = lagrange ;
};
//
if (volume1 == volume0) 
{
   lagmin = lagrange ;
   lagmax = lagrange ;
};

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

inddico=0;
while ((abs(1.-volume/volume0)) > 0.000001)
{
   lagrange = (lagmax+lagmin)*0.5 ;
   theta = ( sqrt(densite*alpha*beta*(beta-alpha)/lagrange) - alpha )/(beta-alpha);
   theta = min(theta,1) ;
   theta = max(theta,0) ;
   volume=int2d(Th)(theta) ;
   inddico=inddico+1;
   if (volume < volume0) 
      {lagmin = lagrange ;} ;
   if (volume > volume0)
      {lagmax = lagrange ;} ;
};

//cout<<"nombre d'iterations de la dichotomie: "<<inddico<<endl;

// Résolution de l'equation
conduction;

//Calcul de la compliance
compliance=int2d(Th)(g*u);
compli << compliance  << endl ;
cout<<"Compliance: "<<compliance<<" Volume: "<<volume<<" Lagrange: "<<lagrange<<endl;

//////////////////////////////////////////////////////
// On affiche la proportion actuelle de phase alpha //
//////////////////////////////////////////////////////

legende="Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
plot(Th,theta,fill=1,value=0,viso=vviso,cmm=legende,wait=0,ps=sauve+iter+".eps"); 

};


//On sauve la forme finale
legende="Forme finale, Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
plot(Th,theta,fill=1,value=0,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);