real pi=4*atan(1.);
real alpha = 7*pi/4;     
real omega = pi/alpha;
 
int N=8;     // Le premier pas du maillage est h=1/N
int Nit = 5; // Nbre de divisions par 2 du pas de maillage (a augmenter au maximum des possibilites de l'ordinateur)

//Definition de la frontiere du domaine, label=0 correspond a la condition de Neumann, label=1 a la condition de Dirichlet.
border D1(t=0,1){x=t; y=0; label=1;};
border D2(t=1,0){x=t*cos(alpha); y=t*sin(alpha); label=1;};
border C(t=0 , alpha){x=cos(t); y=sin(t); label=0;};

real[int] errH10(Nit), errL2(Nit), ordredeconvergence(Nit - 1), tabh(Nit);  // tableaux d'erreur ,d'ordre de convergence, de pas de maillage

// Definitions de fonctions utiles : r,theta (coordonnees polaires), phi et ses derivees, u0 et ses derivees et f.
func r=sqrt(x^2 + y^2) ;
func phi = (r<0.5) 
             + (r>=0.5)*0.5*(1 + cos(2*pi*(r - 0.5)));
func phiprime = -pi*(r>=0.5)*sin(2*pi*(r - 0.5));
func phiseconde =  (-2*pi^2)*(r>=0.5)*cos(2*pi*(r - 0.5));

func rnot0 = r+exp(-100000*r);                      // Pour eviter les divisions par 0 en r=0 dans le calcul de theta.
func theta = acos(x/rnot0) * (y>=0) + (2*pi - acos(x/rnot0))*(y<0);

func u0 = phi*r^omega * sin(omega*theta);
func GRADru0 = r^(omega - 1) *(phiprime*r + omega*phi)*sin(omega*theta);
func GRADthetau0 = (omega/rnot0)*phi*r^omega*cos(omega*theta);
func dxu0 = cos(theta)*GRADru0 - sin(theta)*GRADthetau0;
func dyu0 = sin(theta)*GRADru0 +cos(theta)*GRADthetau0;

func f = -r^(omega - 1)*((1 + 2*omega)*phiprime + r*phiseconde)*sin(omega*theta);

//////////////////////////////   BOUCLE SUR LE  PAS DE MAILLAGE  ///////
for(int i=0; i<Nit; i++ ){

//////////// Creation du maillage Th + Espace d'Elements Finis Vh //////////
mesh Th=buildmesh( C(alpha*N)  +D1(N) + D2(N));  //plot(Th, wait = 1);
fespace Vh(Th,P1);     Vh uh,vh,err;

////// Definition et resolution de la Formulation Variationnelle //
solve Lap(uh,vh) = int2d(Th)(dx(uh)*dx(vh) + dy(uh)*dy(vh) ) + int2d(Th)(-f*vh) + on(1,uh=0);
plot(uh,wait=1,fill=1,value=1,grey=1, cmm="SOLUTION APPROCHEE",boundary=0);

////////////// Calculs d'erreurs //////
err= uh - u0;   plot(err,wait=1,fill=1,value=1,bw=1, cmm="ERREUR", boundary=1);
errL2[i] = sqrt(int2d(Th)( (uh - u0)^2));
errH10[i] = sqrt(int2d(Th)((dx(uh) - dxu0)^2 + (dy(uh) - dyu0)^2 ));
tabh[i] =1./N;

N = 2*N;  
}
// //////////////////////////// FIN  BOUCLE SUR LE  PAS DE MAILLAGE  //////////

////////////// Estimations de l'ordre de convergence  ///////
for(int i=0; i<Nit - 1; i++) {ordredeconvergence[i] = log(errL2(i)/errL2(i + 1))/(log(2.));}
//////////////   AFFICHAGE ////////////////////////
cout  << endl << "pas des maillages : " << tabh << endl;
cout << endl << "erreurs H10 : " << errH10(0:Nit - 1) << endl;
cout << endl << "erreurs L2 : " << errL2(0:Nit - 1) << endl;
cout << endl << "estimations de l'ordre : " << ordredeconvergence(0:Nit - 2) << endl << endl;