//  Discussion:
//
//    This example solves the Poisson equation repeatedly, using each
//    solution to adapt the mesh before the next solution is computed.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/poisson_adaptive/poisson_adaptive.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 June 2015
//
//  Reference:
//
//    Frederic Hecht,
//    Freefem++,
//    Third Edition, version 3.22
//
cout << "\n";
cout << "poisson_adaptive\n";
cout << "  FreeFem++ version\n";
cout << "  Solve the Poisson equation in the L-shaped region\n";
cout << "  using adaptive refinement to reduce the error.\n";
//
//  Define the border.
//
border aaa (t=0,1)   {x=t;y=0;};
border bbb (t=0,0.5) {x=1;y=t;};
border ccc (t=0,0.5) {x=1-t;y=0.5;};
border ddd (t=0.5,1) {x=0.5;y=t;};
border eee (t=0.5,1) {x=1-t;y=1;};
border fff (t=0,1)   {x=0;y=1-t;};
//
//  Define Th, the triangulation of the domain to the left of the border.
//
mesh Th = buildmesh ( aaa(6) + bbb(4) + ccc(4) + ddd(4) + eee(4) + fff(6) );
//
//  Plot the initial mesh.
//
plot ( Th, wait = true, ps = "poisson_adaptive_mesh_0.ps" );
//
//  Define Vh, the finite element space defined over Th, using P1 basis functions.
//
fespace Vh ( Th, P1 );
//
//  Define U and V as elements of the space VH.
//  Moreover, U is the zero function.
//
Vh u = 0.0;
Vh v;
//
//  Define F, the PDE right hand side function, and
//  G, used to define the Dirichlet boundary condition.
//
func f = 1.0;
func g = 0.0;
//
//  Define the weak form of the problem.
//
problem Problem1 ( u, v, solver = CG, eps = -1.0e-6 ) 
  = int2d (Th) ( dx(u)*dx(v) + dy(u)*dy(v) ) 
  + int2d (Th) ( v * f ) 
  + on ( aaa, bbb, ccc, ddd, eee, fff, u = g );
//
//  Carry out an iteration in which the mesh is adaptively refined.
//
real error = 0.1;
real coef = 0.1^(1.0/5.0);

for ( int it = 0; it < 10; it++ )
{   
  Problem1;
  plot ( u, wait = true, ps = "poisson_adaptive_uh_"+it+".ps" );
  Th = adaptmesh ( Th, u, inquire = 1, err = error );
  plot ( Th, wait = true, ps = "poisson_adaptive_mesh_"+(it+1)+".ps" );
  error = error * coef;
}
//
//  Terminate.
//
cout << "\n";
cout << "poisson_adaptive:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );