Incorrect convergence rate for "Element_P4" in 2d

General Discussion
1 
load "Element_P3"
load "Element_P4"
int  k  = 1*2*2;
mesh Th = square(10*k, 10*k);
//plot(Th);
func Pk = P4;
fespace Vhc ( Th, Pk );
Vhc uc,ucX,ucY,fC,fCh,ucEX,err;
//--exact solution, partial derivatives and rhs----
ucEX = exp(x^2/4+y^2/4);
ucX  = exp(x^2/4+y^2/4)*x/2;
ucY  = exp(x^2/4+y^2/4)*y/2;
fC   = exp(x^2/4+y^2/4)*(-x^2/4-y^2/4);
//------ get the matrice and right hand side ------
varf ahc(u,v) = int2d(Th)( dx(u)*dx(v)+dy(u)*dy(v) +u*v );
varf rhsC(unused,v) = int2d(Th)(fC*v)
                     +int1d(Th,1)( (N.x*ucX+N.y*ucY)*v )
                     +int1d(Th,2)( (N.x*ucX+N.y*ucY)*v )
                     +int1d(Th,3)( (N.x*ucX+N.y*ucY)*v )
                     +int1d(Th,4)( (N.x*ucX+N.y*ucY)*v ); 
//----------------solve the linear system----------
matrix Ahc;
Ahc = ahc(Vhc,Vhc);
set(Ahc, solver=CG, eps=1e-20);
fCh[] = rhsC(0,Vhc);
uc[]  = Ahc^-1*fCh[];
//--------------- get the error -------------------
err[] = uc[]-ucEX[];
cout.scientific<<"Infinity error:"
<<max(abs(err[].min), abs(err[].max))<<endl;
real l2er;
l2er  = ( int2d(Th)( (ucEX-uc)^2) );
cout.scientific<<"L2 error: "<<sqrt(l2er)<<endl;

Question: in this example, I verified the L2 convergence rate for P3 element as k varies which is correct.
However, the convergence rate for P4 element is incorrect. Am I wrong somewhere?

Thanks a lot for your help.
P4Element.edp (1.1 KB)Preformatted text

2

You need to use an appropriate qforder, see Quadrature formulae or this similar topic Default quadrature rule.

3

Thank you so much, Pierre.
It works!

4

Dear Newhat,

may I ask, which qforder did you use?

5

qforder=5, the quadrature formulae should have higher (at least not less) order than the element.