Why is the projection of the gradient of a constant function with quadratic elements nonzero?

Question

I am trying to compute the gradient of a function defined on a CG function space of degree 2. When the function is constant, I am finding that the result is nonzero. I don't understand why. Here is an example code snippet:

from dolfin import *

mesh = UnitSquareMesh( 16, 16 )

V = FunctionSpace(mesh, "CG", 2)
v = Function(V)

v.vector()[:] = 10.

w = project(-grad( v ), VectorFunctionSpace(mesh, "CG", 2) )

print w.vector()[:].min(), w.vector()[:].max()

The values of w are small, but nonzero. When I change V to be of degree 1, all values of w are exactly zero. Could someone explain this discrepancy? Any suggestions as to how I could keep V with degree 2 and get the correct gradient values?

Answer 1

This is a numerical roundoff in the generated code (for the projection forms). Actually it is not that bad. One should maybe measure the error as w.vector()[:].norm('linf')/v.vector().norm('linf'). One can also tweak form compiler parameters to get get some improvement.

>>> w = project(-grad( v ), VectorFunctionSpace(mesh, "CG", 2),  form_compiler_parameters={})
>>> print w.vector()[:].norm('linf')/v.vector().norm('linf')
3.71706040551e-13
>>> w = project(-grad( v ), VectorFunctionSpace(mesh, "CG", 2), form_compiler_parameters={'precision': 17})
>>> print w.vector()[:].norm('linf')/v.vector().norm('linf')
3.30199902966e-13
>>> w = project(-grad( v ), VectorFunctionSpace(mesh, "CG", 2), form_compiler_parameters={'precision': 17, 'epsilon': 1e-20})
>>> print w.vector()[:].norm('linf')/v.vector().norm('linf')
1.3903806388e-13
>>> w = project(-grad( v ), VectorFunctionSpace(mesh, "CG", 2), form_compiler_parameters={'precision': 17, 'epsilon': 1e-20, 'representation': 'uflacs'})
>>> print w.vector()[:].norm('linf')/v.vector().norm('linf')
1.88949719446e-14