//  Discussion:
//
//    Solve the steady Poisson equation in the unit square.
//
//      - uxx - uyy = f(x,y)
//
//    with boundary condition
//
//      u = 0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    12 February 2021
//
//  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 << "poisson_2d():\n";
cout << "  FreeFem++ version.\n";
cout << "  Solve the Poisson equation in a square region.\n";
//
//  n: number of nodes in each spatial direction.
//
int nx = 11;
int ny = 11;
//
//  h: spacing between nodes.
//
real hx = 1.0 / nx;
real hy = 1.0 / ny;
//
//  Th: the triangulation of the unit square.
//
mesh Th = square ( nx, ny );
//
//  Plot the mesh.
//
plot ( Th, 
  ps = "poisson_2d_mesh.ps", 
  cmm = "poisson_2d: the mesh" );
//
//  Define Vh, the finite element space.
//  Here, "P1" requests piecewise linear functions.
//
fespace Vh ( Th, P1 );
//
//  uh: the trial function that approximates the solution.
//
Vh uh;
//
//  vh: the test functions that multiply the equation.
//
Vh vh;
//
//  f: the right hand side function f(x,y).
//
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));
//
//  g: the boundary condition function g(x,y).
//
func g = 0.0;
//
//  define "poisson()" to be the problem to be solved.
//
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 for uh.
//
poisson;
//
//  Plot the computed solution.
//
plot ( uh, fill = true, value = true, 
  ps = "poisson_2d_fem.ps", 
  cmm = "poisson_2d: FEM solution" );
//
//  For this problem, we know the exact solution.
//  u: the exact solution, projected into finite element space.
//
Vh u = sin ( 5 * pi * x * ( 1 - x ) ) * sin ( 4 * pi * y * ( 1 - y ) );
//
//  Plot the exact solution.
//
plot ( u, fill = true, value = true,  
  ps = "poisson_2d_exact.ps", 
  cmm = "poisson_2d: exact" );
//
//  Compute the L2 error.
//
real l2error = sqrt ( int2d ( Th ) ( ( u - uh )^2 ) );
cout << endl;
cout << "  Hx = " << hx << endl;
cout << "  Hy = " << hy << endl;
cout << "  L2 error = " << l2error << endl;
//
//  Terminate.
//
cout << "\n";
cout << "poisson_2d():\n";
cout << "  Normal end of execution.\n";

exit ( 0 );