//  Discussion:
//
//    The Laplacian operator is applied on a square with a reentrant corner.
//
//    The geometry is defined by an internal angle PI < OMEGA <= 2PI.
//    
//    -uxx - uyy = f in the region.
//    u = g on the boundary.
//
//    F(X,Y) = 0
//    ALPHA = PI / OMEGA
//    R = sqrt ( X^2 + Y^2 )
//    THETA = arctan ( Y / X )
//    G(X,Y) = R^ALPHA * SIN ( ALPHA * THETA)
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_02/mitchell_02d.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    24 May 2020
//
//  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_02d\n";
cout << "  FreeFem++ version\n";
cout << "  Reentrant corner\n";

real omega = 7.0 * pi / 4.0;

cout << "  Omega = " << omega << "\n";
int p = 5;

border b1 ( t = 0.0, +1.0 )   { x = t;    y = 0.0; label = 1; }
real l1 = 1.0;
int n1 = ( round ) ( l1 * p + 1 );
cout << "  N1 = " << n1 << "\n";
border b2 ( t = 0.0, +1.0 )   { x = +1.0; y = t;   label = 1; }
real l2 = 1.0;
int n2 = ( round ) ( l2 * p + 1 );
cout << "  N2 = " << n2 << "\n";
border b3 ( t = 1.0, -1.0 )   { x = t;    y = 1.0; label = 1; }
real l3 = 2.0;
int n3 = ( round ) ( l3 * p + 1 );
cout << "  N3 = " << n3 << "\n";
//
//  What we do depends on OMEGA!
//
real xc = - cos ( omega ) / sin ( omega );
real yc = -1.0;
real rc = -1.0 / sin ( omega );

border b4 ( t = 1.0, -1.0 ) { x = -1.0; y = t;   label = 1; }
real l4 = 2.0;
int n4 = ( round ) ( l4 * p + 1 );
border b5 ( t = -1.0, xc  ) { x = t; y = yc; label = 1; }
real l5 = xc + 1.0;
int n5 = ( round ) ( l5 * p + 1 );
border b7 ( t = rc, 0.0 ) { x = t * cos ( omega ); y = t * sin ( omega ); label = 1; }
real l7 = rc;
int n7 = ( round ) ( l7 * p + 1 );

mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) 
                    + b5 ( n5 ) + b7 ( n7 ) ); 
//
//  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 right hand side function F.
//
func f = 0.0;
//
//  Define the boundary condition function G.
//
real alpha = pi / omega;
func g = pow ( x * x + y * y, alpha / 2.0 ) * sin ( alpha * atan2 ( y, x ) );
//
//  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_02d_u.eps" );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "mitchell_02d_mesh.eps" );
//
//  Save the mesh file.
//
savemesh ( Th, "mitchell_02d.msh" );
//
//  Terminate.
//
cout << "\n";
cout << "mitchell_02d\n";
cout << "  Normal end of execution.\n";

exit ( 0 );