Incorrect convergence rate for "Element_P4" in 2d

load "Element_P3"
load "Element_P4"
int  k  = 1*2*2;
mesh Th = square(10*k, 10*k);
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

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

Thank you so much, Pierre.
It works!