The high-order finite element Duffy de Rham complex and low-order-refined preconditioning

Abstract: In this work, we construct high-order finite element spaces for the $L^2$ de Rham complex on triangular meshes amenable to low-order-refined preconditioning. The spaces are constructed using the Duffy transformation, by pulling back appropriately chosen polynomial spaces defined on the unit square; in addition to piecewise polynomials, these spaces also contain certain rational functions, and they reduce to the standard Lagrange, Nédélec, and discontinuous finite elements in the lowest-order case. We establish spectral equivalence, independent of the polynomial degree, of the stiffness matrices defined on these spaces with the lowest-order stiffness matrices defined on refined meshes, constructed using a Gauss-Lobatto triangular lattice. Spectral equivalence of the operators is a consequence of norm equivalences in Jacobi-weighted $L^2$ norms, which are established by proving stability of the Jacobi-Gauss-Lobatto interpolation operator in shifted norms. The low-order-refined preconditioners can also be used to precondition the standard piecewise polynomial finite element spaces using a fictitious space approach. The low-order-refined system can in turn be preconditioned effectively using algebraic multigrid methods. The analytical estimates are confirmed by numerical results on a variety of high-order problems, including on mixed meshes and surface meshes.

• The paper introduces high-order Duffy-mapped finite element spaces for triangulated L2 de Rham complexes, enabling degree-robust preconditioning.
• It rigorously demonstrates uniform spectral equivalence between high-order and low-order-refined operators using detailed Jacobi-weighted norm analysis.
• Numerical experiments confirm the efficiency and scalability of matrix-free implementations and robust preconditioners across mesh and polynomial refinements.

High-Order Finite Element Duffy de Rham Complex and Low-Order-Refined Preconditioning

Introduction and Motivation

The paper "The high-order finite element Duffy de Rham complex and low-order-refined preconditioning" (2604.00148) introduces a novel class of high-order finite element spaces for the $L^2$ de Rham complex on triangular meshes, leveraging the Duffy transformation to enable efficient and degree-robust preconditioning. The aim is to generalize the low-order-refined (LOR) preconditioning methodology, well established in tensor-product (quadrilateral and hexahedral) contexts, to simplicial (triangular) meshes while preserving operator spectral equivalence independently of the polynomial degree.

The current state of preconditioning in high-order finite elements relies heavily on the equivalence of high-order discretizations with low-order operators on refined meshes, enabling efficient matrix-free iterative solvers. However, non-tensorial geometries like triangles present significant challenges: standard transfer operators become non-trivial, and ensuring that discrete differential operators retain suitable properties requires careful construction of both spaces and degrees of freedom (DOFs). This work addresses this gap through the design of new finite element complexes and the analysis of their spectral equivalence properties.

Construction of High-Order Duffy-Mapped Spaces

The central contribution is the formulation and analysis of finite element spaces for the $L^2$ de Rham complex on triangles via pullbacks under the Duffy transformation. The element mapping from the reference square to the reference triangle introduces (a) rational extensions to classical polynomial spaces and (b) specific Jacobi-weighted structures.

Letting $\varphi: [0,1]^2 \to \Delta$ denote the Duffy map, local spaces on the triangle are defined as pullbacks of structured polynomial or tensor-product polynomial spaces on the square, mapped with appropriate Piola-type transformations for each function space in the complex. The resulting spaces, which include standard Lagrange, Nédélec, and discontinuous elements in the lowest-order case, are generalizations containing both standard polynomials and Duffy-induced rational functions.

Degrees of freedom for these spaces are provided by mapping tensor-product Gaussian or Gauss–Lobatto point lattices from the square to the triangle, resulting in a "collapsed lattice" on the triangle. These DOFs induce element-level global finite element spaces with precisely defined continuity and compatibility properties, enabling the construction of global spaces that respect the $L^2$ finite element de Rham complex on the mesh.

Crucially, the discrete differential operators (grad, curl) have matrix representations coinciding with the incidence matrices of the triangular lattice, yielding structurally sparse and algebraically simple operators well-suited for multigrid and domain decomposition methods.

Low-Order-Refined Preconditioning Theory

Low-order-refined preconditioning shows that, for high-order discretizations on a mesh, there exists a suitably refined mesh on which corresponding low-order (piecewise-linear or lowest-order) discretizations are spectrally equivalent to the original high-order operators, independent of the polynomial order $N$. In this paper, the refined mesh is obtained by superimposing the Duffy-transformed lattice on each triangle, decomposing each triangle into subcells corresponding to the mapped lattice structure.

