CG2 for poisson equation should be 3rd order L2 convergence but for some reason I am getting near 4th order convergence. If I use DG2 elements, I get the 3rd order convergence as theory suggests.

Why is this happening? Below is the snippet of my code

mesh = UnitSquareMesh(20,20,quadrilateral=False)

Q = FunctionSpace(mesh,CG,2)


...


x = SpatialCoordinate(Q.mesh())

u_exact = interpolate(sin(2*x[0]**2*math.pi)*sin(2*x[1]**2*math.pi),Q)

f = interpolate( 16*math.pi**2*sin(2*math.pi*x[0]**2)*sin(2*math.pi*x[1]**2)*(x[0]**2+x[1]**2) 

      - 4*math.pi*(cos(2*math.pi*x[0]**2)*sin(2*math.pi*x[1]**2) + sin(2*math.pi*x[0]**2)*cos(2*math.pi*x[1]**2)),Q)


...


L2_error = errornorm(u_exact,u_h)


I never had this "problem" earlier this year when I ran CG2 experiments (though I didn't use the UFL format, but I don't think that's the problem here).

Justin