# Interpolation of low-order polynomial not so accurate

I found that in dolfinx, the following two equivalent expressions

\begin{align} \int_\Omega div(\nabla v) v, \qquad \text{and} \qquad - \int_{ \Omega} \nabla v \cdot \nabla v + \int_{\partial \Omega} v\nabla v \cdot n \end{align}

are numerically not the same, i.e. their difference is quite “large”. The following assumes that v is a second order polynomial and computes the two integrals:

from mpi4py import MPI
import ufl
from basix.ufl import element
from dolfinx.fem import (Function, form, functionspace,assemble_scalar)
from dolfinx.mesh import *
from ufl import div, dx, ds, grad, inner, FacetNormal

# mesh, spaces
Nx=50
msh = create_unit_square(MPI.COMM_WORLD, Nx, Nx)
P1 = element("Lagrange", "triangle", 3)
XV = functionspace(msh, P1)
x = ufl.SpatialCoordinate(msh)
n = FacetNormal(msh)

# function 1+x^2*y^2
def func(x):
return 1. + x[0]**2 * x[1]**2
pex = Function(XV)
pex.interpolate(func)

# INTEGRAL BEFORE AND AFTER INT. BY PARTS
print(assemble_scalar(I1)-assemble_scalar(I2))


On my installation (dolfinx version 0.7, conda-forge, virtual environment in linux), using third-order FEM polynomials, the difference is around 10^{-4}. Is this to be expected?

Bonus question: since when is the double dollar sign syntax  for mathmode not working anymore? For some reason, in this post it was necessary to use single dollar signs for the mathjax compiler.

1 Like

Your function v \in V is not satisfied by the regularity requirement of the polynomial space where you’ve chosen v = 1 + x^2 y^2 \notin \mathcal{P}^3. You can either ensure that v \in V by making p \geq 4 or simplify the function, e.g., v = x^2 + y^2. Correcting for this I see equivalence to machine precision.

It’s more straightforward to consider your problem in terms of the question: what are the regularity requirements of \sigma such that

\int_\Omega \nabla \cdot \sigma \; v \; \mathrm{d} \vec{x} = - \int_\Omega \sigma \cdot \nabla v \; \mathrm{d} \vec{x} + \int_{\partial\Omega} \sigma \cdot \vec{n} \; v \; \mathrm{d} s

holds?

This has lead to the development of finite element bases which are H(\mathrm{div}) conforming. See, for example the mixed Poisson demo in legacy dolfin.

You can investigate this further in terms of the research going into divergence free (thereby H(\mathrm{div}) conforming) incompressible flow velocity approximations.

See for example:

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 4, 4)
n = ufl.FacetNormal(mesh)

def assemble_system(sigma):
# arbitrary test function
v = dolfinx.fem.Function(dolfinx.fem.FunctionSpace(mesh, ("DG", 2)))
v.interpolate(lambda x: x[0])

lhs = ufl.inner(ufl.div(sigma), v) * ufl.dx
rhs = - ufl.inner(sigma, ufl.grad(v)) * ufl.dx + ufl.inner(ufl.dot(sigma, n), v) * ufl.ds

lhs = dolfinx.fem.assemble_scalar(dolfinx.fem.form(lhs))
rhs = dolfinx.fem.assemble_scalar(dolfinx.fem.form(rhs))
return lhs - rhs

def try_scalar_element(family, p):
V = dolfinx.fem.FunctionSpace(mesh, (family, p))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: 1.0 + x[0]**2 * x[1]**2)
residual = assemble_system(sigma)
print(f"element: {V.ufl_element()}, residual = {residual:.3e}")
try_scalar_element("CG", 1)
try_scalar_element("CG", 2)
try_scalar_element("CG", 3)
try_scalar_element("CG", 4)

def try_vector_element(family, p):
V = dolfinx.fem.FunctionSpace(mesh, (family, p))
sigma = dolfinx.fem.Function(V)
sigma.interpolate(lambda x: 1 + np.stack((x[0]**2*x[1]**2, x[0]*x[1]**2)))
residual = assemble_system(sigma)
print(f"element {V.ufl_element()}: residual = {residual:.3e}")
try_vector_element("RT", 1)
try_vector_element("BDM", 1)


Which gives me:

element: Basix element (P, triangle, 1, gll_warped, unset, False), residual = -2.500e-01
element: Basix element (P, triangle, 2, gll_warped, unset, False), residual = 4.036e-02
element: Basix element (P, triangle, 3, gll_warped, unset, False), residual = -7.986e-03
element: Basix element (P, triangle, 4, gll_warped, unset, False), residual = -1.522e-13
element: Basix element (RT, triangle, 1, legendre, unset, False): residual = -1.221e-15
element: Basix element (BDM, triangle, 1, legendre, unset, False): residual = -1.887e-15

1 Like

Thank you very much for you answer!