Arnold-Boffi-Falk spaces using Basix

General
1

Hello everyone,

I want to use the Arnold-Boffi-Falk H(div)-conforming elements (see DefElement: Arnold–Boffi–Falk) in my FEniCSx scripts. As far as I know, they are not implemented by default, but I could define them using Basix.

I have no experience using Basix, but I found this helpful demo about the lowest-order RT element on triangles.

https://docs.fenicsproject.org/basix/v0.5.1/python/demo/demo_custom_element.py.html#degree-1-ravairt-thomas-element

Since the Arnold-Boffi-Falk elements are somewhat similar to the RT elements on quadrilaterals (see DefElement: Q H(div)), I figured that a good first step would be adapting the demo to generate the lowest-order RT element on quadrilaterals. Here are my efforts.

import basix
import numpy as np
from basix import CellType, MapType, PolynomialType, LatticeType, SobolevSpace

# Coefficient matrix
# ------------------

wcoeffs = np.zeros((4, 8))

pts, wts = basix.make_quadrature(CellType.quadrilateral, 2)
poly = basix.tabulate_polynomials(PolynomialType.legendre, CellType.quadrilateral, 1, pts)
x = pts[:, 0]
y = pts[:, 1]
for i in range(4):
    wcoeffs[0, 4+i] = sum((1-y)*poly[i, :]*wts)
    wcoeffs[1, i] = sum((x-1)*poly[i, :]*wts)
    wcoeffs[2, i] = sum(-x*poly[i, :]*wts)
    wcoeffs[3, 4+i] = sum(y*poly[i, :]*wts)

# Interpolation
# -------------

pts, wts = basix.make_quadrature(CellType.interval, 1)

# The points associated with each edge are calculated by mapping the quadrature points to each edge.

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

M = [[], [], [], []]
for _ in range(4):
    M[0].append(np.zeros((0, 2, 0, 1)))
for normal in [[0, 1], [-1, 0], [-1, 0], [0, 1]]:
    mat = np.empty((1, 2, len(wts), 1))
    mat[0, 0, :, 0] = normal[0]*wts
    mat[0, 1, :, 0] = normal[1]*wts
    M[1].append(mat)
M[2].append(np.zeros((0, 2, 0, 1)))

# Creating the element
# --------------------

custom_element = basix.create_custom_element(
    CellType.quadrilateral, [2], wcoeffs, x, M, 0, MapType.contravariantPiola, SobolevSpace.HDiv,
    False, 0, 1)

I don’t get any errors when I run this code, but the resulting custom element leads to different results when compared to the RT element implemented in UFL (RTCF of degree 1).

Can someone help me figure out where the problem is?

2

I fixed a bug in how DOLFINx treats custom elements yesterday (Fix DOF transformations for custom elements (#2588) · FEniCS/dolfinx@d9620f2 · GitHub). This bug would affect the element you made, so you’ll get better results if you rebuild DOLFINx with this patch.

FEniCSx can’t yet exactly support these elements, as they use integrals of the divergence in their functionals, and elements with derivatives in their definition currently don’t work with FFCx and DOLFINx. But you might be able to cook up some alternative definitions of the functionals that work. (I’d have to read the ABF paper again to work out how important the derivatives being involved is.)

3

Hi Matthew, thanks for your reply!

I will try my code again using the new patch.

About the ABF elements, I don’t know a way to circumvent the div in their definitions. But I wonder if it would be possible to define them as a mixed element RT + some bubble functions. Please let me know if you get any updates on the subject.