//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mesh_adaptive/mesh_adaptive.edp
//
//  Discussion:
//
//    Demonstrate the adaptmesh facility.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    26 January 2015
//
//  Reference:
//
//    John Chrispell, Jason Howell,
//    Finite Element Approximation of Partial Differential Equations
//    Using FreeFem++, or, How I Learned to Stop Worrying and Love
//    Numerical Analysis.
//
//    Frederic Hecht,
//    New development in FreeFem++,
//    Journal of Numerical Mathematics,
//    Volume 20, Number 3-4, 2012, pages 251-265.
//
cout << "\n";
cout << "mesh_adaptive\n";
cout << "  FreeFem++ version\n";
cout << "  Repeatedly adapt a mesh according to some error indicator.\n";
//
//  Define Th, a mesh on the square [-1,+1]x[-1,+1].
//
mesh Th = square ( 5, 5, [2*x-1, 2*y-1] );
//
//  Define F, a mesh weight function.
//
func f = 10.0 * x^3 + y^3 + atan2 ( 0.0001, sin(5.0*y)-2.0*x );
//
//  Define Vh, a piecewise linear finite element space.
//
fespace Vh ( Th, P1 );
//
//  Define fh to be the interpolant in Vh of the function f.
//
Vh fh = f;
//
//  Make a variable to hold the plotfile names.
//
string plotfile;
//
//  Display the initial mesh.
//
int i = 0;
plotfile = "mesh_adaptive_" + i + ".ps";
plot ( Th, fh, wait = true, ps = plotfile );
//
//  Repeatedly adapt the mesh.
//
for ( int i = 1; i <= 3; i++ )
{
  Th = adaptmesh ( Th, fh );
//
//  Update fh to the new mesh.
//
  fh = f;
//
//  Plot this mesh.
//
  plotfile = "mesh_adaptive_" + i + ".ps";
  plot ( Th, fh, wait = true, ps = plotfile );
}
//
//  Terminate.
//
cout << "\n";
cout << "mesh_adaptive\n";
cout << "  Normal end of execution.\n";

exit ( 0 );