//  Discussion:
//
//    - uxx - uyy = f on the unit square.
//    u = g on the boundary.
//
//    f = - gxx - gyy
//    g = atan ( alpha * (sqrt((x-xc)^2+(y-yc)^2)-r0) )
//
//    Suggested parameter values:
//    * alpha =   20, xc = -0.05, yc = -0.05, r0 = 0.7
//    * alpha = 1000, xc = -0.05, yc = -0.05, r0 = 0.7
//    * alpha = 1000, xc =  1.5,  yc =  0.25, r0 = 0.92
//    * alpha =   50, xc =  0.5,  yc =  0.5,  r0 = 0.25
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_09/mitchell_09b.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    22 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.
//
cout << "\n";
cout << "mitchell_09b\n";
cout << "  FreeFem++ version\n";
cout << "  The wave front problem.\n";

border bottom ( t = 0.0, 1.0 )   { x = t;   y = 0.0;  label = 1; }
border right  ( t = 0.0, 1.0 )   { x = 1.0; y = t;    label = 1; }
border top    ( t = 1.0, 0.0 )   { x = t;   y = 1.0;  label = 1; }
border left   ( t = 1.0, 0.0 )   { x = 0.0; y = t;    label = 1; }
//
//  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 problem parameters.
//
real alpha = 1000.0;
real xc = - 0.05;
real yc = - 0.05;
real r0 = 0.7;
//
//  Define the right hand side function F.
//
func f = 
  ( 
    - alpha + pow ( alpha, 3 ) 
    * ( - r0 * r0 + pow ( x - xc, 2 ) + pow ( y - yc, 2  ) )
  )
  /
  ( 
    sqrt ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) 
    * pow ( 
      -1.0 + alpha * alpha * 
      ( 
        - r0 * r0 - pow ( x - xc, 2 ) - pow ( y - yc, 2 )
        + sqrt ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) 
      )
    , 2 ) 
  );
//
//  Define the boundary function G.
//
func g = atan ( alpha * ( sqrt ( pow ( x - xc, 2 ) + pow ( y - yc, 2 ) ) - r0 ) );
//
//  Solve the variational problem.
//
solve Laplace ( u, v ) 
  = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  - int2d ( Th ) ( f * v ) 
  + on ( 1, u = g );
//
//  Plot the solution.
//
plot ( u, wait = true, fill = true, ps = "mitchell_09b_u.ps" );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "mitchell_09b_mesh.ps" );
//
//  Save the mesh file.
//
savemesh ( Th, "mitchell_09b.msh" );
//
//  Terminate.
//
cout << "\n";
cout << "mitchell_09b\n";
cout << "  Normal end of execution.\n";

exit ( 0 );