///////////////////////////////////////////////////////
//  Level set method for the geometric optimization  //
//        of a minimal compliance cantilever         //
//                                                   //
//         Ecole Polytechnique, MAP 562              //
//    Copyright G. Allaire, C. Dousset 2007          //
///////////////////////////////////////////////////////
// La fonction ligne de niveaux est appelée phi. La 
// forme est caractérisée par phi<0. La fonction phi
// est convectée par un schéma caractéristique implicite
// (cf. J. Strain, JCP 151, 1999).
// The level set function is denoted by phi. The shape 
// is characterized by phi<0. The function phi is advected 
// by a caracteristic implicit scheme (cf. J. Strain, 
// JCP 151, 1999).
///////////////////////////////////////////////////////
int n=2;	//mesh size parameter
real mu=8;	//Lame coefficient
real lambda=1;	//Lame coefficient 
real g1=0;	//horizontal component of the applied force on the boundary
real g2=-5;	//vertical component of the applied force on the boundary
real lag=0.5;	//Lagrange multiplier for the weight
real T=.1;	//descent step or time step for the advection
real TT=0.01;	//time step for the re-initialization
real weak=0.0001;	//multiplicative coefficient for the "weak" material mimicking void
real eps1=0.001;	//small parameter avoiding division by zero in the normal computation
real eps2=0.001;	//regularization parameter for the velocity
real tol=0.0001; 	//tolerance for the increase of the objective function
real compliance;
real weight;
real objective;
real objectiveinitial;
real compliancem;
real weightm;
real objectivem;
int kx=4;	//number of holes in the horizontal direction of the initial shape
int ky=3;	//number of holes in the vertical direction of the initial shape
int N=30;	//number of iterations for the gradient algorithm
int Ninit=5;	//number of iterations for the re-initialization
int iter;
int iterinit;
string caption;
string save="cantilever-ls";
ofstream object("cantilever-ls.obj");

//Definition of the boundary and of the mesh

// right boundary of the mesh
border bd1(t=-0.5,-0.05)	{ x=1;  y=t; label=1; };
border bd2(t=-0.05,0.05)	{ x=1;  y=t; label=3; };
border bd3(t=0.05,0.5)		{ x=1;  y=t; label=1; };
//left  boundary of the mesh
border bg(t=0.5,-0.5)		{ x=-1; y=t; label=2; };
// upper boundary of the mesh
border bs(t=1,-1)		{ x=t;  y=0.5; label=1; };
// lower boundary of the mesh
border bi(t=-1,1)		{ x=t;  y=-0.5; label=1; };
// label meanings: 1=Neumann, 3=applied force, 2=Dirichlet

mesh Th= buildmesh (bd1(9*n)+bd2(2*n)+bd3(9*n)+ bg(20*n)+bs(40*n)+bi(40*n));

//Definition of finite element spaces

fespace Vh(Th,[P1,P1]);
fespace Vh1(Th,P1);

Vh [u,v],[w,s],[um,vm];
Vh1 N1, N2, phi, V, phim, M1, M2, S, phiinit, nabla, Xapprox, X, 
Vreg, vv, Xapproxm, Xm;

// Definition of the initial shape (perforated by holes according to kx and ky)

func phi0=-0.1+sin(pi*kx*x)*sin(pi*ky*(y-0.5));
phi=phi0;
plot(phi,wait=0,fill=1,value=1,ps=save+".init.eps");

// Re-initialization of phi such that it becomes the signed distance  
// function to the shape boundary. Implicit linearized solving of :
// (d phi)/(d t) + sign(phi) ( |gradphi)| - 1 ) = 0

Vh1 h1=hTriangle;
real h=h1[].max; //maximal size of a mesh triangle
cout<<"maximal size of a mesh triangle h="<<h<<endl;

for (iterinit=1;iterinit< 20;iterinit=iterinit+1)
{
	nabla=(dx(phi))^2+(dy(phi))^2;
	S=phi/(sqrt(phi^2+h*h*nabla)); //approximation of sign(phi)
	M1=dx(phi)/(sqrt(nabla+eps1^2)); //approximation of the normal
	M2=dy(phi)/(sqrt(nabla+eps1^2)); //approximation of the normal
	phiinit=convect([-S*M1,-S*M2],TT,phi);
	phiinit=phiinit+TT*S;
	phi=phiinit;
};
	
//plot(phi,wait=1,fill=1,value=1,ps=save+".init-re.eps");

// Computation of the elastic displacement

X= (phi<0); //characteristic function of the shape
Xapprox = X + weak*(1-X) ; //approximation of the shape characteristic function

