// p-Laplace problem, p>2, Galerkin approach
int nint = 100;
real p = 4.0; // p>2 otherwise NLCG input "too" exact
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 u,v,phi,delta;
real tol = 1e-8;
func f = 1.0; 
func uexact = (p-1.0)/p*(0.5)^(1/(p-1.0))*(1-(x^2+y^2)^(p/(2*p-2.0)));
real Ju = 2*pi*(0.5^(p/(p-1.0))*(p-1.0)/(3*p-2.0)/p-
		(p-1.0)/p*0.5^(1.0/(p-1.0))*(0.5-(p-1.0)/(3*p-2.0)));
func real J(real[int] & u)
{
  Vh uloc;
  uloc[] = u;
  return int2d(Th)((dx(uloc)^2+dy(uloc)^2)^(p/2)/p-f*uloc);
}
// initial solution given by Poisson (p=2)
problem Poisson(u,v) = int2d(Th)((dx(u)*dx(v)+dy(u)*dy(v)))
  -int2d(Th)(f*v)
  +on(1,u=0);
Poisson; // now u is the solution of Poisson
// newton
problem Newtonit(delta,phi) = int2d(Th)((p-2)*(dx(u)^2+dy(u)^2)^((p-4.0)/2)*
				      (dx(u)*dx(delta)+dy(u)*dy(delta))*
				      (dx(u)*dx(phi)+dy(u)*dy(phi)))
  +int2d(Th)((dx(u)^2+dy(u)^2)^((p-2.0)/2)*
  	     (dx(delta)*dx(phi)+dy(delta)*dy(phi)))
  +int2d(Th)((dx(u)^2+dy(u)^2)^((p-2.0)/2)*(dx(u)*dx(phi)+dy(u)*dy(phi)))-int2d(Th)(f*phi)
  +on(1,delta=0);
Newtonit;
cout << "delta " << delta[].linfty << endl;
while (delta[].linfty > tol)
  {
    u = u+delta;
    Newtonit;
    cout << "delta " << delta[].linfty << endl;
  }
cout << "L2 error " << sqrt(int2d(Th)((uexact-u)^2)) << endl;
cout << "J error  " << abs(J(u[])-Ju) << endl;