/****************************************************************************/
/* This file is part of FreeFEM.                                            */
/*                                                                          */
/* FreeFEM is free software: you can redistribute it and/or modify          */
/* it under the terms of the GNU Lesser General Public License as           */
/* published by the Free Software Foundation, either version 3 of           */
/* the License, or (at your option) any later version.                      */
/*                                                                          */
/* FreeFEM is distributed in the hope that it will be useful,               */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of           */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the            */
/* GNU Lesser General Public License for more details.                      */
/*                                                                          */
/* You should have received a copy of the GNU Lesser General Public License */
/* along with FreeFEM. If not, see <http://www.gnu.org/licenses/>.          */
/****************************************************************************/

verbosity=0;

// Parameters
real theta = 4.*pi/3.;
real a = 1., b = 1.; // the length of the semimajor axis and semiminor axis
func f = -4*(cos(x^2+y^2-1) - (x^2+y^2)*sin(x^2+y^2-1));

// Mesh
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);}

// Exact solution
func phiexact = sin(x^2+y^2-1);

// Solve loop
real[int] L2error(2);
for(int n = 0; n < 2; n++) {
  mesh Th = buildmesh(Gamma1(40*(n+1)) + Gamma2(20*(n+1)));
  fespace Vh(Th, P2);
   Vh phi, w;

  solve laplace(phi, w)
    =int2d(Th)(
        dx(phi)*dx(w)
      + dy(phi)*dy(w)
    )
    - int2d(Th)(f*w)
    - int1d(Th,Gamma2)(2*w)
    + on(Gamma1, phi=0);
  plot(Th, phi, wait=true, ps="membrane.eps"); //Plot Th and phi

  L2error[n] = sqrt(int2d(Th)((phi - phiexact)^2));
}

for(int n = 0; n < 2; n++)
  cout << " L2error " << n << " = " <<  L2error[n] << endl;

cout << " convergence rate = " << log(L2error[0]/L2error[1])/log(2.) << endl;
