Wrong interior facet integral value BDM element

dolfinx
1

Today I tried integrating BDM elements of order 1 https://defelement.org/elements/brezzi-douglas-marini.html along interior edges of a 2D domain and encountered a seemingly wrong value.
In the finite element space there are two degrees of freedom on each facet, so I expect the integral of a BDM-function with its outwards normal to be the average between the two DOF values. However my code only returns the value of one singular DOF, which (in my opinion) should not be the case.
MWE:

import numpy as np
from dolfinx.mesh import create_unit_square
from dolfinx.fem import functionspace, Function, assemble_scalar
from mpi4py import MPI
import ufl

mesh = create_unit_square(MPI.COMM_WORLD, 1, 1)
V = functionspace(mesh, ("BDM", 1))

bdm_function = Function(V)
rand_init = np.random.rand(len(bdm_function.x.array))
bdm_function.x.array[:] = rand_init

n = ufl.FacetNormal(mesh)
dS = ufl.Measure("dS", domain=mesh) 

form_final = ufl.inner(bdm_function('+'), n('+')) * dS
result = assemble_scalar(fem.form(form_final))

print("Rand init: ", rand_init)
print("Computed integral:", result)
print("DOF value at index 2: ", rand_init[2])

On my end the result is always equal to this one parameter at DOF 2. Changing quadrature rules or orders didn’t impact the result. The BDM space has the expected dimension, 10 and I expected that the integral uses the 2 DOFs on the interior edge instead of only one.
I tried the (presumably) same example in legacy fenics, where it returned the expected value:

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "BDM", 1)
bdm_function = Function(V)
random_vector = np.random.rand(len(bdm_function.vector()[:]))
bdm_function.vector()[:] = random_vector

n = FacetNormal(mesh)
dS = Measure("dS", domain=mesh) 

form = inner(bdm_function('-'), n('-')) * dS
result = assemble(form)

print("Computed integral: ", result)
print("Predicted result: ", (random_vector[2] + random_vector[3])/2)

I got the indices by using np.arange before, changing directions of the normal as expected only flips the sign of the result.

These two code snippets were done in disjoint conda envivorments, for the first I installed version 0.9.0 of FEniCSx from conda-forge and it uses fenics-ufl version 2024.2.0 with python 3.12.9. For the second I also installed FEniCS from conda-forge, so everything has version 2019.1.0 , again in python, version 3.10.16

2

Hello,

The BDM elements in Basix are actually a variant of those shown on DefElement: on DefElement, the DOFs are defined as integrals against 1-s and s, but in Basix they are defines as integrals against the Legendre polynomials, ie 1 and \sqrt3(2s-1) for lowest degree BDM. This is so that these elements are more stable when you go to high degree.

This change means that the integral of the basis functions over the edge are all equal to zero except for the first basis function associated with that edge. If you were to instead multiply by a non-constant polynomial then integrate, you should see some different non-zero values on the the edges.

In case it’s helpful, the basis functions for degree 1 BDM in Basix come out as:

  • (-x, -y)
  • \sqrt3(x, -y)
  • (x - 1, y)
  • \sqrt3(1-x - 2y, y)
  • (-x, 1 - y)
  • \sqrt3(-x, 2x + y - 1)

I computed these using this code using symfem:

import symfem

element = symfem.create_element("triangle", "BDM", 1, variant="legendre")
print(element.get_basis_functions())
3 Likes
3

Hello,
thanks for the quick reply on Friday. You are of course correct and that fully explains my confusion in the question, I didn’t read the “Variants” section on DefElement .
Unfortunately I need something to represent \mathcal{DG}_1 on the interior facets of my mesh and were therefore interested in the other variant with integrals taken against the Lagrange basis functions.
In case anyone else is interested, I defined them manually with a basix.ufl.custom_element. You can find the code below. It seems to work for my needs with the versions stated in my original question, but I of course cannot guarantee correctness.

#custom BDM element
import basix
import basix.ufl
from dolfinx import default_real_type

wcoeffs = np.eye(6)#*0.5

pts, wts = basix.make_quadrature(CellType.interval, 1, rule=basix.QuadratureType.gll)
#GLL quadrature of degree 1 corresponds to trapezoidal rule

x = [[], [], [], []]
for _ in range(3):
    x[0].append(np.zeros((0, 2)))
x[1].append(np.array([[1 - p[0], p[0]] for p in pts]))
x[1].append(np.array([[0, p[0]] for p in pts]))
x[1].append(np.array([[p[0], 0] for p in pts]))
x[2].append(np.zeros((0, 2)))

#interpolation matrix
M = [[], [], [], []]
for _ in range(3):
    M[0].append(np.zeros((0, 2, 0, 1))) #No DOFs on vertices
for normal in [[-1, -1], [-1, 0], [0, 1]]:
    mat = np.empty((2, 2, len(wts), 1)) #2 DOFs per edge, function 2D, no extra gradients
    mat[0, 0, :, 0] = normal[0] * np.array([1,0])
    mat[0, 1, :, 0] = normal[1] * np.array([1,0])
    mat[1, 0, :, 0] = normal[0] * np.array([0,1])
    mat[1, 1, :, 0] = normal[1] * np.array([0,1])
    M[1].append(mat)
M[2].append(np.zeros((0, 2, 0, 1)))

custom_bdm_element = basix.ufl.custom_element(
    basix.CellType.triangle,
    [2],
    wcoeffs,
    x,
    M,
    0,
    basix.MapType.contravariantPiola,
    basix.SobolevSpace.HDiv,
    False,
    0,
    1,
    dtype=default_real_type)

This then lets me define V = fem.functionspace(msh, custom_bdm_element) which I can use in forms etc.

1 Like