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


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

ofstream compli("pylon3dSIMP.obj") ;
string save="pylon3dSIMP";		// Prefix of backup files
int niter=200;			// Number of iterations
int n=40;			// 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=0;
func g3=-100;
real penexp=5.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.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 pylon thickness //
/////////////////////// /////////////////////////////
real hmin=0.001;
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

real x0,x1;
real y0,y1;

border abd(t=0,1)  {x=1;y=-0.5+t; label=1;};
border abt(t=0,1)  {x=1-2*t;y=0.5;  label=1;};
border abg(t=0,1)  {x=-1;y=0.5-t; label=1;};
border abi(t=0,1)  {x=-1+2*t;y=-0.5;label=1;};

border bbd(t=0,1)  {x=0.5;y=-1.5+3*t; label=2;};
border bbt(t=0,1)  {x=0.5-t;y=1.5;  label=2;};
border bbg(t=0,1)  {x=-0.5;y=1.5-3*t; label=2;};
border bbi(t=0,1)  {x=-0.5+t;y=-1.5;label=2;};

border cbd(t=0,1)  {x=1;y=-1.5+3*t; label=1;};
border cbt(t=0,1)  {x=1-2*t;y=1.5;  label=1;};
border cbg(t=0,1)  {x=-1;y=1.5-3*t; label=1;};
border cbi(t=0,1)  {x=-1+2*t;y=-1.5;label=1;};


border dbd(t=0,1)  {x=2;y=-1+2*t; label=1;};
border dbt(t=0,1)  {x=2-4*t;y=1;  label=1;};
border dbg(t=0,1)  {x=-2;y=1-2*t; label=1;};
border dbi(t=0,1)  {x=-2+4*t;y=-1;label=1;};

border ebd(t=0,1)  {x=1;y=-1+2*t; label=1;};
border ebt(t=0,1)  {x=1-2*t;y=1;  label=1;};
border ebg(t=0,1)  {x=-1;y=1-2*t; label=1;};
border ebi(t=0,1)  {x=-1+2*t;y=-1;label=1;};

func XX=x;
func YY=y;
func ZZ=z; 

x0=-1.0;x1=1.0;
y0=-0.5;y1=0.5;
n=20;m=10;
mesh postb=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

x0=-1.0;x1=1.0;
y0=-0.5;y1=0.5;
n=20;m=10;
mesh posth=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

x0=-0.5;x1=0.5;
y0=-1.5;y1=1.5;
n=10;m=30;
mesh postgd=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

x0=-1.0;x1=1.0;
y0=-1.5;y1=1.5;
n=20;m=30;
mesh postdd=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

x0=-2.0;x1=2.0;
y0=-1.0;y1=1.0;
n=40;m=20;
mesh crossb=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

x0=-1.0;x1=1.0;
y0=-1.0;y1=1.0;
n=20;m=20;
mesh crossgd=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

x0=-2.0;x1=2.0;
y0=-1.0;y1=1.0;
n=40;m=20;
mesh crossddh=square(n,m,[x0+(x1-x0)*x,y0+(y1-y0)*y]);

int[int] refpb=[0,111];
int[int] refph=[0,112];
int[int] refpdev=[0,121];
int[int] refpder=[0,122];
int[int] refpg=[0,131];
int[int] refpd=[0,132];

int[int] refcb=[0,211];
int[int] refch=[0,212];
int[int] refcdev=[0,221];
int[int] refcder=[0,222];
int[int] refcg=[0,231];
int[int] refcd=[0,232];

mesh3 Th3pb=movemesh23(postb,transfo=[XX,YY,0.0],refface=refpb,orientation=-1);
mesh3 Th3pd=movemesh23(postgd,transfo=[1.0,XX,YY+1.5],refface=refpd,orientation=1);
mesh3 Th3pg=movemesh23(postgd,transfo=[-1.0,XX,YY+1.5],refface=refpg,orientation=-1);
mesh3 Th3pdev=movemesh23(postdd,transfo=[XX,0.5,YY+1.5],refface=refpdev,orientation=-1);
mesh3 Th3pder=movemesh23(postdd,transfo=[XX,-0.5,YY+1.5],refface=refpder,orientation=1);
mesh3 Th3ph=movemesh23(postb,transfo=[XX,YY,3.0],refface=refph,orientation=1);

