//////////////////////////////////////
//                                  //
// Thickness or sizing optimization //
//      of a 3-d elastic plate      //
//                                  //
//        Ecole Polytechnique       //
//            ANR FF2A3             //
//   Copyright G. Allaire, A. Kelly //
//           January 2010           //
//////////////////////////////////////

load "msh3"
//load "tetgen"
load "medit"

ofstream compli("plate3d.obj") ;
string save="plate3d";		// Prefix of backup files
int niter=50;			// Number of iterations
int n=20;			// Size of the mesh
int m=20;
real lagrange=1.;		// Lagrange multiplier for the volume constraint
real lagmin, lagmax ;		// Bounds for the Lagrange multiplier
int inddico ;
real compliance;		// Compliance
real volume0;			// Initial volume 
real volume,volume1;		// Volume of the current shape
real exa=1;			// Coefficient for the shape deformation
string caption;			// Caption for the graphics
real E=100;			// Young modulus
real nu=0.3;			// Poisson coefficient (between -1 and 1/2)
func g1=0;  			// Applied forces 
func g2=0;
func g3=-100;
real[int] vviso(21);
real[int] vvviso(21);
for (int i=0;i<21;i++)
{
vviso[i]=0.80+i*0.010 ;
vvviso[i]=i*0.05;
};
//////////////////////////////////////
// Computation of Lame coefficients //
//////////////////////////////////////
real lambda=E*nu/((1.+nu)*(1.-2.*nu));
real mu=E/(2.*(1.+nu));

////////////////////////////////////////////////////
// Lower and upper bounds for the plate thickness //
/////////////////////// /////////////////////////////
real hmin=0.1;
real hmax=1.0;
real hmoy=0.25;
func h0=hmoy ;

/////////////////////// 			1:Dirichlet boundary condition
// Domain definition // 			2:Neumann boundary or free boundary condition
/////////////////////// 			3:Non-homogeneous Neumann or applied load
border a(t=0,2*pi){x=0.25*cos(t); y=0.25*sin(t);label=101;};

border bd(t=0,1)  {x=1;y=-1+2*t; label=1;};
border bs(t=0,1)  {x=1-2*t;y=1;  label=1;};
border bg(t=0,1)  {x=-1;y=1-2*t; label=1;};
border bi(t=0,1)  {x=-1+2*t;y=-1;label=1;};

border bdi(t=0,0.95)  {x=0.95;y=-0.95+2*t; label=2;};
border bsi(t=0,0.95)  {x=0.95-2*t;y=0.95;  label=2;};
border bgi(t=0,0.95)  {x=-0.95;y=0.95-2*t; label=2;};
border bii(t=0,0.95)  {x=-0.95+2*t;y=-0.95;label=2;};
 
mesh Th2=buildmesh(bd(m)+bs(m)+bg(m)+bi(m)+bdi(n)+bsi(n)+bgi(n)+bii(n)+a(n));
savemesh(Th2,save+"2D.msh");

fespace Vh22(Th2,P2);
Vh22 hcut;
Vh22 reg=region;
plot(Th2,reg,wait=1,value=true,fill=1);

///////////////////////
// Building the mesh //
///////////////////////

int[int] rup=[0,200,1,201,5,205],  rdown=[0,100,1,101,5,105];
real zmax=0.5;
real zmin=0.0;
int nlayers=30;

cout<<"mesh in"<<endl;
mesh3 Th3=buildlayers(Th2,nlayers,zbound=[zmin,zmax],reffaceup=rup,reffacelow=rdown);
//mesh3 Dh;
savemesh(Th3,save+".mesh");
cout<<"meshed"<<endl;
savemesh(Th3,save+".msh");
cout<<"mshed"<<endl;
/////////////////////////////////////////////
// Definition of the finite element spaces //
/////////////////////////////////////////////
fespace Vh0(Th3,P03d);
fespace Vh2(Th3,P23d);
fespace Vh1(Th3,P13d);

Vh2 u,v,w,p,q,r;
Vh1 h,hold,density;

h = h0 ;
///////////////////////
// Elasticity system //
//////////////////////
problem elasticity([u,v,w],[p,q,r],solver=CG) =
    int3d(Th3)(
	      2.*h*mu*(dx(u)*dx(p)+dy(v)*dy(q)+dz(w)*dz(r)+((dx(v)+dy(u))*(dx(q)+dy(p)))/2.
		       +((dx(w)+dz(u))*(dx(r)+dz(p)))/2.+((dy(w)+dz(v))*(dy(r)+dz(q)))/2.)
              +h*lambda*(dx(u)+dy(v)+dz(w))*(dx(p)+dy(q)+dz(r))
	)
	      -int2d(Th3,200)(g1*p+g2*q+g3*r)	           
	      +on(105,u=0,v=0,w=0)
;

////////////////////
// Initial volume //
////////////////////
volume0=int3d(Th3)(h);

////////////////////////
// Initial compliance //
////////////////////////
cout<<"elasticity in"<<endl;
elasticity;
cout<<"elasticity out"<<endl;
compliance=int2d(Th3,200)(g1*u+g2*v+g3*w);
compli << compliance  << endl ;
cout<<"Initialization. Compliance: "<<compliance<<" Volume: "<<volume0<<endl;

//////////////////////////////////
// Plot of the mesh deformation //
//////////////////////////////////
real coef = 0.2 ;
mesh3 Dh = movemesh3 (Th3,transfo=[x+coef*u,y+coef*v,z+coef*w]);
//plot(Dh,wait=0,ps=save+"initdef.eps");
cout<<"Application of Medit"<<endl;
medit("deformed", Dh);

////////////////////////////////
//     Optimisation loop      //
////////////////////////////////
 
int iter;

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

hold = h ;

///////////////////////////////////////////////
// Optimality condition for the thickness    //
// in terms of the density of elastic energy //
///////////////////////////////////////////////

density = 2*h*mu*(dx(u)^2+dy(v)^2+dz(w)^2+(dx(v)+dy(u))^2/2.+(dx(w)+dz(u))^2/2.+(dy(w)+dz(v))^2/2.)+h*lambda*(dx(u)+dy(v)+dz(w))^2;

h = sqrt(density*hold/lagrange) ;
h = min(h,hmax) ;
h = max(h,hmin) ;

/////////////////////////////////////////////////////
// Bracketing the value of the Lagrange multiplier //
/////////////////////////////////////////////////////
volume=int3d(Th3)(h);
volume1 = volume ;
//
if (volume1 < volume0)
{
   lagmin = lagrange ;
   while (volume < volume0)
{ 
      lagrange = lagrange/2 ;
      h = sqrt(density*hold/lagrange) ;
      h = min(h,hmax) ;
      h = max(h,hmin) ;
      volume=int3d(Th3)(h) ; 
};
   lagmax = lagrange ;
};
//
if (volume1 > volume0)
{
   lagmax = lagrange ;
   while (volume > volume0)
{
      lagrange = 2*lagrange ;
      h = sqrt(density*hold/lagrange) ;
      h = min(h,hmax) ;
      h = max(h,hmin) ;
      volume=int3d(Th3)(h) ;
};
   lagmin = lagrange ;
};
//
if (volume1 == volume0) 
{
   lagmin = lagrange ;
   lagmax = lagrange ;
};

///////////////////////////////////////////////////////
// Dichotomy on the value of the Lagrange multiplier //
///////////////////////////////////////////////////////

inddico=0;
while ((abs(1.-volume/volume0)) > 0.000001)
{
   lagrange = (lagmax+lagmin)*0.5 ;
   h = sqrt(density*hold/lagrange) ;
   h = min(h,hmax) ;
   h = max(h,hmin) ;
   volume=int3d(Th3)(h) ;
   inddico=inddico+1;
   if (volume < volume0) 
      {lagmin = lagrange ;} ;
   if (volume > volume0)
      {lagmax = lagrange ;} ;
};

//cout<<"number of iterations of the dichotomy: "<<inddico<<endl;

// Solving the elasticity system
elasticity;

//Computation of the compliance
compliance=int2d(Th3,200)(g1*u+g2*v+g3*w);
compli << compliance  << endl ;
cout<<"Compliance: "<<compliance<<" Volume: "<<volume<<" Lagrange: "<<lagrange<<endl;

//////////////////////////////////////////////
// Plot the thickness of the current design //
//////////////////////////////////////////////

caption="Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
plot(Th3,h,dim=3,fill=1,viso=vviso,cmm=caption,wait=0,ps=save+iter+".eps");

{ ofstream file(save+iter+".bb");
	file << h[].n << " \n";
	int j;
	for (j=0;j<h[].n ; j++)  
	file << h[][j] << endl;  }

//for(int i=1;i<10;i++)
//{
// real yy=i/20.; // compute yy.
//  // do 3d -> 2d interpolation.
// hcut= h(x,y,yy);
// caption="Design CUT z= "+yy+" , Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
// plot(hcut,cmm=caption,fill=1,value=true,viso=vvviso,wait= 1);
//}

//plot(Th3,h); 
//medit("Current Design", Th3, h, order=1);

/////////////////////
// End of the loop //
/////////////////////
};

//Plot the final design
caption="Final design, Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
//plot(Th,h,fill=1,value=1,viso=vviso,cmm=caption,ps=save+".eps");
medit("Sol", Th3, h, order=2);

for(int i=1;i<10;i++)
{
 real yy=i/20.; // compute yy.
  // do 3d -> 2d interpolation.
 hcut= h(x,y,yy);
 caption="Final design CUT z= "+yy+" , Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
 plot(hcut,cmm=caption,fill=1,value=true,viso=vvviso,ps=save+"FINALz"+i+".eps",wait= 1);
}


//////////////////////////////////
// Plot of the mesh deformation //
//////////////////////////////////
Dh = movemesh3 (Th3,transfo=[x+coef*u,y+coef*v,z+coef*w]);
//plot(Dh,wait=0,ps=save+"def.eps"); 
medit("deformed", Dh);


{ ofstream file(save+".bb");
	file << h[].n << " \n";
	int j;
	for (j=0;j<h[].n ; j++)  
	file << h[][j] << endl;  }