//  Discussion:
//
//    Coupled elasticity equations on a square with a slit.
//
//    a * u1xx + b * u1yy + c * u2xy = f1
//    b * u2xx + a * u2yy + c * u1xy = f2
//
//    f1 = 0
//    f2 = 0
//    u1 = g1 on boundary
//    u2 = g2 on boundary
//
//    a = - E * ( 1 - nu^2 ) / ( 1 - 2 * nu )
//    b = - E * ( 1 - nu^2 ) / ( 2 - 2 * nu )
//    c = - E * ( 1 - nu^2 ) / ( 1 - 2 * nu ) / ( 2 - 2 * nu )
//
//    E = 1
//    nu = 0.3
//
//    kappa = 3 - 4 * nu
//    g = 0.5 * e / ( 1 + nu )
//
//    c1 = ( kappa - q * ( lambda + 1 ) * cos ( lambda * theta )
//    c2 = lambda * cos ( ( lambda - 2.0 ) * theta )
//    g1(r,theta) = r^lambda * ( c1 - c2 ) / 2 / g
//    s1 = ( kappa - q * ( lambda + 1 ) * sin ( lambda * theta )
//    s2 = lambda * sin ( ( lambda - 2.0 ) * theta )
//    g2(r,theta) = r^lambda * ( s1 - s2 ) / 2 / g
//
//    Suggested Parameter Values:
//
//    lambda = 0.5444837367825, q =  0.5430755788367
//    lambda = 0.9085291898461, q = -0.2189232362488
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/mitchell_03/mitchell_03a.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 February 2015
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    Mark Ainsworth, Bill Senior,
//    Aspects of an adaptive hp-finite element method:
//    Adaptive strategy for conforming approximation and efficient solvers,
//    Computer Methods in Applied Mechanics and Engineering,
//    Volume 150, 1997, pages 65-87,
//
//    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_03a:\n";
cout << "  FreeFem++ version\n";
cout << "  Mitchell linear elasticity problem.\n";
//
//  Set parameters.
//
real nu = 0.3;
real e = 1.0;
real lambda = 0.5444837367825;
real q = 0.5430755788367;

cout << "\n";
cout << "  Lambda = " << lambda << "\n";
cout << "  Q = " << q << "\n";
//
//  Specify OMEGA, the angle for the slit.
//  We can't use omega = 2pi, or even 2pi - 0.05 because the mesher fails.
//
real omega = 8.0 * pi / 4.0 - 0.10;
cout << "  Omega = " << omega << "\n";
//
//  Specify P, the relative number of points on each line.
//  We can't use P = 3 or the mesher fails.
//
int p = 5;
cout << "  Border precision P = " << p << "\n";
//
//  Specify the lines the form the border of the region.
//
border b1 ( t = 0.0, +1.0 )   { x = t;    y = 0.0; label = 1; }
real l1 = 1.0;
int n1 = ( round ) ( l1 * p + 1 );
border b2 ( t = 0.0, +1.0 )   { x = +1.0; y = t;   label = 1; }
real l2 = 1.0;
int n2 = ( round ) ( l2 * p + 1 );
border b3 ( t = 1.0, -1.0 )   { x = t;    y = 1.0; label = 1; }
real l3 = 2.0;
int n3 = ( round ) ( l3 * p + 1 );
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, +1.0  ) { x = t; y = -1.0; label = 1; }
real l5 = 2.0;
int n5 = ( round ) ( l5 * p + 1 );
real xc = 1.0;
real yc = tan ( omega );
real rc = 1.0 / cos ( omega );
border b6 ( t = -1.0, yc  ) { x = 1.0; y = t; label = 1; }
real l6 = yc + 1.0;
int n6 = ( round ) ( l6 * 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 );
//
//  Create the mesh.
//
mesh Th = buildmesh ( b1 ( n1 ) + b2 ( n2 ) + b3 ( n3 ) + b4 ( n4 ) 
                    + b5 ( n5 ) + b6 ( n6 ) + b7 ( n7 ) );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "mitchell_03a_mesh.eps" );
//
//  Save the mesh file.
//
savemesh ( Th, "mitchell_03a.msh" );
//
//  Define Vh, the finite element space defined over Th, using 
//  a vector of two P1 basis functions.
//
fespace Vh ( Th, [ P1, P1 ] );
//
//  Define U and V, piecewise continuous functions over Th.
//
Vh [ u1, u2 ];
Vh [ v1, v2 ];
//
//  Define the right hand side function F.
//
func f1 = 0.0;
func f2 = 0.0;
//
//  Define the boundary condition function G.
//
//  While I would have MUCH PREFERRED to write G1 and G2 as standard C++
//  functions, I could not discover the appropriate way to do so, and so
//  I had to pack the whole formula for each into a one-liner.
//
real kappa = 3.0 - 4.0 * nu;
real g = 0.5 * e / ( 1.0 + nu );

func g1 = pow ( sqrt ( x * x + y * y ), lambda ) *
  ( ( kappa - q * ( lambda + 1.0 ) ) * cos ( lambda * atan2 ( y, x ) ) 
  - lambda * cos ( ( lambda - 2.0 ) * atan2 ( y, x ) ) ) / 2.0 / g;

func g2 = pow ( sqrt ( x * x + y * y ), lambda ) *
  ( ( kappa - q * ( lambda + 1.0 ) ) * sin ( lambda * atan2 ( y, x ) ) 
  - lambda * sin ( ( lambda - 2.0 ) * atan2 ( y, x ) ) ) / 2.0 / g;
//
//  Set coefficients for the variational problem.
//
real a = + e * ( 1.0 - nu * nu ) / ( 1.0 - 2.0 * nu );
real b = + e * ( 1.0 - nu * nu ) / ( 2.0 - 2.0 * nu );
real c = + e * ( 1.0 - nu * nu ) / ( 1.0 - 2.0 * nu ) / ( 2.0 - 2.0 * nu );
//
//  Solve the variational problem.
//  u1 = x, u2 = y is accepted as a boundary condition,
//  but not u1 = g1, u2 = g2.
//
solve Laplace ( u1, u2, v1, v2 ) 
  = int2d ( Th ) ( a*dx(u1)*dx(v1) + b*dy(u1)*dy(v1) + c*dx(u2)*dy(v1) )
  + int2d ( Th ) ( b*dx(u2)*dx(v2) + a*dy(u2)*dy(v2) + c*dx(u1)*dy(v2) )
  + on ( 1, u1 = g1, u2 = g2 );
//
//  Plot the solution.
//
plot ( [ u1, u2 ], wait = true, fill = true, ps = "mitchell_03a_u.eps" );
//
//  Terminate.
//
cout << "\n";
cout << "mitchell_03a:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );