Expression degree and quadrature rules

Question

I'm trying to understand what happens when we specify the flag degree
in an Expression. In particular, I'd like to understand what impact it has
on convergence and how to interpret it in terms of quadrature rules.

For this purpose, I've experimented with the following code, which
estimates FE convergence in $H^1$. When constructing the expression of rhs,
it is sufficient to set its degree to be $\geq 2$ to maintain cubic convergence.
On the other hand, the order of convergence decreases to 2 when degree
is set it to 0 or 1.

Bonus question: how to manually impose the order of the quadrature rule to
be employed when assembling the stiffness matrix and the rhs in the variational
form?

from dolfin import *
import numpy as np
import sympy as sy
from sympy.printing import ccode

#build reference solution
x_ = sy.Symbol('x[0]')
y_ = sy.Symbol('x[1]')
uref_ = sy.sin(x_ + 2*y_)
rhs_ = - sy.diff(uref_, x_, x_) - sy.diff(uref_, y_, y_)

uref = Expression(ccode(uref_), degree=10)
rhs = Expression(ccode(rhs_), degree=10)

err_h1 = []
iMax = 4
nElem = [20*2**ii for ii in range(iMax)]
for N in nElem:
    mesh = UnitSquareMesh(N, N)

    U = FunctionSpace(mesh, "CG", 3)

    u = Function(U)
    v = TestFunction(U)

    F = inner(grad(v), grad(u))*dx - rhs*v*dx
    bc = DirichletBC(U, uref, "on_boundary")
    solve(F == 0, u, bc)
    err_h1.append(errornorm(uref, u, 'H1'))

print "H1 error: ", err_h1

#estimate rate of convergence
errh1 = np.array(err_h1) 
quotient = errh1[:-1]/errh1[1:]
rateh1 = np.log(quotient)/np.log(2)
print "rate: ", rateh1

Answer 1

An Expression will be interpolated into a finite element space before the quadrature is carried out. The degree argument to the Expression constructor is the degree of the FE space. Alternatively, you can specify the element, e.g.

P1 = FiniteElement("Lagrange", triangle, 1)
f = Expression("sin(x[0])", element = P1)

The quadrature degree is by default automatically determined by dolfin, such that the quadrature is exact for highest order polynomial in the form being assembled. You can manually set the quadrature degree by passing additional arguments to assemble. For example,

A = assemble(form, form_compiler_options = {"quadrature_degree", 2})

Also note that if the Expression is replaced with a ufl function of the coordinates it will not be interpolated. For example, you can replace the code for your reference solution with the following.

x = SpatialCoordinate(mesh)
uref = sin(x[0] + 2 * x[1])
rhs = -div(grad(uref))

Comments

I noticed that the name of the flag should probably be

form_compiler_parameters

and the syntax should be

form_compiler_parameters = {"quadrature_degree": 2}

Thank you once again