//  Location:
//
//    https://people.sc.fsu.edu/~jburkardt/freefem_src/error_norms/error_norms.edp
//
//  Discussion:
//
//    Error calculation for Poisson equation in a square.
//
//      - uxx - uyy = f
//
//    Assuming we know the exact solution, compute the L2 and H1 errors.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    25 May 2020
//
//  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 << "error_norms:\n";
cout << "  FreeFem++ version\n";
cout << "  Compute L2 and H1 error norms of an approximate solution.\n";
//
//  n = 8 corresponds to 9 nodes in X and Y directions, and h = 1/8.
//
int n = 8;
//
//  Define Th, the triangulation of the square.
//
mesh Th = square ( n, n );
//
//  Display the mesh.
//
plot ( Th, ps = "mesh.ps" );
//
//  Define Vh, the finite element space, using piecewise linears.
//
fespace Vh ( Th, P1 );
//
//  Define uh, the trial function.
//
Vh uh;
//
//  Define vh, test function.
//
Vh vh;
//
//  Define F(X,Y), the right hand side function.
//
func f = 
  +      sin(5*pi*x*(1-x)) *      sin(4*pi*y*(1-y)) * pow(5*pi*(1-x)-5*pi*x,2)
  + 10 * cos(5*pi*x*(1-x)) * pi * sin(4*pi*y*(1-y))
  +      sin(5*pi*x*(1-x)) *      sin(4*pi*y*(1-y)) * pow(4*pi*(1-y)-4*pi*y,2) 
  + 8  * sin(5*pi*x*(1-x)) * pi * cos(4*pi*y*(1-y));
//
//  Define G(X,Y), for the boundary condition.
//
func g = 0.0;
//
//  Define the problem.
//
problem poisson ( uh, vh ) 
  = int2d ( Th ) 
  ( 
    dx(uh) * dx(vh) 
  + dy(uh) * dy(vh)
  )
  - int2d ( Th )
  (
    f * vh 
  )
  + on ( 1, 2, 3, 4, uh = g );
//
//  Solve the problem.
//
poisson;
//
//  Plot the approximate solution.
//
plot ( uh, fill = true, value = true, ps = "approximate.ps" );
//
//  Define the true solution.
//
func u =                      sin(5*pi*x*(1-x))                    *sin(4*pi*y*(1-y));
func ux = (5*pi*(1-x)-5*pi*x)*cos(5*pi*x*(1-x))                    *sin(4*pi*y*(1-y));
func uy =                     sin(5*pi*x*(1-x))*(4*pi*(1-y)-4*pi*y)*cos(4*pi*y*(1-y));
//
//  Project true solution.
//
Vh up = sin ( 5 * pi * x * ( 1 - x ) ) * sin ( 4 * pi * y * ( 1 - y ) );
//
//  Plot the projected exact solution.
//
plot ( up, fill = true, value = true, ps = "exact.ps" );
//
//  Compute and print the L2 error.
//
real ul2 = sqrt 
( 
  int2d ( Th ) 
  ( 
    ( u - uh )^2 
  ) 
);
cout << "L2 error = " << ul2 << endl;
//
//  Compute and print the H1 seminorm error.
//
real uh1 = sqrt 
( 
  int2d ( Th )
  (
    ( ux - dx(uh) )^2 + ( uy - dy(uh) )^2
  )
);
cout << "H1 seminorm error = " << uh1 << endl;
//
//  Terminate.
//
cout << "\n";
cout << "error_norms:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );