//  Discussion:
//
//    Solve a time-dependent boundary value problem on the unit square.
//
//    Use the backward Euler 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
//
//    The backward Euler method allows us to take a time step dt,
//    using the solution uold at the previous step, replacing the
//    time derivative by (u-uold)/dt.  As a finite element expression,
//    the problem is to solve the following integral equation for u:
//
//      Integral ( u - uold ) v + dt ( ux vx + uy vy - f(x,y,t,u) v ) = 0
//
//     where v is any finite element space basis vector.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 March 2021
//
//  Author:
//
//    John Burkardt
//
cout << endl;
cout << "backward_euler:" << endl;
cout << "  FreeFem++ version." << endl;
cout << "  Solve a time-dependent boundary value problem in a square." << endl;
cout << "  Use the backward Euler 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 tmin = 0.0;
real tmax = 2.0;

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 quadratics.
//
  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 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:
//
//  u - uold 
//  ---------  -  uxx - uyy = f(t,u)
//     dt
//
//    ( uh - uoldh ) vh + dt * uhx vhx + dt * uhy vhy - dt f vh = 0
//
  problem backwardeulerstep ( uh, vh ) =
     int2d ( Th ) 
     ( 
       uh * vh 
     )
     - int2d ( Th )
     (
       uoldh * vh
     )   
     + int2d ( Th )
     (
       + dt * dx ( uh ) * dx ( vh )
       + dt * dy ( uh ) * dy ( vh ) 
     )
    - int2d ( Th )
    (
      dt * f * vh
    )
    + on ( 1, 2, 3, 4, uh = 0.0 );
//
//  Perform nt timesteps.
//
  cout << endl;
  cout << "  T  ||Exact||  ||Approx||   ||Error||  ||Error||/||Exact||" << endl;
  cout << endl;

  for ( int it = 0; it <= nt; it++ )
  {
    t = ( ( nt - it ) * tmin + it * tmax ) / nt;

    if ( it == 0 )
    {
      uh = uexact;
    }
    else
    {
      uoldh = uh;
      backwardeulerstep;
    }
//
//  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 = "backward_euler_" + nt + ".ps", 
    cmm = "Backward Euler solution, nt = " + nt );
}
//
//  Terminate.
//
cout << "\n";
cout << "backward_euler:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );