///////////////////////////////////////////
// MAP 562 (G. Allaire, T. Wick)	 //
// Winter 2016/2017			 //
// Feb 7, 2017				 //
//					 //	
// Modified from the previous script	 //
// by G. Allaire and O. Pantz 		 //
//					 //
// TP 6 				 //
///////////////////////////////////////////



////////////////////////////////////////
// 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;							// Penalisation of the volume
real GeoPrec=0.1;							// scale of the geometric mesh
real step=8.;								// initial descent step			
string name="cantilever";
bool restart=false,adapt=false;						// restart from previous computations / Adapt mesh or not
macro xh [P1,P1]							// finite element space for convection field
macro vh [P2,P2]							// finite element space for the displacement


/////////////////////////////////////////////
////  Definition of two initial meshes   ////
////  Sh = D        (big box) 	 	 ////
////  Th = \Omega_0 (true object)        ////
/////////////////////////////////////////////
int Free,Neumann,Dirichlet,Dirichlet0; 		// Labels for optimized boundaries, Neumann and Dirichlet conditions
int FreezeX,FreezeY,FreezeXY; 			// Freeze degree of freedom of theta
mesh Sh;					// The Full mesh

// Definition of the full mesh
{
	{Free=1;Neumann=2;Dirichlet=3;Dirichlet0=7;};
	{FreezeX=4;FreezeY=5;FreezeXY=6;}

	if(restart) {Sh=readmesh(name+"-Sh.msh");} else
	{
	
		// The definition of the initial mesh could be modified as well //
		// -- a cantilever occupying the domain \Omega_0 --

			real D=1,N=0.1,L=4,D0=0.5,D00=0.1,deltax=0.5,deltay=1.;
			int n=5;

			// On the small boundaries Left0 and Left4 we fix
			// everything to avoid singularities in the solution
			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 denoted with D 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));
	}
}



// Definition of the truncated mesh of the true object
mesh Th;					// The mesh of the object
int ShapeLabel;					// Region label of the shape
fespace Zh(Sh,P0); 
Zh reg=region;					// Region labels

// Definition of the shape and the inner label
{
	if(restart)
	  {
	    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);

}


// Visualization of the two initial meshes
plot(Th,wait=1,cmm="initial mesh of the shape");
plot(Sh,wait=1,cmm="initial mesh of the domain");



////////////////////////////////////////
// Finite element spaces
////////////////////////////////////////

// Define finite element space
fespace Vh(Th,vh), Xh(Sh,xh);		// finite element spaces for the displacements and the convection field
int quadraticOrder=3;



////////////////////////////////////////
// Define macros for short hand notation
////////////////////////////////////////
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
macro Id [Id1,Id2]		// The identity map
macro Convect(theta,step,u) [convect(theta,step,u[0]),convect(theta,step,u[1])] // convection of a vector field




////////////////////////////////////////
// Defining the underlying equations
////////////////////////////////////////


// Define elasticity (state equation)
Vh u,v;					// displacement
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);

// For the theta calculation, 
// we use the FreeFem function to define the jump 
// across an element edge.
// qft assignes explicitely a Gauss quadrature formula.
// Here, of order 2.

Xh theta,dtheta;			// convection field
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);
	

// Define the derivative required to calculate the new mesh
varf dJ(theta,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);


////////////////////////////////////////
// Optimization loop
////////////////////////////////////////

// As usual J is the cost function
real J;				


for(int iter=0;iter<Niter;iter++)
{	
		// --------------- 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));

		// Solve the state equation
		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 the convection field (in which direction the
		// shape/nodes are moved
		plot(theta,cmm=iter+"; J="+J+"; step="+step,coef=1./normetheta);
			
			
		
		// ------------           Move the mesh and adapt the "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 Id=[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,Id);	// 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 for 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;}
					//  Check for 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 !
};


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