//  Discussion:
//
//    -eps * Laplace(u) + 2 du/dx + du/dy = f on [-1,+1]x[-1,+1]
//    u = g on the boundary.
//
//    f = -eps*(g'')+2 dg/dx+dg/dy
//    g = (1-exp(-(1-x)/eps)*(1-exp(-(1-y)/eps)*cos(pi(x+y))
//
//    Suggested parameter values:
//    * eps = 0.1
//    * eps = 0.001
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test06b.edp
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    18 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 = -1.0, +1.0 )   { x = t;    y = -1.0; label = 1; }
border right  ( t = -1.0, +1.0 )   { x = +1.0; y = t;    label = 1; }
border top    ( t = +1.0, -1.0 )   { x = t;    y = +1.0; label = 1; }
border left   ( t = +1.0, -1.0 )   { x = -1.0; y = t;    label = 1; }
//
//  Define Th, the triangulation of the "left" side of the boundaries.
//
int n = 40;
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 the problem parameter.
//  EPS = 0.001 is desired, but even eps = 0.01 fails to compute.
//
real eps = 0.01;
//
//  Define the right hand side function F.
//
func f = exp ( - 2.0 / eps ) / pow ( eps, 2 ) * 
  (
    - exp ( ( 1.0 - y ) / eps ) * 
    ( 
      -2.0 * exp ( 1.0 / eps ) * pow ( pi * eps, 2 )
     + exp ( x / eps ) * ( 1.0 + 2.0 * pow ( pi * eps, 2 ) ) 
    ) 
    * cos ( pi * ( x + y ) ) - 
    ( 
      3.0 * exp ( 2.0 / eps ) - exp ( ( 1.0 - x ) / eps )
      - exp ( ( 1.0 - y ) / eps ) - exp ( ( x + y ) / eps ) 
    )
    * eps * pi * sin ( pi * ( x + y ) ) 
  );
//
//  Define the Dirichlet boundary condition function G.
//
func g = ( 1.0 - exp ( - ( 1.0 - x ) / eps ) ) * 
         ( 1.0 - exp ( - ( 1.0 - y ) / eps ) ) * 
         cos ( pi * ( x + y ) );
//
//  Solve the variational problem.
//
solve Laplace ( u, v ) 
  = int2d ( Th ) ( eps*dx(u)*dx(v) + eps*dy(u)*dy(v) + 2*dx(u)*v + dy(u)*v )
  - int2d ( Th ) ( f * v ) 
  + on ( 1, u = g );
//
//  Plot the solution.
//
plot ( u, wait = 1, fill = true, ps = "test06b_u.eps" );
//
//  Plot the mesh.
//
plot ( Th, wait = 1, ps = "test06b_mesh.eps" );
//
//  Save the mesh file.
//
savemesh ( Th, "test06b.msh" );