ElementQuadRT0

[gdmcbain]

Implementing the lowest order Raviart–Thomas quadrilateral was raised in #824 , which had begun by using to_meshtrion ex17.mesh. Reference there was made to Brezzi & Fortin, but this element was already described in §5 of the original paper (1977).

The basis functions on the canonical unit square are such that their normal components at the midpoints are the Kronecker delta. R & T didn't give them explicitly but by handcalculation, I think this means

phi = [
    [0, y-1],
    [x, 0],
    [0, y],
    [x-1, 0]]

with the divergence dphi being unity for each.

[gdmcbain]

Draft PR: #832

[gdmcbain]

A little script to document my working, showing where that phi came from.

import numpy as np
from numpy.polynomial import Polynomial

import skfem as fe

x = fe.Basis(fe.MeshQuad(), fe.ElementQuadRT0() * fe.ElementQuad0()).elem.elems[0].doflocs
unnormalized_normals = x - x.mean(axis=0)
normals = unnormalized_normals / np.linalg.norm(unnormalized_normals, axis=1)[:, None]
monomials = np.hstack([np.ones((x.shape[0], 1)), x, np.prod(x, axis=1)[:, None]])
A = normals @ np.tile(np.eye(2), (1, 2))  #  [[1, 0, x, 0], [0, 1, 0, y]], with (x, y) to be put in next
A[:, 2:] *= x
Ainv = np.linalg.inv(A)
phi_polys = np.array([[Polynomial(c) for c in coeffs] for coeffs in [[[1], [0], [0, 1], [0]], [[0], [1], [0], [0, 1]]]] @ Ainv).T
phi = [[list(p.coef) for p in pp] for pp in phi_polys]
print(phi.shape)

The output is

[[[0.0], [-1.0, 1.0]], [[0.0, 1.0], [0.0]], [[0.0], [0.0, 1.0]], [[-1.0, 1.0], [0.0]]]

which is interpreted as a list of pairs of polynomials (in x, and y, respectively), i.e.

• 0, −1 + y
• x, 0
• &c., as above.