Can't reproduce integration of a function using P1 finite elements.

Question

Hello,

I am trying to debug a domain decomposition (DD) code that uses FEniCS to assemble the subdomain problems. The DD method I have coded converges at a different rate and with a different optimality parameter than the same algorithm in FreeFEM++ and my own finite element implementation.

To get to the bottom of this, I have stored the assembled right hand side vector $\int f \cdot v$ in all three cases. I was able to get the FreeFEM++ code to match my attempt at the integration using one-point Gaussian quadrature, but I cannot get the FEniCS code to match the other two. Below is the FEniCS code forcing the particular quadrature rule. Am I doing something wrong?

from dolfin import *

parameters["form_compiler"]["representation"] = "quadrature"
parameters["form_compiler"]["quadrature_degree"] = 2

n = 9 
x = np.linspace(0,1,n)
y = np.linspace(0,1,n)

mesh = UnitSquareMesh(n-1,n-1)

f = Expression('(-8.*pi*pi) * ( sin((2.0*pi*x[0]) - (0.75*pi)) * \
                 sin((2.0*pi*x[1]) - (0.75*pi)) )')

# Function spaces
V = FunctionSpace(mesh, 'CG', 1)
v2d = vertex_to_dof_map(V)
v = TestFunction(V)
L = assemble(inner(f,v)*dx)

Answer 1

FEniCS will, by default, interpolate the Expression f in a finite element space. Try specifying the element to control the accuracy, e.g.

f = Expression('(-8.*pi*pi) * ( sin((2.0*pi*x[0]) - (0.75*pi)) * \
             sin((2.0*pi*x[1]) - (0.75*pi)) )', degree=4)

to interpolate f with a polynomial of degree 4.

The present default behaviour is not ideal. See discussion at https://bitbucket.org/fenics-project/dolfin/issues/355

Comments

Thank you for this! I was able to match the integrations by not only using degree=6, but by setting the global parameters:

parameters["form_compiler"]["representation"] = "quadrature"
parameters["form_compiler"]["quadrature_degree"] = 1

This was a surprise, because I was comparing to a qforder=2 in FreeFEM. I guess the difference is that in FEniCS the "quadrature_degree" refers to the number of points used in the scheme and in FreeFEM the "qforder" actually refers to the order of the scheme.

This might be something worth noting in the documentation (unless I missed it). FreeFEM docs have a table (p185) explaining the different choices.