////////////////////////  PARAMETERS ////////////////////////////
verbosity=0;

macro g [0.,-1.]							// Applied surface load 
real Y=15;								// Young Modulus
real nu=0.35;								// Poisson Ratio
real lambda=Y*nu/((1.+nu)*(1.-2.*nu));					// Lame Moduli
real mu=Y/(2.*(1.+nu));
int Niter=50;								// Number of iterations
real penVolume=0.05;							// Penalization of the volume
real GeoPrec=0.1;							// scale of the geometric mesh
real step=8.;								// initial descent step			
string name="cantilever";
bool poursuit=false,adapt=false;					// Restart at previous computed state ? / Adapt mesh wrt the displacement ?

///////////////////////////////////////////
////  Definition  of the initial mesh  ////
///////////////////////////////////////////
	int Free,Neumann,Dirichlet,Dirichlet0; 				// Labels for optimized boundaries, Neumann and Dirichlet conditions
		{Free=1;Neumann=2;Dirichlet=3;Dirichlet0=7;};
	int FreezeX,FreezeY,FreezeXY; 					// Freeze degree of freedom of theta
		{FreezeX=4;FreezeY=5;FreezeXY=6;}
	mesh Sh;							// The Total mesh
	if(poursuit) {Sh=readmesh(name+"-Sh.msh");} else
	{
		// The definition of the initial mesh could be modified as will //
		// -- a cantilever --
		real D=1,N=0.1,L=4,D0=0.5,D00=0.1,deltax=0.5,deltay=1.;
		int n=5;
		border Left0(t=D,D-D00){x=0;y=t;label=Dirichlet0;};
		border Left1(t=D-D00,D-D0){x=0;y=t;label=Dirichlet;};
		border Left2(t=D-D0,-D+D0){x=0;y=t;label=Free;};
		border Left3(t=-D+D0,-D+D00){x=0;y=t;label=Dirichlet;};
		border Left4(t=-D+D00,-D){x=0;y=t;label=Dirichlet0;};
		border Bottom(t=0,L){x=t;y=t*(D-N)/L-D;label=Free;};
		border Right(t=-N,N){x=L;y=t;label=Neumann;};
		border Top(t=L,0){x=t;y=t*(N-D)/L+D;label=Free;};
		
		real r=0.1;
		border Hole1(t=2.*pi,0.){x=r*cos(t)+L/2.;y=r*sin(t);label=Free;}
				
		// -- BIG BOX containing the shape --//
		border DLeft(t=D+deltay,-D-deltay){x=-deltax;y=t;label=FreezeX;};
		border DBottom(t=-deltax,L+deltax){x=t;y=-D-deltay;label=FreezeY;};
		border DRight(t=-D-deltay,D+deltay){x=L+deltax;y=t;label=FreezeX;};
		border DTop(t=L+deltax,-deltax){x=t;y=D+deltay;label=FreezeY;};
		
		Sh=buildmesh(
			  DLeft(2*D*n)+DBottom(L*n)+DRight(2*D*n)+DTop(L*n)
			+Left0(D0*n)+Left1(D0*n)+Left2(2*D*n)+Left3(D0*n)+Left4(D0*n)+Bottom(L*n)+Right(2*D*n)+Top(L*n)
			+Hole1(-10));
	}
	
	mesh Th;					// The mesh of the shape
	int ShapeLabel;					// Region label of the shape
	fespace Zh(Sh,P0); Zh reg=region;		// Region labels
	if(poursuit)
		{
		
			ifstream Label(name+"-label"); Label>>ShapeLabel;
		}
	else
	{
		real px=0.1; real py=0;		// a point of the shape
		ShapeLabel=reg(px,py);
	}
	Th=trunc(Sh,reg==ShapeLabel);
	plot(Th,wait=1);
	
////////////////////////////////////////
//     end of initial mesh definition //
////////////////////////////////////////

////////////////////// END OF PARAMETERS LIST ///////////////////
fespace Xh(Sh,[P1,P1]);							// finite element space for the convection field
fespace Wh(Th,[P2,P2]);							// finite element space for the displacement 
int quadraticOrder=3;
// Some macro
	macro u [u1,u2]							// displacement
	macro v [v1,v2] 						// test function for elasticity
	macro e(u) [dx(u[0]),(dx(u[1])+dy(u[0]))/sqrt(2.),dy(u[1])] 	// deformation tensor
	macro div(u) (dx(u[0])+dy(u[1]))				// divergence operator
	macro theta [theta1,theta2]					// descent direction
	macro ptheta [ptheta1,ptheta2]					// previous descent direction
	macro dtheta [dtheta1,dtheta2]					// test function fot descent direction
	macro Theta [Theta1,Theta2]					// interpolated descent direction