mesh3 Th3cb=movemesh23(crossb,transfo=[XX,YY,3.0],refface=refcb,orientation=-1);
mesh3 Th3ch=movemesh23(crossddh,transfo=[XX,YY,5.0],refface=refch,orientation=1);
mesh3 Th3cd=movemesh23(crossgd,transfo=[2.0,XX,YY+4.0],refface=refcd,orientation=1);
mesh3 Th3cg=movemesh23(crossgd,transfo=[-2.0,XX,YY+4.0],refface=refcg,orientation=-1);
mesh3 Th3cdev=movemesh23(crossddh,transfo=[XX,1.0,YY+4.0],refface=refcdev,orientation=-1);
mesh3 Th3cder=movemesh23(crossddh,transfo=[XX,-1.0,YY+4.0],refface=refcder,orientation=1);

///////////////////////
// Building the mesh //
///////////////////////
mesh3 Th33=Th3pb+Th3pd+Th3pg+Th3pdev+Th3pder+Th3cb+Th3ch+Th3cd+Th3cg+Th3cdev+Th3cder;
medit("pylon",Th33);

mesh3 Th333=Th3pb+Th3ph+Th3pd+Th3pg+Th3pdev+Th3pder;

real[int] domain=[0.,0.,1.,66,0.001];
mesh3 Th3=tetg(Th33,switch="paAAQYVVV",nbofregions=1,regionlist=domain);

medit("pylontets",Th3);
savemesh(Th3,save+".mesh");

cout<<"========================================"<<endl;
cout<<"meshed"<<endl;
cout<<"========================================"<<endl;

border prof1(t=0,1)  {x=-1.0+2.0*t;y=0.0;};
border prof2(t=0,1)  {x=1.0;y=0.0+3.0*t;};
border prof3(t=0,1)  {x=1.0+t;y=3.0;};
border prof4(t=0,1)  {x=2.0;y=3.0+2.0*t;};
border prof5(t=0,1)  {x=2.0-4.0*t;y=5.0;};
border prof6(t=0,1)  {x=-2.0;y=5.0-2.0*t;};
border prof7(t=0,1)  {x=-2.0+t;y=3.0;};
border prof8(t=0,1)  {x=-1.0;y=3.0-3.0*t;};

cout<<"========================================"<<endl;
cout<<"profile borders"<<endl;
cout<<"========================================"<<endl;

mesh hprofile=buildmesh(prof1(n)+prof2(n)+prof3(n)+prof4(n)+prof5(n)+prof6(n)+prof7(n)+prof8(n));
//plot(hprofile);
/////////////////////////////////////////////
// Definition of the finite element spaces //
/////////////////////////////////////////////
fespace Vh0(Th3,P03d);
fespace Vh1(Th3,P13d);
fespace Vh2(Th3,P23d);

cout<<"========================================"<<endl;
cout<<"3Dspaces"<<endl;
cout<<"========================================"<<endl;

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

Vh2 loadcore=((y)^(2.0)+(z-4.0)^(2.0)-0.01<=0.0);

fespace Vh02(hprofile,P0);
Vh02 hcut;
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,231)(loadcore*(g1*p+g2*q+g3*r))-int2d(Th3,232)(loadcore*(g1*p+g2*q+g3*r))	           
	      +on(111,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<<"========================================"<<endl;
cout<<"elasticity in"<<endl;
cout<<"========================================"<<endl;

elasticity;
cout<<"========================================"<<endl;
cout<<"elasticity out"<<endl;
cout<<"========================================"<<endl;

compliance=int2d(Th3,231)(loadcore*(g1*u+g2*v+g3*w))+int2d(Th3,232)(loadcore*(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;

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 = 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 = 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 = 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;
cout<<"volume =  "<<volume<<endl;
cout<<"volume0 =  "<<volume0<<endl;
real shit=abs(1.0);
while ((abs(1.-volume/volume0)) > 0.000001)
//while (1.0 > 0.000001)
{
   lagrange = (lagmax+lagmin)*0.5 ;
   h = hold-(mov*gradJ+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,231)(loadcore*(g1*u+g2*v+g3*w))+int2d(Th3,232)(loadcore*(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=vvviso,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=-1.0+i/5.; // compute yy.
  // do 3d -> 2d interpolation.
 hcut= h(x,yy,y);
 caption="Design CUT z= "+yy+" , Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
 plot(hcut,cmm=caption,fill=1,value=true,viso=vvviso,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=-1.0+i/5.; // compute yy.
  // do 3d -> 2d interpolation.
 hcut= h(x,yy,y);
 caption="Final Design CUT z= "+yy+" , Iteration "+iter+", Compliance "+compliance+", Volume="+volume;
 plot(hcut,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;  }