//  Discussion:
//
//    -del p del(u) = 0  in [-1,+1]x[-1,+1].
//                u = g  on the boundary.
//
//    p = a1*I in Q1 and Q3
//      = a2*I in Q2 and Q4
//    f = 0
//    g = r^gamma * mu(theta)
//
//    mu(theta) = 
//    Q1: cos ( gamma * ( pi/2-sigma) ) * cos ( gamma * ( theta-pi/2+rho )    )
//    Q2: cos ( gamma *          rho  ) * cos ( gamma * ( theta-pi+sigma )    )
//    Q3: cos ( gamma *        sigma  ) * cos ( gamma * ( theta-pi-rho )      )
//    Q4: cos ( gamma * ( pi/2 - rho) ) * cos ( gamma * ( theta-3pi/2-sigma ) )
//
//    Parameter values:
//
//      a1 = 161.4476387975881
//      a2 = 1.0
//      gamma = 0.1
//      rho = pi / 4
//      sigma = -14.92256510455152
//
//    Note that the formulation in the Mitchell reference uses different
//    names for variables, a somewhat different formulation for mu(theta),
//    some typographical errors, and does not seem to produce the result
//    displayed in the figure.  The data and formulation from the Morin
//    reference seems more satisfactory.
//
//    I CANNOT get this problem to give me anything like the results 
//    displayed by Mitchell, nor can I explain the lack of symmetry.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_11/mitchell_11.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    09 February 2015
//
//  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.
//
//    Pedro Morin, Ricardo Nochetto, Kunibert Siebert,
//    Data oscillation and convergence of adaptive FEM,
//    SIAM Journal on Numerical Analysis,
//    Volume 38, Number 2, 2000, pages 466-488.
//
cout << "\n";
cout << "mitchell_11\n";
cout << "  FreeFem++ version\n";
cout << "  Solve the intersecting interface problem.\n";
//
//  Set parameters.
//
real a1 = 161.4476387975881;
real a2 = 1.0;
real gam = 0.1;
real rho = pi / 4.0;
real sigma = -14.92256510455152;
//
//  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 = 20;
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 exact solution function G = r^gam * mu(theta).
//
//func g {
func real g ( real xx, real yy )
{
  real r;
  real rg;
  real theta;
  real value;

  r = sqrt ( xx * xx + yy * yy );
  rg = pow ( r, gam );

  theta = atan2 ( yy, xx );

  if ( 0.0 <= xx && 0.0 <= yy )
  {
    value = rg * cos ( gam * ( 0.5 * pi - sigma ) ) 
               * cos ( gam * ( theta - 0.5 * pi + rho ) );
  }
  else if ( xx < 0.0 && 0.0 <= yy )
  {
    value = rg * cos ( gam * rho ) 
               * cos ( gam * ( theta - pi + sigma ) );
  }
  else if ( xx < 0.0 && yy < 0.0 )
  {
    value = rg * cos ( gam * sigma ) 
               * cos ( gam * ( theta - pi - rho ) );
  }
  else if ( 0.0 <= xx && yy < 0.0 )
  {
    value = rg * cos ( gam * ( 0.5 * pi - rho ) ) 
               * cos ( gam * ( theta - 1.5 * pi - sigma ) );
  }
  else
  {
    cout << "WHAT THE HELL? X = " << xx << "  Y = " << yy << "\n";
  }
  return value;
}
//
//  Define right hand side function F.
//
func f = 0.0;
//
//  Define coefficient function P.
//
func real p ( real xx, real yy )
{
  real value;
  if ( 0.0 <= xx * yy )
  {
    value = a1;
  }
  else
  {
    value = a2;
  }
  return value;
}
//
//  Solve the variational problem.
//
solve Laplace ( u, v ) 
  = int2d ( Th ) ( p(x,y) * ( dx(u)*dx(v) + dy(u)*dy(v) ) )
  - int2d ( Th ) ( f * v ) 
  + on ( 1, u = g(x,y) );
//
//  Plot the solution.
//
plot ( u, wait = true, fill = true, ps = "mitchell_11_u.ps" );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "mitchell_11_mesh.ps" );
//
//  Save the mesh file.
//
savemesh ( Th, "mitchell_11.msh" );
//
//  Terminate.
//
cout << "\n";
cout << "mitchell_11\n";
cout << "  Normal end of execution.\n";

exit ( 0 );