int m = 5; real[int] L2error(m+1); real[int] H1error(m+1); load "Element_P3"; load "Element_P4"; int n = 1; func f = -(2*n+1)^2*(x^2+y^2)^(n-0.5); func uexact = sqrt(x^2+y^2)^(2*n+1); // in H^(2*n+1) func dxuexact = (2*n+1)*(x^2+y^2)^(n-0.5)*x; func dyuexact = (2*n+1)*(x^2+y^2)^(n-0.5)*y; int n1 = 4; // number of intervals for each edge of the square for (int i=0;i<=m;i++) { mesh Th = square(n1,n1,[-1+2*x,-1+2*y]); cout << "number of sites: " << Th.nv << endl; cout << "number of triangles: " << Th.nt << endl; fespace Vh(Th,P1); int r = 1; int qo = 6; //fespace Vh(Th,P2); int r = 2; int qo = 6; //fespace Vh(Th,P3); int r = 3; int qo = 6; //fespace Vh(Th,P4); int r = 4; int qo = 6; // int qo = 8; Vh u,phi; cout << "degree of freedom: " << Vh.ndof << endl; solve Poisson(u,phi)=int2d(Th,qforder=qo)(dx(u)*dx(phi)+dy(u)*dy(phi))-int2d(Th,qforder=qo)(f*phi) +on(1,u=sqrt(x^2+1)^(2*n+1))+on(2,u=sqrt(1+y^2)^(2*n+1)) +on(3,u=sqrt(x^2+1)^(2*n+1))+on(4,u=sqrt(1+y^2)^(2*n+1)); real L22 = int2d(Th)((u-uexact)^2); real H12 = int2d(Th)((dx(u)-dxuexact)^2+(dy(u)-dyuexact)^2); H1error[i] = sqrt(L22+H12); // order s, s = min(r,2n) L2error[i] = sqrt(L22); // order s+1, s = min(r,2n) n1 = n1*2; } for (int i=0;i