//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/hecht_freefem++/adapt.edp
//
//  Discussion:
//
//    The domain is the L shaped region.
//
//    PDE inside L:
//
//      - uxx - uyy = f(x,y) = x - y
//
//    BC on the boundary of Omega:
//
//      du/dn = 0
//
//    Because this is a full Neumann problem, it is singular.
//    To avoid singularity, we want to add the condition:
//
//      Integral u = 0
//
//    This equation is incorporated into the variational problem
//    with a tiny multiplier.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    21 January 2015
//
//  Reference:
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//

//
//  Although we go to some trouble to specify that all sides of the
//  region will have a label of "1", we never actually do anything with
//  this label.
//
int[int] lab = [1,1,1,1];
//
//  Create a mesh for a square, using 6 nodes on each side.
//
mesh Th = square ( 6, 6, label = lab );
//
//  Restrict the mesh to the L-shaped subregion of the square.
//  The newly created borders also get the label "1".
//
Th = trunc ( Th, x < 0.5 | y < 0.5, label = 1 );
//
//  Our finite element space will be P1 (piecewise linear Lagrange functions)
//
fespace Vh ( Th, P1 );
//
//  Let the symbols u and v stand for sample functions in this space.
//
Vh u, v;
//
//  Define the problem, and call it "Poisson".
//    Integral eps * u*v + ux*vx + uy*vy = Integral (x-y)*v
//
//  Notice there are no boundary conditions specified.
//
//  The additional eps*u*v term removes the singularity.
//
problem Poisson ( u, v, solver = CG, eps = 1.0e-6 ) =
    int2d ( Th, qforder=2 ) ( 1.0e-10*u*v + dx(u)*dx(v) + dy(u)*dy(v) ) 
  + int2d ( Th, qforder=2 ) ( (x-y)*v ); 
//
//  TOL is a tolerance for the adaptive mesh.
//  
real tol = 0.01;
//
//  Repeatedly solve the problem.
//
for ( int i = 0; i < 4; i++ )
{
  Poisson;
  plot ( u, wait = 1 );
  Th = adaptmesh ( Th, u, err = tol );
  plot ( Th, wait = 1 );
//
//  The next statement transfers U from the old to the new mesh.
//
  u = u;
//
//  Solve the problem again, with a reduced error tolerance.
//
  tol = tol / 2.0;
};