///////////////////////////////////////
//                                   //
//     Optimisation de la forme      //
// d'un radiateur de refroidissement //
//                                   //
// fonction cout = compliance =      //
// énergie de dissipation thermique  //
//                                   //
// Ecole Polytechnique, MAP 562      //
// Copyright G. Allaire, 2004        //
//                                   //
///////////////////////////////////////
ofstream objectif("radiateur.obj") ;
string sauve="radiateur";	// Nom du fichier de sauvegarde
int niter=50;			// Nombre d'itérations
int npen= 30;			// Début des itérations avec pénalisation
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;			// Légende pour les sorties graphiques
real pi = 4*atan(1) ;
func g=1;  			// Flux de chaleur appliqué

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

/////////////////////////////////
// Conduction des deux phases  //
/////////////////////////////////
real alpha=0.01;
real beta=1.0;
func theta0=(x>0.1)*(x<0.9)*(y<0.9) ;

///////////////////////////// 			1:Condition de Dirichlet
// Définition du domaine   // 			2:Condition de Neumann nulle
///////////////////////////// 			3:Flux de chaleur impose
mesh Th ;

border bd(t=0,1)   { x=1; y=t;label=2; };   // bord droit de la forme
border bg(t=1,0)   { x=0; y=t;label=2; };   // bord gauche de la forme
border bs(t=1,0)   { x=t; y=1;label=3; };   // bord supérieur de la forme
border b1(t=0,0.1)   { x=t; y=0;label=1; };   // bord inférieur de la forme
border b2(t=0.1,0.9)   { x=t; y=0;label=2; };   // bord inférieur de la forme
border b3(t=0.9,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)+b1(1*n)+b2(8*n)+b3(1*n));
plot(Th,wait=0,ps=sauve+".maillage.eps"); 
savemesh(Th,sauve+".msh");

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

Vh2 u,v;
Vh0 theta,densite, thetb;

//////////////////////////////////
// Proportions des deux phases  //
//////////////////////////////////
theta = theta0 ;
thetb = 1-theta ;
///////////////////////////
// Définition du système //
///////////////////////////
problem conduction(u,v) =
    int2d(Th)(
               (alpha*theta+beta*(1-theta))*(dx(u)*dx(v)+dy(u)*dy(v))
	)
    -int1d(Th,3)(g*v)	           
    +on(1,u=0)
;

//////////////////////////////
// Calcul du volume initial //
//////////////////////////////
volume0=int2d(Th)(1-theta0);

/////////////////////////
// Compliance initiale //
/////////////////////////
conduction;
compliance=int1d(Th,3)(g*u);
objectif << compliance  << endl ;
cout<<"Initialisation. Compliance = "<<compliance<<" Volume = "<<volume0<<endl;
thetb=1-theta ;
int iter;
iter=0;
legende="Iteration "+iter+", Compliance "+compliance+", Volume="+volume0;
plot(Th,thetb,fill=1,value=0,viso=vviso,cmm=legende,wait=0,ps=sauve+iter+".eps"); 
{ ofstream file(sauve+iter+".bb");
	file << thetb[].n << " \n";
	int j;
	for (j=0;j<thetb[].n ; j++)  
	file << thetb[][j] << endl;  }


///////////////////////////////////
//     Boucle d'optimisation     //
///////////////////////////////////
 
for (iter=1;iter< niter;iter=iter+1)  
{
cout <<"Iteration " <<iter <<" ----------------------------------------" <<endl;


/////////////////////////////////////////
// Calcul de la densité d'énergie      //
// Formule pour la proportion optimale //
/////////////////////////////////////////

densite = (((alpha*theta+beta*(1-theta)))^2)*(dx(u)^2+dy(u)^2);
theta = ( beta - sqrt(densite*(beta-alpha)/(lagrange)) )/(beta-alpha);
theta = min(theta,1) ;
theta = max(theta,0) ;
//////////////////////////
//     Pénalisation     //
//////////////////////////
if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; };

///////////////////////////////////////////////////////////
// Calcul d'un encadrement du multiplicateur de Lagrange //
///////////////////////////////////////////////////////////
volume=int2d(Th)(1-theta);
volume1 = volume ;
//
if (volume1 < volume0)
{
   lagmin = lagrange ;
   while (volume < volume0)
{ 
      lagrange = lagrange/2 ;
      theta = ( beta - sqrt(densite*(beta-alpha)/(lagrange)) )/(beta-alpha);
      theta = min(theta,1) ;
      theta = max(theta,0) ;
//////////////////////////
//     Pénalisation     //
//////////////////////////
if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; };
      volume=int2d(Th)(1-theta) ; 
};
   lagmax = lagrange ;
};
//
if (volume1 > volume0)
{
   lagmax = lagrange ;
   while (volume > volume0)
{
      lagrange = 2*lagrange ;
      theta = ( beta - sqrt(densite*(beta-alpha)/(lagrange)) )/(beta-alpha);
      theta = min(theta,1) ;
      theta = max(theta,0) ;
//////////////////////////
//     Pénalisation     //
//////////////////////////
if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; };
      volume=int2d(Th)(1-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 = ( beta - sqrt(densite*(beta-alpha)/(lagrange)) )/(beta-alpha);
   theta = min(theta,1) ;
   theta = max(theta,0) ;
//////////////////////////
//     Pénalisation     //
//////////////////////////
if (iter > npen) { theta = ( 1 - cos(pi*theta) )/2 ; };
   volume=int2d(Th)(1-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'équation
conduction;

//Calcul de la compliance
compliance=int1d(Th,3)(g*u);
objectif << compliance  << endl ;
cout<<"Compliance: "<<compliance<<" Volume: "<<volume<<" Lagrange: "<<lagrange<<endl;

////////////////////////////////////////////////////
// On affiche la conductivité du radiateur actuel //
////////////////////////////////////////////////////
thetb=1-theta ;
legende="Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
plot(Th,thetb,fill=1,value=0,viso=vviso,cmm=legende,wait=0,ps=sauve+iter+".eps"); 
{ ofstream file(sauve+iter+".bb");
	file << thetb[].n << " \n";
	int j;
	for (j=0;j<thetb[].n ; j++)  
	file << thetb[][j] << endl;  }

};


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