// p-Laplace problem, p>2, Rayleigh-Ritz
// Lecture notes, 10.3
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,w,phi;
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] nablaJR(real[int] & ucoeff)
// gradient of J^R applied to each phi function
{
  Vh uloc;
  uloc[] = ucoeff;
  varf F(unused,phi) = int2d(Th)((dx(uloc)^2+dy(uloc)^2)^((p-2)/2)*
                                 (dx(uloc)*dx(phi)+dy(uloc)*dy(phi))-f*phi)
    +on(1,unused=0); // it takes value 0 at points on the boundaries
  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)))-int2d(Th)(f*v)
  +on(1,u0=0);
// u0 is the solution of Poisson's problem
// simple solution, no preconditioner
u = u0;
NLCG(nablaJR,u[],nbiter=150,eps=-1e-8);
cout << "L2 error " << sqrt(int2d(Th)((uexact-u)^2)) << endl;
cout << "J error  " << abs(J(u[])-Juexact) << endl;
// we try with a preconditioner
// now we construct a bilinear simplification of HJR
// notice that u0 is Poisson's solution
varf HtildeJR(w,v) = int2d(Th)((dx(u0)^2+dy(u0)^2)^((p-2)/2)*
                               (dx(w)*dx(v)+dy(w)*dy(v)))
  +on(1,w=0);
// we compute the corresponding matrix and factorize once and for all
matrix H = HtildeJR(Vh,Vh,solver=Cholesky,factorize=1);
func real[int] applyprec(real[int] & vec)
{
  real[int] r = H^(-1)*vec;
  return r;
}
// one preconditioner
u = u0;
NLCG(nablaJR,precon=applyprec,u[],nbiter=150,eps=-1e-8);
cout << "L2 error " << sqrt(int2d(Th)((uexact-u)^2)) << endl;
cout << "J error  " << abs(J(u[])-Juexact) << endl;
// after a certain (small) number of NLCG iterations, it would be possible
// to compute a new preconditioner with its solution and then restart