//  Discussion:
//
//    -uxx-uyy=f on the square [0,1]x[0,1];
//    u = g on boundary.
//
//    f = -gxx-gyy
//    g = 2^(4p) x^p (1-x)^p y^p (1-y)^p
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test01.edp
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    19 December 2014
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Frederic Hecht,
//    Freefem++,
//    Third Edition, version 3.22
//
//    William Mitchell,
//    A collection of 2D elliptic problems for testing adaptive
//    grid refinement algorithms,
//    Applied Mathematics and Computation,
//    Volume 220, 1 September 2013, pages 350-364.
//
border bottom ( t = 0.0, 1.0 )   { x = t;    y = 0.0; }
border right  ( t = 0.0, 1.0 )   { x = 1.0;  y = t; }
border top    ( t = 1.0, 0.0 )   { x = t;    y = +1.0; }
border left   ( t = 1.0, 0.0 )   { x = 0.0;  y = t; }
//
//  Define Th, the triangulation of the "left" side of the boundaries.
//
int n = 10;
mesh Th = buildmesh ( bottom ( n ) + right ( n ) + top ( n ) + left ( n ) );
//
//  Define Vh, the finite element space defined over Th, using P1 basis functions.
//
fespace Vh ( Th, P1 );
//
//  Define U, V, and F, piecewise continuous functions over Th.
//
Vh u;
Vh v;
//
//  Define function F.
//
int p = 10;
func f = - 2 ^ ( 4 * p ) * 
  (
    (       p * ( p - 1 ) * x ^ ( p - 2 ) * ( 1 - x ) ^ ( p     )
      - 2 * p *   p       * x ^ ( p - 1 ) * ( 1 - x ) ^ ( p - 1 )
      +     p * ( p - 1 ) * x ^ ( p     ) * ( 1 - x ) ^ ( p - 2 ) )
    *                       y ^ ( p )     * ( 1 - y ) ^ ( p )
    +
    (       p * ( p - 1 ) * y ^ ( p - 2 ) * ( 1 - y ) ^ ( p     )
      - 2 * p *   p       * y ^ ( p - 1 ) * ( 1 - y ) ^ ( p - 1 )
      +     p * ( p - 1 ) * y ^ ( p     ) * ( 1 - y ) ^ ( p - 2 ) )
    *                       x ^ ( p )     * ( 1 - x ) ^ ( p )
  );
//
//  Define the variational problem.
//  Solve the variational problem.
//
solve Laplace ( u, v ) 
  = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  - int2d ( Th ) ( f * v ) 
  + on ( bottom, u = 0.0 )
  + on ( right,  u = 0.0 )
  + on ( top,    u = 0.0 )
  + on ( left,   u = 0.0 );
//
//  Plot the solution.
//
plot ( u, wait = 1, fill = true, ps = "test01_u.eps" );
//
//  Plot the mesh.
//
plot ( Th, wait = 1, ps = "test01_mesh.eps" );
//
//  Save the mesh file.
//
savemesh ( Th, "test01.msh" );