We provide an implementation of the \(\mathrm {P}_p\)  and \(\mathrm {DP}_{p-1}\)  elements with the
FDM basis functions on the interval in FIAT [1]}. The
extension to quadrilaterals and hexahedra is achieved by taking tensor products
of the one-dimensional elements with FInAT [2]}. Code for the
sum-factorized evaluation of the residual is automatically generated by
Firedrake [3]}, [4]}, implementing
a Gauß–Lobatto quadrature rule with \(3(p+1)/2\)  points along each direction.
The sparse preconditioner discretizing the auxiliary operator is implemented as
a PETSc [5]} preconditioner as \(\texttt {firedrake.FDMPC}\) .
The preconditioners are additive on each level and multiplicative between levels.
The Cholesky factorization of the patch matrices is computed using
CHOLMOD [6]} and the ICC factorization is done with PETSc's own
implementation.
Most of our computations were performed on a single node
of the ARCHER2 system, with two 64-core AMD EPYC 7742 CPUs (2.25 GHz) and 512 GiB of RAM.