Wh u,v;
problem elasticity(u,v)=
	 int2d(Th,qforder=quadraticOrder)(2.*mu*(e(u)'*e(v))+lambda*div(u)*div(v))
	-int1d(Th,Neumann)(g'*v)
	+on(Dirichlet,u1=0,u2=0);

Xh theta,dtheta;

problem Computetheta(theta,dtheta)=
	int2d(Sh,qft=qf1pT)(2.*mu*(e(theta)'*e(dtheta))+lambda*div(theta)*div(dtheta))
	+int1d(Sh,Free)(jump(reg)*(2.*mu*(e(u)'*e(u))+lambda*div(u)^2-penVolume)*([N.x,N.y]'*dtheta))	
	+on(Dirichlet0,Neumann,theta1=0,theta2=0)
	+on(Dirichlet,theta1=0)
	+on(FreezeX,theta1=0)
	+on(FreezeY,theta2=0);
	

varf dJ(theta,dtheta)=	// used for step adaptation ...
	int1d(Sh,Free)(jump(reg)*(2.*mu*(e(u)'*e(u))+lambda*div(u)^2-penVolume)*([N.x,N.y]'*dtheta))	
	+on(Dirichlet0,Neumann,theta1=0,theta2=0)
	+on(Dirichlet,theta1=0)
	+on(FreezeX,theta1=0)
	+on(FreezeY,theta2=0);

real J;

for(int iter=0;iter<Niter;iter++)
{// Optimization loop	
		// --------------- Solve Elasticity problem and evaluation of the cost function -------------//
		cout<<"mesh adaptation"<<endl;
		Sh=adaptmesh(Sh,GeoPrec,IsMetric=1); 
		cout<<"solve elasticity"<<endl;
		reg=(region==ShapeLabel);Th=trunc(Sh,(reg==1));
		elasticity;
		// adaptation of the mesh ?
		if(adapt) {Th=adaptmesh(Th,u,err=5.e-3,hmax=GeoPrec/3.); elasticity;}
		J=int1d(Th,Neumann)(g'*u)+penVolume*int2d(Th)(1.);
		cout<<"Iteration "<<iter<<"  --- J="<<J<<endl;
		
		// ---------------                  Compute a descent direction                -------------//
		Computetheta;
		real normetheta=sqrt(int2d(Th)(theta'*theta));
		plot(theta,cmm=iter+"; J="+J+"; step="+step,coef=1./normetheta);
			
			
		
		// ------------           Move the mesh and optimize "time" step                -------------//
		{
			real ps0,ps;mesh pSh=Sh;Xh ptheta=theta; 			// Used for time step adaptation
			{
				real[int] DJ(Xh.ndof); DJ=dJ(0,Xh);
				ps0=((DJ)'*theta1[]);
			};
			bool OK=false;
			while(!OK)	// OK=false if "oscillations" are detected... in this case, we retry with a smaller step
			{
				// Move Sh along Theta during time dt=step. 
				{
					Xh [Id1,Id2]=[x,y];  
					Xh [X,Y];
					Sh=pSh;theta=ptheta;
					bool moveOK=false;
					while(!moveOK)
					{
						// Try to move the mesh along theta during time step
						cout<<"convection"<<endl;
						[X,Y]=[convect(theta,step,Id1),convect(theta,step,Id2)];	// Can bug here for hight velocity
						cout<<"convection DONE"<<endl;
						real minT0= checkmovemesh(pSh,[x,y]);	// min triangle area
						if (checkmovemesh(pSh,[X,Y])>minT0/100) {
							Sh=movemesh(pSh,[X,Y]);moveOK=true;
						} else {
								step/=2;
								cout<< "                            Triangle reversed = step/2 ="<<step<<endl;
							}
						// ... if Triangle are reversed, we try again with a smaller time step
					}
				}	
				// Interpolation of theta on the new mesh
				{
					Xh [inter1,inter2]=[0,0]; inter1[]=theta1[]; theta=[0,0]; theta1[]=inter1[];
				}
			
				// Check if oscillations... if so, step is decreased, otherwise OK is set to true
				{
					// solve elasticity
					reg=(region==ShapeLabel);Th=trunc(Sh,(reg==1));elasticity;
					if(adapt){Th=adaptmesh(Th,u,err=5.e-3,hmax=GeoPrec/3.); elasticity;}
					//  oscillations ?
					{
						real[int] DJ(Xh.ndof); DJ=dJ(0,Xh);
						theta=theta;ps=((DJ)'*theta1[]);
						cout<<"step adaptation"<<endl;
						if (ps>0) {
							step/=2.;
							cout<<"GO BACK !!!!  step="<<step<<endl;
							} 
						else OK=true;
					}
				}
			}
					
			// Everything fine... let's try to increase time step if possible
			cout<<"step="<<step<<endl;
			real coef=ps0/(ps0-ps); coef=min(coef,2.);
			if (coef>1.) step*=(1.+coef)/2.;
		} // mesh moved and time step optimized... next iteration !
}// Optimization loop

// save result
savemesh(Sh,name+"-Sh.msh");
{ofstream Label(name+"-label"); Label<<ShapeLabel;}
plot(Th,cmm="Final Shape J="+J);
cout<<"J="<<J<<endl;