//  Discussion:
//
//    -uxx - uyy = f on the upside down L shaped region
//    u = g on the boundary.
//
//    f = - gxx - gyy
//    g = awful formula
//
//    Suggested parameter values:
//      omega = 3 * pi / 2
//      xw =  0.0
//      yw = -0.75
//      r0 =  0.75
//      alphaw = 200.0
//      xp = sqrt(5)/4
//      yp = -0.25
//      alphap = 1000
//      epsilon = 0.01
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/examples/mitchell_freefem++/test12.edp
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    26 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.
//

//
//  Set parameters.
//

//
//  Set the borders of the L-shaped region.
//
border e1 ( t = -1.0,  0.0 ) { x = t;    y = -1.0; label = 1; }
border e2 ( t = -1.0,  0.0 ) { x = 0.0;  y = t;    label = 1; }
border e3 ( t =  0.0,  1.0 ) { x = t;    y = 0.0;  label = 1; }
border e4 ( t =  0.0,  1.0 ) { x = 1.0;  y = t;    label = 1; }
border e5 ( t =  1.0, -1.0 ) { x = t;    y = 1.0;  label = 1; }
border e6 ( 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;
int n2 = 2 * n;
mesh Th = buildmesh ( e1 ( n ) + e2 ( n ) + e3 ( n ) + e4 ( n ) + e5 ( n2 ) + e6 ( n2 ) );
//
//  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 right hand side function F.
//
func f
//
//  Define boundary condition function G.
//
func g
//
//  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 = 1, fill = true, ps = "test12_u.eps" );
//
//  Plot the mesh.
//
plot ( Th, wait = 1, ps = "test12_mesh.eps" );
//
//  Save the mesh file.
//
savemesh ( Th, "test12.msh" );