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

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

ofstream compli("LongCantilever3dSIMP.obj") ;
string save="LongCantilever3dSIMP";		// Prefix of backup files
int niter=50;			// Number of iterations
int n=5;			// 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 complianceold;
real densityint;
real densityold;
real volume0;			// Initial volume 
real volume,volume1;		// Volume of the current shape
real del=0.10;			// Filtering Length Scale
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=-100;
func g3=0;
real penexp=3.0;			//Material Penalization Exponent
real mov=0.0;
real move=1.0;
real[int] vviso(21);
real[int] vvviso(21);
for (int i=0;i<21;i++)
{
vviso[i]=0.90+i*0.005 ;
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 material density //
/////////////////////// /////////////////////////////
real hmin=0.001;
real hmax=1.0;
real hmoy=0.40;
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.05*cos(t); y=0.05*sin(t);label=101;};

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

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

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

int[int] rup=[0,200,1,201],  rdown=[0,100,1,101];
real zmax=3.0;
real zmin=0.0;
int nlayers=60;

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;

border cutbd(t=0,1)  {x=3.0;y=-1+2*t; label=2300;};
border cutbs(t=0,1)  {x=3.0-3.0*t;y=1;  label=2300;};
border cutbg(t=0,1)  {x=0.0;y=1-2*t; label=2300;};
border cutbi(t=0,1)  {x=3.0*t;y=-1;label=2300;};
mesh XCUT=buildmesh(cutbd(m)+cutbs(nlayers)+cutbg(m)+cutbi(nlayers));

fespace VhX2(XCUT,P2);
fespace VhX0(XCUT,P0);

VhX0 hxcut;
/////////////////////////////////////////////
// Definition of the finite element spaces //
/////////////////////////////////////////////
fespace Vh0(Th3,P03d);
fespace Vh1(Th3,P13d);
fespace Vh2(Th3,P23d);


Vh2 u,v,w,p,q,r;
Vh0 h,hold,density,gradJ,gradJraw,mt;
real normgradJ;

h = h0 ;
gradJraw=0.0;
///////////////////////
// Elasticity system //
///////////////////////
problem elasticity([u,v,w],[p,q,r],solver=CG) =
    int3d(Th3)(
	      2.*h^(penexp)*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^(penexp)*lambda*(dx(u)+dy(v)+dz(w))*(dx(p)+dy(q)+dz(r))
	      )
	      -int2d(Th3,200)(g1*p+g2*q+g3*r)	           
	      +on(101,u=0,v=0,w=0)+on(100,u=0,v=0,w=0)
;

/////////////////////////////////
//  Filtering (Regularization) //
/////////////////////////////////

problem regularization(gradJ,mt) =
    int3d(Th3)(
	      -del^(2.0)*(dx(gradJ)*dx(mt)+dy(gradJ)*dy(mt)+dz(gradJ)*dz(mt))+gradJ*mt
	      )
   +int3d(Th3)(
	      -gradJraw*mt
	      );


////////////////////
// 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)  
{
if (iter>0)
{
penexp=3.0;
};
cout <<"Iteration " <<iter <<" ----------------------------------------" <<endl;

complianceold = compliance;
hold = h ;

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

density = 2*(h)^(penexp)*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)^(penexp)*lambda*(dx(u)+dy(v)+dz(w))^2;
densityint = int3d(Th3)(density);
densityold = densityint;

gradJraw= -2*(penexp)*(h)^(penexp-1)*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.)-(penexp)*(h)^(penexp-1)*lambda*(dx(u)+dy(v)+dz(w))^2;

/////////////////////////////////////
// Filtering the movement function //
/////////////////////////////////////
regularization;
//gradJ=gradJraw;
normgradJ=sqrt(int3d(Th3)(gradJ*gradJ));
mov=move/normgradJ;

h = hold-(mov*gradJ+lagrange);
//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 = hold-(mov*gradJ+lagrange);
//      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 = hold-(mov*gradJ+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 = hold-(mov*gradJ+lagrange);
//   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;

density = 2*(h)^(penexp)*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)^(penexp)*lambda*(dx(u)+dy(v)+dz(w))^2;
densityint = int3d(Th3)(density);

//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;
cout<<"Energy Density Integral: "<<densityint<<" Volume: "<<volume<<" Lagrange: "<<lagrange<<endl;


if (compliance<=complianceold)
{
  cout<<"Step Accepted"<<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);

{ 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=-0.5+i/10.; // compute yy.
  // do 3d -> 2d interpolation.
// hcut= h(x,y,yy);
hxcut=h(yy,y,x);
 caption="Design CUT z= "+yy+" , Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
 plot(hxcut,cmm=caption,fill=1,value=false,viso=vvviso,ps=save+iter+"CUT"+i+".eps",wait= 0);
}

};

if (compliance>complianceold)
{
  cout<<"Step Rejected"<<endl;
  move=move/2.;
  cout<<"new step size = "<<move<<endl;
  h = hold;
  compliance=complianceold;
  elasticity;
};
if(move<=1.0e-6)
{
iter=niter;
};
/////////////////////
// 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=-0.5+i/10.; // compute yy.
  // do 3d -> 2d interpolation.
// hcut= h(x,y,yy);
hxcut=h(yy,y,x);
 caption="Final Design CUT z= "+yy+" , Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
 plot(hxcut,cmm=caption,fill=1,value=true,viso=vvviso,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;  }