Central results establish norm equivalence estimates across polynomial degree, even in Jacobi-weighted $L^2$ and $H^1$ seminorms, and uniform spectral equivalence of the high-order and low-order-refined stiffness/mass matrices. The proof leverages intricate 1D interpolation stability results for Gauss–Lobatto interpolation in shifted Jacobi-weighted Sobolev norms. Ultimately, for stiffness matrices $A$ on the high-order space and $A_0$ on the low-order-refined space, the condition number $\kappa(A_0^{-1} A)$ is shown to be bounded independent of $L^2$0.

The preconditioners acquired by replacing the high-order system with the low-order-refined system are highly effective because:

• The LOR matrices are extremely sparse (comparable sparsity to standard lowest-order finite element systems).
• The identification of DOFs between high- and low-order-refined spaces leads to transfer operators that are either trivial or easily computable, further simplifying the preconditioning algorithm.

The mass matrices, however, do not generally enjoy uniform spectral equivalence under this framework; their conditioning with respect to the LOR system grows linearly with $L^2$1 for standard Gauss–Lobatto point sets, an explicit consequence of the weighted norm effects induced by the Duffy transformation.

Numerical Results

Numerical experiments confirm the theoretical findings:

• Uniform spectral equivalence is observed in numerically computed condition numbers for all high-order Laplacian, curl–curl, and $L^2$2 mass operators relative to the LOR systems over a wide range of polynomial degrees ($L^2$3 up to 128).
• The sparsity of the LOR preconditioners matches well with expectations: the number of nonzeros per row in their matrices approaches 9 for the $L^2$4 case, regardless of $L^2$5, even as the high-order matrices' per-row density grows as $L^2$6.
• Realistic unstructured meshes and mixed meshes (with triangles and quadrilaterals) validate that LOR preconditioning is robust to mesh geometry and anisotropy, and interfaces between element types.
• Matrix-free implementations enabled by the Duffy mapping and sum factorization demonstrate clear asymptotic efficiency gains, especially at high polynomial orders, compared to assembled-matrix-based methods.

The preconditioning methodology proved robust not only under order ($L^2$7-) refinement but also under mesh ($L^2$8-) refinement—iteration counts for preconditioned iterative solvers remain essentially constant under both.

Furthermore, fictitious space preconditioning is demonstrated: the high-order Duffy-mapped spaces serve as auxiliary (fictitious) spaces to precondition the standard $L^2$9 Lagrange space via elliptic local projections, yielding degree-robust preconditioners for classical polynomial finite element spaces.

Theoretical and Practical Implications

This work has both theoretical and practical implications:

• Theoretically, it completes an essential extension of LOR preconditioning to simplicial meshes, previously only available for tensor-product meshes.
• The construction supports mixed-mesh compatibility, allowing seamless coupling of triangles and quadrilaterals, and paves the way for robust LOR preconditioning in hybridized and high-order discontinuous Galerkin methods on unstructured grids.
• The detailed Jacobi-weighted norm analysis extends interpolation theory, establishing new stability and norm-equivalence results for high-order interpolants, which may find broader application in numerical analysis of spectral methods.
• The methodology provides a foundation for highly efficient, matrix-free, high-order operator evaluation and solvers, which is critical for extreme-scale scientific computing.

Limitations remain in the mass matrix preconditioning aspect for certain weightings and in extending these constructions to 3D (especially tetrahedra and prisms), which is noted as a future direction.

Conclusion

The construction of high-order finite element spaces for the $\varphi: [0,1]^2 \to \Delta$0 de Rham complex via the Duffy transformation, together with rigorous analysis and practical demonstration of low-order-refined preconditioning, constitutes a substantial advance in robust, scalable high-order finite element technology for triangular meshes. The spectral equivalence results, confirmed in challenging computational regimes, solidify the theoretical underpinnings of LOR preconditioning outside the tensor-product paradigm and open new avenues for high-performance implementation of high-order methods on general unstructured meshes.

Future research is expected to target the extension of this framework to three-dimensional elements, further optimization of mass matrix preconditioning (including optimal node selection and weighting), and development of efficient matrix-free local solvers for fictitious space transfer operators. The approaches here are likely to inform the development of next-generation finite element software and preconditioning strategies for high-fidelity PDE simulation at scale.