//  Discussion:
//
//    Solve a time-dependent boundary value problem on the unit square.
//
//    Use the midpoint method to handle the time derivative.
//
//    Our problem is posed on the unit square, with zero boundary conditions,
//    and having the exact solution
//
//      u = sin(pi x) sin(pi y) e^(-t)
//
//    The boundary value problem is posed as
//
//      du    d^2 u   d^2 u
//      -- -  ----- - ----- = (2 pi^2 - 1 ) sin(pi x) sin(pi y) e^(-t)
//      dt    dx^2    dy^2
//
//    with initial condition at t=0:
//
//      u(x,y,0) = sin(pi x) sin(pi y)
//
//    and boundary condition
//
//      u(x,y,t) = 0
//
//    A step of the midpoint method can be carried out by first taking
//    a backward Euler step of size dt/2, followed by a forward Euler
//    step of size dt/2.  
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 March 2021
//
//  Author:
//
//    John Burkardt
//
cout << endl;
cout << "midpoint:" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Solve a time-dependent boundary value problem in a square." << endl;
cout << "  Use the midpoint method to handle the time variation." << endl;
cout << "  Report L2 norm, L2 error norm, and L2 relative error norm." << endl;

real a;
real dt;
real h;
real hold;
int n;
int nlog2;
int nt;
real t;
real tm;
real tmin = 0.0;
real tmax = 2.0;
real told;

for ( nt = 10; nt <= 40; nt = nt * 2 )
{
  dt = ( tmax - tmin ) / nt;

  nlog2 = 6;
//
//  n: number of elements in each spatial direction.
//
  n = pow ( 2.0, nlog2 );
//
//  h: "width" of each element.
//
  h = 1.0 / n;
  cout << "  N =     " << n << endl;
  cout << "  H =     " << h << endl;
  cout << "  NT=     " << nt << endl;
  cout << "  DT =    " << dt << endl;
//
//  Th: the triangulation of the square.
//
  mesh Th = square ( n, n );
//
//  Define Vh, the finite element space.  "P2" = piecewise quadratic.
//
  fespace Vh ( Th, P2 );
//
//  uh: the finite element function that will approximate the solution.
//  uoldh: the solution at the previous time step.
//
  Vh uh;
  Vh umh;
  Vh uoldh;
  func uexact = sin ( pi * x ) * sin ( pi * y ) * exp ( - t );
//
//  vh: the test functions used to project the error.
//
  Vh vh;
//
//  f: the right hand side.
//
  func f = ( 2.0 * pi^2 - 1.0 ) * sin ( pi * x ) * sin ( pi * y ) * exp ( - t );
//
//  Define the problem:
//
//    We want to take a backward Euler step using a half timestep,
//    and then use a forward Euler step for the second half timestep.
//
//    umh will be the backward Euler solution estimate at the half time.
//
//    We have to be careful to evaluate the derivative at the half time "tm".
//
  problem halfbackwardeulerstep ( umh, vh ) =
     int2d ( Th ) 
     ( 
       umh * vh 
     )
     - int2d ( Th )
     (
       uoldh * vh
     )   
     + int2d ( Th )
     (
       + 0.5 * dt * dx ( umh ) * dx ( vh )
       + 0.5 * dt * dy ( umh ) * dy ( vh ) 
     )
    - int2d ( Th )
    (
      0.5 * dt * ( 2.0 * pi^2 - 1.0 ) * sin ( pi * x ) * sin ( pi * y ) * exp ( - tm ) * vh
    )
    + on ( 1, 2, 3, 4, umh = 0.0 );
//
//  Perform nt timesteps.
//
  cout << endl;
  cout << "  T  ||Exact||  ||Approx||   ||Error||  ||Error||/||Exact||" << endl;
  cout << endl;

  for ( int it = 0; it <= nt; it++ )
  {    
    if ( it == 0 )
    {
      t = tmin;
      uh = uexact;
    }
    else
    {
      told = t;
      t = ( ( nt - it ) * tmin + it * tmax ) / nt;
      tm = ( told + t ) / 2.0;
      uoldh = uh;
//
//  Take a backward Euler step with dt/2.
//
      halfbackwardeulerstep;
//
//  Follow this with a forward Euler step with dt/2.
//
      uh = 2.0 * umh - uoldh;
    }
//
//  Compute error.
//
    real uhnorm = sqrt 
    ( 
      int2d ( Th ) 
      ( 
        ( uh )^2 
      ) 
    );

    real uexactnorm = sqrt 
    ( 
      int2d ( Th ) 
      ( 
        ( uexact )^2 
      ) 
    );

    a = sqrt 
    ( 
      int2d ( Th ) 
      ( 
        ( uh - uexact )^2 
      ) 
    );

    cout << "  " << t 
         << "  " << uexactnorm 
         << "  " << uhnorm 
         << "  " << a 
         << "  " << a / uexactnorm << endl;
  }
//
//  Plot the computed solution at the final time.
//
  plot ( uh, fill = true, value = true, 
    ps = "midpoint_" + nt + ".ps", 
    cmm = "Midpoint solution, nt = " + nt );
}
//
//  Terminate.
//
cout << "\n";
cout << "midpoint:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );