//  Discussion:
//
//    Solve a simple boundary value problem on the unit square.
//    Repeat the process after halving h.
//    Report the errors and error convergence rates.
//
//  Location:
//
//    https://people.sc.fsu.edu/~jburkardt/freefem_src/convergence/convergence.edp  
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    11 March 2021
//
//  Author:
//
//    John Burkardt
//
cout << endl;
cout << "convergence:" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Solve a boundary value problem in a square." << endl;
cout << "  Repeatedly refine the mesh, and report the L2 error norm." << endl;
//
//  Declare some convergence variables here.
//
real a;
real aold;
real arate;
real b;
real bold;
real brate;
real c;
real cold;
real crate;
real h;
real hold;
//
//  Loop for h = 1/8, 1/16, ...
//
for ( int log2 = 3; log2 <= 7; log2++ )
{
//
//  n: number of elements in each spatial direction.
//
  int n = pow ( 2.0, log2 );
//
//  h: characteristic length of each element.
//
  if ( 3 < log2 )
  {
    hold = h;
  }
  h = 1.0 / n;
//
//  Th: the triangulation of the region.
//
  mesh Th = square ( n, n );
//
//  Plot the mesh the first time.
//
  if ( log2 == 3 )
  {
    plot ( Th, wait = false, 
      cmm = "Initial mesh of the region",
      ps = "convergence_mesh.ps" );
  }
//
//  Define Vh, the finite element space.  "P1" = piecewise linears.
//
  fespace Vh ( Th, P1 );
//
//  ch: the finite element function that will approximate the solution.
//
  Vh ch;
//
//  vh: the test functions used to project the error.
//
  Vh vh;
//
//  ce: the exact solution, projected into the finite element space.
//
  Vh ce = ( x - x^2 ) * sin ( pi * y );
//
//  f: the right hand side.
//
  func f = 
      ( 1.0 - 2.0 * x )        * sin ( pi * y )
    + ( x - x^2 )       * pi   * cos ( pi * y )
    + 2.0                      * sin ( pi * y )
    + ( x - x^2 )       * pi^2 * sin ( pi * y );
//
//  Define the problem.
//
  problem bvp ( ch, vh ) 
    = int2d ( Th )
    (
      ( dx(ch) + dy(ch) ) * vh
      + dx(ch) * dx(vh) + dy(ch) * dy(vh) 
    )
    + int2d ( Th )
    (
      - f * vh 
    )
    + on ( 1, 2, 3, 4, ch = ce );
//
//  Solve the problem.
//
  bvp;
//
//  Compute error and convergence rates.
//
  if ( 3 < log2 )
  {
    aold = a;
    bold = b;
    cold = c;
  }

  a = sqrt 
  ( 
    int2d ( Th ) 
    ( 
      ( ch - ce )^2
    ) 
  );

  b = sqrt
  ( 
    int2d ( Th ) ( ( dx(ch) - dx(ce) )^2 ) 
  + int2d ( Th ) ( ( dy(ch) - dy(ce) )^2 )
  );

  c = sqrt ( a^2 + b^2 );

  if ( 3 < log2 )
  {
    arate = ( log ( aold ) - log ( a ) ) / ( log ( hold ) - log ( h ) );
    brate = ( log ( bold ) - log ( b ) ) / ( log ( hold ) - log ( h ) );
    crate = ( log ( cold ) - log ( c ) ) / ( log ( hold ) - log ( h ) );
  }

  cout << endl;
  cout << "  H =     " << h << endl;
  if ( 3 == log2 )
  {
    cout << "  L2 =    " << a << endl;
    cout << "  H0 =    " << b << endl;
    cout << "  H1 =    " << c << endl;
  }
  else
  {
    cout << "  L2 =    " << a << "  L2rate = " << arate << endl;
    cout << "  H0 =    " << b << "  H0rate = " << brate << endl;
    cout << "  H1 =    " << c << "  H1rate = " << crate << endl;
  }
}
//
//  Terminate.
//
cout << "\n";
cout << "convergence:\n";
cout << "  Normal end of execution.\n";