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.

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?