// membrane_error.edp
//
//  Discussion:
//
//    Solve the membrane problem using a manufactured solution.
//
//    We use a circular domain C with boundary Gamma1 + Gamma2.
//
//    We have 
//      u = sin ( x^2 + y^2 - 1 )
//    for which we need to specify a nonhomogeneous Neumann condition on Gamma2.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/membrane_error/membrane_error.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    28 July 2015
//
//  Reference:
//
//    Frederic Hecht,
//    Freefem++,
//    Third Edition, version 3.22
//
cout << "\n";
cout << "membrane_error:\n";
cout << "  FreeFem++ version.\n";
//
//  Define Gamma1 and Gamma2, the boundaries.
//  A and B are the lengths of the semimajor and semiminor axes:
//  THETA specifies the angle at which we switch from GAMMA1 to GAMMA2.
//
real a = 1.0;
real b = 1.0;
real theta = 4.0 * pi / 3.0;
border Gamma1 ( t = 0, theta )    { x = a * cos(t); y = b*sin(t); }
border Gamma2 ( t = theta, 2 * pi ) { x = a * cos(t); y = b*sin(t); }
//
// uexact: the exact solution.
//
func uexact = sin ( x^2 + y^2 - 1.0 );
//
//  Make space for the errors.
//
real[int] L2error(2);

int n;

for ( n = 0; n < 2; n++ )
{
//
//  Define Th, the triangulation of the "left" side of the boundaries.
//
  mesh Th = buildmesh ( Gamma1(20*(n+1)) + Gamma2(10*(n+1)) );
//
//  Define Vh, the finite element space defined over Th, using P2 basis functions.
//
  fespace Vh ( Th, P2 );
//
//  Define U and V, piecewise P2 continuous functions over Th.
//
  Vh u;
  Vh v;
//
//  Define F, the right hand side.
//
  func f = - 4.0 * ( cos ( x^2 + y^2 - 1.0 ) 
    - ( x^2 + y^2 ) * sin ( x^2 + y^2 - 1.0 ) );
//
//  Solve the problem.
//
  solve Laplace ( u, v ) 
    = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
    - int2d ( Th ) ( f * v ) 
    - int1d ( Th, Gamma2 ) ( 2.0 * v )
    + on ( Gamma1, u = 0.0 );
//
//  Plot the solution.
//
  plot ( u, fill = true, wait = true, ps = "membrane_error_u_"+n+".ps" );
//
//  Plot the mesh.
//
  plot ( Th, wait = true, ps = "membrane_error_mesh_"+n+".ps" );
//
//  Compute the error.
//
  L2error[n] = sqrt ( int2d ( Th ) ( ( u - uexact )^2 ) );
  cout << "  L2error " << n << " = " << L2error[n] << endl;
}

cout << "  Convergence rate = " 
     << log ( L2error[0] / L2error[1] ) / log ( 2.0 ) << endl;
//
//  Terminate.
//
cout << "\n";
cout << "membrane_error:\n";
cout << "  FreeFem++ version.\n";

exit ( 0 );