Interpolation of low-order polynomial not so accurate

dolfinx
1

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
I1 = form(inner(div(grad(pex)),pex)*dx)
I2 = form((-1)* inner(grad(pex),grad(pex))*dx + inner(grad(pex),pex*n)*ds )
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
4

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)
	sigma = ufl.grad(u)
	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
5

Thank you very much for you answer!