// p-Laplace problem, p>2, Rayleigh-Ritz
int nint = 100;
real p = 4.0; // p > 2.0 otherwise NLCG input "too" exact. Do not forget .0
real r = 1.0;
border Gamma (t=0,2*pi) { x = r * cos(t); y = r * sin(t); label = 1; }
mesh Th = buildmesh(Gamma(nint));
fespace Vh(Th,P1);
Vh u0,u,v;
func f = 1.0;
func uexact = (p-1)/p*(0.5)^(1/(p-1))*(1-(x^2+y^2)^(p/(2*p-2)));
real Juexact = 2*pi*(0.5^(p/(p-1.0))*(p-1)/(3*p-2.0)/p-
                      (p-1)/p*0.5^(1.0/(p-1))*(0.5-(p-1)/(3*p-2)));
func real J(real[int] & ucoeff)
{
  Vh uloc;
  uloc[] = ucoeff;
  return int2d(Th)((dx(uloc)^2+dy(uloc)^2)^(p/2)/p-f*uloc);
}
func real[int] deltaJR(real[int] & ucoeff)
// gradient of J^R applied to each phi function
{
  Vh uloc;
  uloc[] = ucoeff;
  varf F(unused,v) = int2d(Th)((dx(uloc)^2+dy(uloc)^2)^((p-2)/2)*
                               (dx(uloc)*dx(v)+dy(uloc)*dy(v))-f*v)
    +on(1,unused=0); // it takes value 0 at points on the boundaries
  /*
  In NLCG, x_m = x_{m-1} + \alpha_{m-1} d_{m-1}
  d_{m-1} is a composition of nablaJR directions
  Therefore, if x_0 has the right boundary conditions, we should not touch them
  */
  real[int] r = F(0,Vh);
  return r;
}
// initial solution given by Poisson (p=2)
solve Poisson(u0,v) = int2d(Th)((dx(u0)*dx(v)+dy(u0)*dy(v))) // bilinear form
  -int2d(Th)(f*v) // linear functional
  +on(1,u0=0); // boundary conditions
/*
The default solver is sparsesolver (multi-frontal LU)
It is possible to specify solver=CG, and possibly eps
*/
// u0 is the solution of Poisson's problem
// simple solution, no preconditioner
u = u0;
NLCG(deltaJR,u[],nbiter=200,eps=-1e-4); // a simple minimizer
cout << "L2 error " << sqrt(int2d(Th)((uexact-u)^2)) << endl;
cout << "J error  " << abs(J(u[])-Juexact) << endl;
plot(u);