solve elasticity([u,v],[w,s]) = int2d(Th)(    
2*mu*Xapprox*( dx(u)*dx(w)+dy(v)*dy(s) + (dx(v)+dy(u))*(dx(s)+dy(w))/2 )   
+ lambda*Xapprox*(dx(u)+dy(v))*(dx(w)+dy(s))   )
- int1d(Th,3)(g1*w+g2*s)
+ on(2,u=0,v=0);

// Computation of the objective function (compliance)

compliance=int1d(Th,3)(g1*u+g2*v);
weight=int2d(Th)(X);
objective=compliance+lag*weight;
objectiveinitial=objective;
object << objective   << endl ;
cout<<"initial objective="<<objectiveinitial<<endl;


// Computation of the shape gradient or advection velocity
V=2*mu*Xapprox*((dx(u))^2+(dy(v))^2+((dy(u)+dx(v))^2)/2)+lambda*Xapprox*(dx(u)+dy(v))^2;
V=V-lag*X;

// Definition of the velocity regularization
problem smoothing (Vreg,vv)=
int2d(Th)( eps2*( dx(Vreg)*dx(vv)+dy(Vreg)*dy(vv) ) + Vreg*vv ) - int2d(Th)(V*vv);
// Regularization of the velocity
smoothing ;
V=Vreg;	

// Computation of the normal
nabla=(dx(phi))^2+(dy(phi))^2;
N1=dx(phi)/(sqrt(nabla+eps1^2));
N2=dy(phi)/(sqrt(nabla+eps1^2));

// Plot of the shape 
plot(X,fill=1);

///////////////////////////////////////////////
// Loop of the gradient algorithm iterations //
///////////////////////////////////////////////
for (iter=1;iter< N;iter=iter+1)
{
	cout << "iteration= "<<iter  << endl;
	cout << "descent step T= "<<T  << endl;
	cout<<"previous objective="<<objective<<endl;

// Shape convection
	phim=convect([-V*N1,-V*N2],T,phi);
// This result is provisional as far as it has not been checked that 
// the objective function decreases. 

// Re-initialization of phim (only if the descent step is not too small).

if (T>0.001)
{
for (iterinit=1;iterinit< Ninit;iterinit=iterinit+1)
{
	nabla=(dx(phim))^2+(dy(phim))^2;
	S=phim/(sqrt(phim^2+h*h*nabla));
	M1=dx(phim)/(sqrt(nabla+eps1^2));
	M2=dy(phim)/(sqrt(nabla+eps1^2));
	phiinit=convect([-S*M1,-S*M2],TT,phim);
	phiinit=phiinit+TT*S;
	phim=phiinit;
};
};

// Computation of the provisional elastic displacement

	Xm= (phim<0);
	Xapproxm = Xm + weak*(1-Xm);

solve elasticitym([um,vm],[w,s]) = int2d(Th)(    
2*mu*Xapproxm*( dx(um)*dx(w)+dy(vm)*dy(s) + (dx(vm)+dy(um))*(dx(s)+dy(w))/2 )   
+ lambda*Xapproxm*(dx(um)+dy(vm))*(dx(w)+dy(s))   )
- int1d(Th,3)(g1*w+g2*s)
+ on(2,um=0,vm=0);

// Computation of the provisional objective function

	compliancem=int1d(Th,3)(g1*um+g2*vm);
	weightm=int2d(Th)(Xm);
	objectivem=compliancem+lag*weightm;
	cout <<"provisional objective="<<objectivem<<endl;

// Test to accept or not the descent step according 
// to the decrease or increase of the objective function.

// Accept the descent step
	if (objectivem<=objective*(1+tol))
	{
		compliance=compliancem;
		weight=weightm;
		objective=objectivem;
		object << objective   << endl ;
		phi=phim;
		[u,v]=[um,vm];
		X=Xm;
		Xapprox=Xapproxm;
		plot(X,fill=1);
		if (T<0.05) T=T*1.5;

		V=mu*Xapprox*(2*(dx(u))^2+2*(dy(v))^2+(dy(u)+dx(v))^2)+lambda*Xapprox*(dx(u)+dy(v))^2;
		V=V-lag*X; //Computation of the shape gradient or advection velocity
		smoothing ; //Regularization of the velocity
		V=Vreg;

		nabla=(dx(phi))^2+(dy(phi))^2; //Computation of the normal
		N1=dx(phi)/(sqrt(nabla+eps1^2));
		N2=dy(phi)/(sqrt(nabla+eps1^2));

		cout <<"descent step accepted"<<endl;
	};

// Reject the descent step
	if (objectivem>objective*(1+tol))
	{
		T=T/2;
		cout <<"descent step rejected"<<endl;
	
	};

// End of the iteration loop
}; 


caption="Final shape, N= "+N+", Kx= "+kx+", Ky= "+ky+", objective= "+objective+",compliance= "+compliance+", Lagrange multiplier="+lag;
plot(Th,X,fill=1,cmm=caption,ps=save+".eps");