//  Discussion:
//
//    Let an ellipse have semimajor axis 2 and semiminor axis 1:
//      (x/2)^2 + (y/1)^2 = 1.
//    Let Gamma1 be the ellipse boundary from 0 to 4pi/3, and
//    Gamm2 the boundary from 4pi/3 to 2pi.
//
//    Solve:
//      -uxx - uyy = f in Gamma1 + Gamma2,
//      u = x on Gamma1.
//      du/dn = 0, natural boundary condition on Gamma2.
//
//  Location:
//
//    http://people.sc.fsu.edu/~jburkardt/freefem_src/membrane/membrane.edp
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    18 June 2015
//
//  Reference:
//
//    Frederic Hecht,
//    Freefem++,
//    Third Edition, version 3.22
//
cout << "\n";
cout << "membrane:\n";
cout << "  freefem++ version\n";
cout << "  Model the deflection of an elastic membrane inside an ellipse..\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 = 2.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); }
//
//  Th: the triangulation of the "left" side of the boundaries.
//
mesh Th = buildmesh ( Gamma1(100) + Gamma2(50) );
//
//  Plot the mesh.
//
plot ( Th, wait = true, ps = "membrane_mesh.ps", cmm = "membrane mesh" );
//
//  Define Vh, the finite element space defined over Th, using P2 basis functions.
//
fespace Vh ( Th, P2 );
//
//  Define U, V, and F, piecewise P2 continuous functions over Th.
//  F is a constant function.
//
Vh u;
Vh v;
Vh f = 1.0;
//
//  Define G, a function for the Dirichlet boundary conditions.
//
func g = x;
//
//  Define the weak form of the PDE.
//
problem poisson ( u, v ) 
  = int2d ( Th ) ( dx(u)*dx(v) + dy(u)*dy(v) )
  - int2d ( Th ) ( f * v ) 
  + on ( Gamma1, u = g );
//
//  Solve the PDE.
//
poisson;
//
//  Plot the solution.
//
plot ( u, fill = true, wait = true, ps = "membrane_u.ps", cmm = "membrane solution" );
//
//  Save the mesh file.
//
savemesh ( Th, "membrane.msh" );
//
//  Terminate.
//
cout << "\n";
cout << "MEMBRANE:\n";
cout << "  Normal end of execution.\n";

exit ( 0 );