//  Discussion:
//
//    -d/dx p dudx - d/dy q dudy = f on [0,8.4]x[0,24]
//    Mixed boundary conditions
//
//    f = given
//    g unknown
//
//    P, Q and F are piecewise constant functions:
//
//    Domain  P     Q       F
//    ------ ----  ----    ----
//    1      25.0  25.0     0.0
//    2       7.0   0.8     1.0
//    3       5.0   0.0001  1.0
//    4       0.2   0.2     0.0
//    5       0.05  0.05    0.0
//
//   Y8  +------------------+
//       | K = 1            |
//   Y7  +---------------+  |
//       | K = 5         |  |
//   Y6  +--------+---+  |  |
//       | K = 2  |   |  |  |
//   Y5  +--------+ K |  |  |
//       | K = 3  | = |  |  |
//   Y4  +--------+ 4 |  |  |
//       | K = 2  |   |  |  |
//   Y3  +--------+   |  |  |
//       | K = 5  |   |  |  |
//   Y2  +--------+---+--+  |
//       |                  |
//   Y1  +------------------+
//      X1       X2  X3 X4 X5
//
//    The rectangular subdomains are characterized by 5 X values:
//      0.0, 6.1, 6.5, 8.0, 8.4
//    and 8 Y values:
//      0.0, 0.8, 1.6, 3.6, 18.8, 21.2, 23.2, 24.0
//
//    The boundary conditions are of the form:
//      p dudx dyds - q dudy dxds + c * u = g
//
//    using the following values:
//
//      side     C    gN    Boundary condition
//      ------- ---   ---   ------------------
//      bottom  3.0   1.0:  - q dudy + 3 u = 1
//      right   2.0   2.0:    p dudx + 2 u = 2
//      top     1.0   3.0:    q dudy +   u = 3
//      left    0.0   0.0:  - p dudx       = 0
//
//    The mesh generated by BAMG will assign the following boundary labels:
//
//      side     label
//      -------  -----
//      bottom       1
//      right        2
//      top          3
//      left         4
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_05/mitchell_05.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 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.
//
cout << "\n";
cout << "mitchell_05\n";
cout << "  FreeFem++ version\n";
cout << "  Mitchel test problem 5, the battery.\n";
//
//  Set parameters.
//
real x1 = 0.0;
real x2 = 6.1;
real x3 = 6.5;
real x4 = 8.0;
real x5 = 8.4;
real y1 = 0.0;
real y2 = 0.8;
real y3 = 1.6;
real y4 = 3.6;
real y5 = 18.8;
real y6 = 21.2;
real y7 = 23.2;
real y8 = 24.0;
//
//  Read the mesh.
//
mesh Th = readmesh ( "battery.msh" );
//
//  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 function F.
//
func real f ( real xx, real yy )
{
  real value;

  value = 0.0;
//
//  Region 1:
//
  if ( x4 <= xx || y7 <= yy || yy <= y2 )
  {
    value = 0.0;
  }
//
//  Region 2:
//
  else if ( y3 <= yy && yy <= y4 && xx <= x2 )
  {
    value = 1.0;
  }
  else if ( y5 <= yy && yy <= y6 && xx <= x2 )
  {
    value = 1.0;
  }
//
//  Region 3;
//
  else if ( y4 <= yy && yy <= y5 && xx <= x2 )
  {
    value = 1.0;
  }
//
//  Region 4:
//
  else if ( y2 <= yy && yy <= y6 && x2 <= xx && xx <= x3 )
  {
    value = 0.0;
  }
//
//  Region 5:
//
  else if ( y6 <= yy && yy <= y7 && xx <= x4 )
  {
    value = 0.0;
  }
  else if ( y2 <= yy && yy <= y7 && x3 <= xx && xx <= x4 )
  {
    value = 0.0;
  }
  else
  {
    cout << "F: Stray values  X = " << x << "  Y = " << y << "\n";
  }
  return value;
}
//
//  Define the boundary condition functions G1, G2, G3, G4:
//
func g1 = 1.0;
func g2 = 2.0;
func g3 = 3.0;
func g4 = 0.0;
//
//  Define the coefficient function P.
//
func real p ( real xx, real yy )
{
  real value;

  value = 0.0;
//
//  Region 1:
//
  if ( x4 <= xx || y7 <= yy || yy <= y2 )
  {
    value = 25.0;
  }
//
//  Region 2:
//
  else if ( y3 <= yy && yy <= y4 && xx <= x2 )
  {
    value = 7.0;
  }
  else if ( y5 <= yy && yy <= y6 && xx <= x2 )
  {
    value = 7.0;
  }
//
//  Region 3;
//
  else if ( y4 <= yy && yy <= y5 && xx <= x2 )
  {
    value = 5.0;
  }
//
//  Region 4:
//
  else if ( y2 <= yy && yy <= y6 && x2 <= xx && xx <= x3 )
  {
    value = 0.2;
  }
//
//  Region 5:
//
  else if ( y6 <= yy && yy <= y7 && xx <= x4 )
  {
    value = 0.05;
  }
  else if ( y2 <= yy && yy <= y7 && x3 <= xx && xx <= x4 )
  {
    value = 0.05;
  }
  else if ( y2 <= yy && yy <= y3 && x1 <= xx && xx <= x2 )
  {
    value = 0.05;
  }
  else
  {
    cout << "P: Stray values:  X = " << x << "  Y = " << y << "\n";
  }
  return value;
}
//
//  Define the coefficient function Q.
//
func real q ( real xx, real yy )
{
  real value;

  value = 0.0;
//
//  Region 1:
//
  if ( x4 <= xx || y7 <= yy || yy <= y2 )
  {
    value = 25.0;
  }
//
//  Region 2:
//
  else if ( y3 <= yy && yy <= y4 && xx <= x2 )
  {
    value = 0.8;
  }
  else if ( y5 <= yy && yy <= y6 && xx <= x2 )
  {
    value = 0.8;
  }
//
//  Region 3;
//
  else if ( y4 <= yy && yy <= y5 && xx <= x2 )
  {
    value = 0.0001;
  }
//
//  Region 4:
//
  else if ( y2 <= yy && yy <= y6 && x2 <= xx && xx <= x3 )
  {
    value = 0.2;
  }
//
//  Region 5:
//  Two regions, one of which is described by two rectangles.
//
  else if ( y6 <= yy && yy <= y7 && xx <= x4 )
  {
    value = 0.05;
  }
  else if ( y2 <= yy && yy <= y7 && x3 <= xx && xx <= x4 )
  {
    value = 0.05;
  }
  else if ( y2 <= yy && yy <= y3 && x1 <= xx && xx <= x2 )
  {
    value = 0.05;
  }
  else
  {
    cout << "Q: Stray values:  X = " << x << "  Y = " << y << "\n";
  }
  return value;
}
//
//  Solve the variational problem with Robin boundary conditions.
//
solve Laplace ( u, v ) 
  = int2d ( Th ) ( p(x,y) * dx(u)*dx(v) + q(x,y) * dy(u)*dy(v) )
  - int2d ( Th ) ( f(x,y) * v ) 
  - int1d ( Th, 1 ) (   3.0 * u * v )
  - int1d ( Th, 1 ) ( - 1.0     * v )
  - int1d ( Th, 2 ) (   2.0 * u * v )
  - int1d ( Th, 2 ) ( - 2.0     * v )
  - int1d ( Th, 3 ) (         u * v )
  - int1d ( Th, 3 ) ( - 3.0     * v );
//
//  Plot the solution.
//
plot ( u, wait = true, fill = true, ps = "mitchell_05_u.eps" );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "mitchell_05_mesh.eps" );
//
//  Save the mesh file.
//
savemesh ( Th, "mitchell_05.msh" );
//
//  Terminate.
//
cout << "\n";
cout << "mitchell_05:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );