// MAP 562 (G. Allaire, T. Wick)
// Winter 2016/2017
// Jan 12, 2017
//
// TD 2
// Exercise 2

// Task 1: Run the program
// Task 2: Understand how the derivative dJ is obtained (recall the 
//   theoretical exercise)
// Task 3: Change the boundary to only Dirichlet conditions
// Task 4: Change the step-size -> How many (gradient) iterations to you observe ?
// Task 5: Change the lower tolerance to stop the gradient algorithm
 
int N=50;	// Nb of edges on the boundary

// Define volume and traction forces, respectively
func f=1;	
func g=-1;	

// Define the p-value (p > 2)
int p=4;

// Define the tolerance at which we want to stop
// the calculation
real lowertolerance =1.e-4;

// Define the polynomial degree
macro vh() P1 //

// Define the mesh and boundaries
mesh Th;
int GammaN,GammaD;
{
	GammaN=1;GammaD=2;
	int[int] labels=[GammaN,GammaD,GammaD,GammaN];
	Th=square(N,N,label=labels);
}

fespace Vh(Th,vh);

// Define the trial and test functions
// Here the trial function is the solution du
// of a projection.
// And phi is the test function
// The solution u^{k} is denoted by u
Vh du,phi,u;		

// Initialize the initial guess u^0 as zero
u = 0;

macro grad(u) [dx(u),dy(u)] //	  

// Implement the functional given in the exercise
macro J(u) int2d(Th)(grad(u)'*grad(u)/2.+(grad(u)'*grad(u))^(p/2)/p) 
         - int2d(Th)(f*u)
	 - int1d(Th,GammaN)(g*u)  //  

// Define the problem
// Here the variational form is obtained by calculating the Gateaux-derivative
// of the functional J
// Rather than evaluating directly the derivative at the previous 
// step u^k, we compute the projection (du,phi) = J'(u)(phi)
// in an appropriate scalar product. Here we work in the H^1 space,
// thus one possibility is \int \nabla du \nabla phi, as shown 
// just below
problem Jlin(du,phi) = -int2d(Th)(grad(du)'*grad(phi))
		       +int2d(Th)((1+sqrt(grad(u)'*grad(u))^(p-2))*(grad(u)'*grad(phi)))
		       -int2d(Th)(f*phi)
		       -int1d(Th,GammaN)(g*phi)
		       +on(GammaD,du=0);


// Solution loop
real funcnorm = 1e+10;
real stepsize = 0.1;

// Counter of the number of gradient steps
int i=0;

// Final loop
while(funcnorm > lowertolerance) 
{
	// Solve the projection problem for obtaining the update du
	Jlin;

	// Compute new iterate u^{k+1}
	u = u - stepsize * du;

	// Compute the norm of the updates, which is 
	// used as criterion to stop the final loop.
	// This criterion is reasonable since the updates du should 
	// tend to zero while approaching the 
	// optimal solution u^*
	funcnorm = sqrt(int2d(Th)(du*du));

	// Evaluate the sought functional value
	real residual = J(u);

	// Output on the terminal and visualization 
	cout<<"#"<<i<<" J="<<residual << "  " << funcnorm<<endl;
	plot(u,cmm="#"+i+" J="+residual,wait=0);
	
	// Increment counter
	i++;
}