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

// Task 1: Run the program
// Task 2: Understand how the derivative dJ is obtained (recall the 
//   theoretical exercise) and how Newton's method is implemented
// Task 3: Change the boundary to only Dirichlet conditions
// Task 4: What is the step-size in this program? 
//   Can we change the step-size and if yes to what values?
//   To what would other step-sizes correspond? (Hint: globalization)
// Task 5: Run gradient method and Newton's method and compare
//   the number of iterations 


// Number of boundary edges 
int N=50;	

// Definition of volume and surface traction
func f=1;
func g=-1;	

// power-law index (p > 2)
int p=4;

// Tolerance for the Newton solver
real lowertolerance = 1.e-10;	

macro vh() P1 //

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);
// Similar to the gradient method:
// du: Newton update
// phi: test function
// u: solution u^k
Vh du, phi, u;		

// Set initial Newton Newton guess u^0 to 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) //


// Implement the Newton problem
// As in chapter 3, page 37 described we now need 
// also the second derivative of the functional J.
// Specifically we need to determine the update du.
// This is done by solving a linear equation system of 
// the form J''(u)(phi)du = J'(u) (i.e., something like Ax = b).
// Such a problem is similar to solving a `standard' FE problem
// as e.g., Poisson's problem
problem Jlin(du,phi)=int2d(Th)((1+sqrt(grad(u)'*grad(u))^(p-2))*(grad(du)'*grad(phi))
		    +(p-2)*sqrt(grad(u)'*grad(u))^(p-4)*(grad(u)'*grad(du))*(grad(u)'*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
// Initialize the norm
real funcnorm = 1.0e+10;
int i=0;

// Final loop
while(funcnorm > lowertolerance)
{
	// Solve the (linear) problem, which yields the Newton update du
	Jlin;
	// Add update to previous solution u
	u=u+du;

	// Compute the norm of the solution
	funcnorm=sqrt(int2d(Th)(du*du));

	// Evaluate the value of the functional w.r.t. to the current solution
	real evalJ = J(u);
	cout<<"#"<<i<<" J ="<<evalJ << " " <<funcnorm<<endl;
	plot(u,cmm="#"+i+" J="+evalJ+" normw = "+funcnorm,wait=0);
	i++;
}