//  Discussion:
//
//    -uxx-uyy = f on [-1,+1]x[-1,+1]
//    u = g on the boundary.
//
//    f = -gxx-gyy.
//    g = cos(pi*y/2)                       for x <= beta*(y+1)
//      = cos(pi*y/2)+(x-beta*(y+1))^alpha  for      beta*(y+1) < x
//
//    Suggested parameter values:
//    * alpha = 2.5, beta = 0.0
//    * alpha = 1.1, beta = 0.0
//    * alpha = 1.5, beta = 0.6
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_10/mitchell_10b.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    23 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_10b\n";
cout << "  FreeFem++ version\n";
cout << "  Solve the interior line singularity problem.\n";
//
//  Set parameters.
//
real alpha = 1.1;
real beta = 0.0;
//
//  Set the borders.
//
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 = 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 G.
//
func real g ( real alpha, real beta )
{
  real value;
  if ( x <= beta * ( y + 1.0 ) )
  {
    value = cos ( pi * y / 2.0 );
  }
  else
  {
    value = cos ( pi * y / 2.0 ) + pow ( x - beta * ( y + 1.0 ), alpha );
  }
  return value;
}
//
//  Define function F.
//
func real f ( real alpha, real beta )
{
  real value;
  if ( x <= beta * ( y + 1.0 ) )
  {
    value = 0.25 * pi * pi * cos ( pi * y / 2.0 );
  }
  else
  {
    value = 0.25 * pi * pi * cos ( pi * y / 2.0 ) 
      + alpha * ( alpha - 1 ) * ( 1.0 + beta * beta )
      * pow ( x - beta * ( y + 1.0 ), alpha - 2 );
  }
  return value;
}
//
//  Solve the variational problem.
//
solve Laplace ( u, v ) 
  = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  - int2d ( Th ) ( f ( alpha, beta ) * v ) 
  + on ( 1, u = g ( alpha, beta ) );
//
//  Plot the solution.
//
plot ( u, wait = true, fill = true, ps = "mitchell_10b_u.ps" );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "mitchell_10b_mesh.ps" );
//
//  Save the mesh file.
//
savemesh ( Th, "mitchell_10b.msh" );
//
//  Terminate.
//
cout << "\n";
cout << "mitchell_10b\n";
cout << "  Normal end of execution.\n";

exit ( 0 );