Quasi-optimal hp estimates for projection-based interpolation on hexahedral elements

The 2006 paper of Rachowicz, Pardo and Demkowicz introduced the first fully automatic three-dimensional hp-adaptive finite element strategy, on hexahedral meshes, and identified as its principal unresolved theoretical issue the missing interpolation theory for the projection-based interpolation (PBI) operators actually implemented: the construction uses scaled $H^1$ edge/face products on 1-irregular meshes, and quasi-optimal hp error estimates were expected by analogy with the tetrahedral theory but not proved. Subsequent work completed the tetrahedral program (Demkowicz–Buffa 2005; Cao–Demkowicz 2005; Melenk–Rojik 2020), produced p-robust polynomial extension operators for the cube (Costabel–Dauge–Demkowicz 2008), and built alternative hp interpolation theories on hexahedra in the Jacobi-weighted framework (Guo–Zhang 2009; Guo–Yi 2013), but the quasi-optimality and commuting properties of the hexahedral PBI itself on variable-order, hanging-node meshes remain unproven. The surviving open problem is to prove or rigorously refute quasi-optimal hp estimates for the hexahedral exact-sequence projection-based interpolation operators.

Background

The hp-version of the finite element method (hp-FEM) approximates a boundary-value problem with piecewise polynomials whose local mesh size $h$ and local polynomial degree $p$ may both vary over the mesh; with a near-optimal distribution of $h$ and $p$ the error can decay exponentially in the number of degrees of freedom, even for solutions with corner and edge singularities. Exploiting this in three dimensions requires an algorithm that decides, automatically and element by element, whether to subdivide an element (h-refinement), raise its polynomial order (p-refinement), or both, possibly anisotropically. Rachowicz, Pardo and Demkowicz (2006) built the first fully automatic 3D hp strategy for elliptic problems on hexahedral meshes: it maintains a coarse and a uniformly hp-refined fine grid, solves on both, and drives every refinement decision by comparing local projection errors of the fine-grid solution onto candidate refinements of each coarse element. The engine of these comparisons is projection-based interpolation (PBI), the object of this problem.

PBI is an element-local, hierarchical interpolation operator. For an $H^1$ function on an element $K$ it first fixes the vertex values, then projects the traces on edges in fractional norms, then on faces, and finally computes the interior degrees of freedom as a minimum-energy correction — a constrained $L^2$-type projection with the boundary data already prescribed. Because the interpolant's trace on any edge or face depends only on the function's trace on that edge or face, interpolants assembled element by element are globally conforming, and the local projection errors serve as computable surrogates for the discretization error. In the 3D hp code the edge and face projections use scaled $H^1$ (energy-type) inner products, and the same construction extends through the exact sequence $H^1(K)\xrightarrow{\nabla}H(\mathrm{curl},K)\xrightarrow{\nabla\times}H(\mathrm{div},K)\xrightarrow{\nabla\cdot}L^2(K)$, where $H(\mathrm{curl})$ and $H(\mathrm{div})$ are the spaces of square-integrable vector fields with square-integrable curl or divergence, discretized by Nédélec edge elements (curl-conforming finite elements whose degrees of freedom live on mesh edges) and Raviart–Thomas-type face elements (divergence-conforming elements with face-based degrees of freedom); the corresponding operators $\Pi^{\mathrm{grad}}$, $\Pi^{\mathrm{curl}}$, $\Pi^{\mathrm{div}}$ are required to commute with the differential operators (the discrete de Rham diagram), which is what ties them to Maxwell applications. The mesh is 1-irregular: neighboring elements may differ by at most one refinement level, so edges and faces can carry hanging nodes — mesh points that must be handled by constrained approximation, coupling an element's projections to those of refined neighbors — and element orders may be anisotropic, i.e. a different polynomial degree per direction on the tensor-product hexahedral space.

The theoretical anchor is quasi-optimality: does the PBI error behave like the best-approximation error from the local polynomial space, up to a constant independent of $h$ and $p$ (at worst up to explicit logarithmic factors in $p$)? For tetrahedral elements this was established in a sequence of works — Demkowicz and Buffa (2005) proved quasi-optimal p-interpolation estimates for the full exact sequence (with the commuting diagram resting on an explicitly stated conjecture on p-robust polynomial extension operators), Cao and Demkowicz (2005) removed the $\varepsilon$-loss at the price of a $(\ln p)^{3/2}$ factor, and Melenk and Rojik (2020) proved commuting, logarithm-free optimal p-estimates on tetrahedra, building on the polynomial extension operators for tetrahedra. For hexahedra the picture is different: the 3D hp code implements the hexahedral PBI with scaled $H^1$ products, the authors expected quasi-optimal hp estimates analogous to the tetrahedral results, but no proof was available — the gap the source paper itself identifies as its principal unresolved theoretical issue.

Two structural obstacles make the hexahedral case harder than the tetrahedral one. First, the trace norms of order $1/2$ and $-1/2$ on the boundary of a cube are non-local, and decomposing them into face- and edge-local pieces is exactly what costs logarithmic factors even on tetrahedra. Second, tensor-product polynomial spaces with independent orders per direction do not interact with the vertex–edge–face hierarchy the way isotropic polynomial spaces on simplices do, and the 1-irregular hanging-node constraints tie each element's operator to its neighbors' refinement state. Moreover, PBI commutes with the element map only when the map is a pure scaling (plus rigid motion), so the usual route from reference-element p-estimates to hp-estimates on physical elements via a Bramble–Hilbert argument (the classical device that transfers approximation estimates from a reference element to physical elements) requires additional care on general hexahedra. Partial building blocks exist: Costabel, Dauge and Demkowicz (2008) constructed p-robust polynomial extension operators for $H^1$, $H(\mathrm{curl})$ and $H(\mathrm{div})$ spaces on the cube (for the full cube, not the variable-order hierarchical setting), and Guo and collaborators developed stable polynomial extensions and Chebyshev projection–interpolation on tetrahedral and hexahedral elements within the Jacobi-weighted Sobolev/Besov framework (function spaces weighted by Jacobi-polynomial weights, tailored to high-order polynomial approximation) — a different operator with best-approximation properties, but not the PBI operator used by the automatic algorithm.

Problem Statement

Prove, or rigorously disprove, quasi-optimality of the hexahedral projection-based interpolation operators — the construction implemented in the fully automatic 3D hp code (hierarchical vertex–edge–face–element projections with scaled $H^1$ edge/face products, variable and possibly anisotropic polynomial orders, applied on 1-irregular hexahedral meshes through constrained approximation) — in the sense established for tetrahedra.

Specifically: (a) for the $H^1$-conforming operator $\Pi^{\mathrm{grad}}$ on an element $K$, establish or refute hp-estimates of the form $\|u-\Pi^{\mathrm{grad}}u\|_{H^1(K)} \le C\,\Big(\frac{h_K}{p}\Big)^{r-1}\,(\ln p)^{\alpha}\,|u|_{H^r(K)}$ for $u\in H^r(K)$ in the natural regularity range of the tetrahedral theory (roughly $r>3/2$ for the $H^1$ operator), on affine images of the reference hexahedron, where $h_K$ is the diameter of $K$, with $C$ independent of $h_K$ and $p$ and $\alpha$ explicit — for anisotropic orders, $p$ is read as the minimal directional order of $K$, or the bound is stated in the corresponding per-direction form — the affirmative target for (a) being the quasi-optimality bound, i.e. boundedness of the PBI error by the best-approximation error from the local polynomial space up to a factor $C(\ln p)^{\alpha}$, of which the displayed rate inequality is a consequence; (b) the analogous estimates for the $H(\mathrm{curl})$- and $H(\mathrm{div})$-conforming operators $\Pi^{\mathrm{curl}}$ and $\Pi^{\mathrm{div}}$ with their natural regularity ranges; (c) the commutation relations $\nabla\,\Pi^{\mathrm{grad}} = \Pi^{\mathrm{curl}}\,\nabla,\quad \nabla\times\,\Pi^{\mathrm{curl}} = \Pi^{\mathrm{div}}\,(\nabla\times),\quad \nabla\cdot\,\Pi^{\mathrm{div}} = P\,(\nabla\cdot)$ with p-uniform stability, where $P$ is the $L^2$ projection; and (d) control of all constants with respect to the 1-irregular refinement configuration, in particular whether the constrained (hanging-node) projections remain stable uniformly in the refinement pattern.

A proof must establish all four sub-claims (a)–(d) together; a rigorous refutation of any one of them equally resolves the question — an explicit configuration (element or small patch, orders, family of functions indexed by $p$ or by refinement depth) in which, respectively, the ratio of the PBI error to the best-approximation error provably grows without bound, one of the commutation relations fails with p-uniform stability, or the constrained (hanging-node) projections are not uniformly stable over the 1-irregular refinement pattern would show that the implemented operator must be redesigned. Results restricted to tetrahedra, to two dimensions, or to the full cube with uniform polynomial degree do not answer the question: the target is the hexahedral, variable-order, hanging-node setting actually used by the automatic 3D hp algorithm.

References

1. W. Rachowicz, D. Pardo, and L. Demkowicz (2006). Fully automatic hp-adaptivity in three dimensions. Computer Methods in Applied Mechanics and Engineering 195, 4816–4842. DOI 10.1016/j.cma.2005.08.022, https://doi.org/10.1016/j.cma.2005.08.022 (preprint: ICES Report 04-22, https://www.oden.utexas.edu/media/reports/2004/0422.pdf)
2. L. Demkowicz and A. Buffa (2005). H1, H(curl) and H(div)-conforming projection-based interpolation in three dimensions. Computer Methods in Applied Mechanics and Engineering 194, 267–296. DOI 10.1016/j.cma.2004.07.007, https://doi.org/10.1016/j.cma.2004.07.007
3. W. Cao and L. Demkowicz (2005). Optimal error estimate of a projection based interpolation for the p-version approximation in three dimensions. Computers & Mathematics with Applications 50, 359–366. DOI 10.1016/j.camwa.2005.04.005, https://doi.org/10.1016/j.camwa.2005.04.005
4. M. Costabel, M. Dauge, and L. Demkowicz (2008). Polynomial extension operators for H1, H(curl) and H(div)-spaces on a cube. Mathematics of Computation 77, 1967–1999. DOI 10.1090/S0025-5718-08-02108-X, https://doi.org/10.1090/S0025-5718-08-02108-X
5. J. M. Melenk and C. Rojik (2020). On commuting p-version projection-based interpolation on tetrahedra. Mathematics of Computation 89, 45–87. DOI 10.1090/mcom/3454, https://doi.org/10.1090/mcom/3454 (preprint arXiv:1802.00197, https://arxiv.org/abs/1802.00197)
6. L. Demkowicz and J. Kurtz (2007). Projection-based interpolation and automatic hp-adaptivity for finite element discretizations of elliptic and Maxwell problems. ESAIM: Proceedings 21, 1–15. DOI 10.1051/proc:072101, https://doi.org/10.1051/proc:072101
7. J. Kurtz and L. Demkowicz (2007). A fully automatic hp-adaptivity for elliptic PDEs in three dimensions. Computer Methods in Applied Mechanics and Engineering 196, 3534–3545. DOI 10.1016/j.cma.2006.10.053, https://doi.org/10.1016/j.cma.2006.10.053
8. L. Demkowicz, J. Kurtz, D. Pardo, M. Paszyński, W. Rachowicz, and A. Zdunek (2008). Computing with hp-Adaptive Finite Elements, Vol. 2: Frontiers: Three Dimensional Elliptic and Maxwell Problems with Applications. Chapman & Hall/CRC. DOI 10.1201/9781420011692, https://doi.org/10.1201/9781420011692
9. B. Guo and L. Yi (2013). Local Chebyshev projection–interpolation operator and application to the h–p version of the finite element method in three dimensions. Mathematical Models and Methods in Applied Sciences 23(5), 777–801. DOI 10.1142/S0218202512500601, https://doi.org/10.1142/S0218202512500601
10. B. Guo and J. Zhang (2009). Stable and compatible polynomial extensions in three dimensions and applications to the p and h-p finite element method. SIAM Journal on Numerical Analysis 47, 1195–1225. DOI 10.1137/070688006, https://doi.org/10.1137/070688006
11. S. Henneking and L. Demkowicz (2023). hp3D User Manual, version 2 (September 2023). arXiv:2207.12211, https://arxiv.org/abs/2207.12211
12. S. Li and W. Zheng (2026). Trace-Preserving hp Interpolation and Polynomial Liftings on Conforming Hexahedral Meshes. arXiv:2606.02125, https://arxiv.org/abs/2606.02125
13. S. Li and W. Zheng (2026). Boundary-preserving hp interpolation and p-robust discrete harmonic extensions on tetrahedral meshes. arXiv:2606.02086, https://arxiv.org/abs/2606.02086
14. S. Li and W. Zheng (2026). hp-Optimal DG Approximation and Robust Schwarz Decompositions on One-Irregular Cubical Meshes. arXiv:2606.27728, https://arxiv.org/abs/2606.27728
15. D. Pardo, L. Demkowicz, C. Torres-Verdín, and M. Paszyński (2007). A self-adaptive goal-oriented hp-finite element method with electromagnetic applications. Part II: Electrodynamics. Computer Methods in Applied Mechanics and Engineering 196, 3585–3597. DOI 10.1016/j.cma.2006.10.016, https://doi.org/10.1016/j.cma.2006.10.016
16. M. Paszyński and L. Demkowicz (2006). Parallel, fully automatic hp-adaptive 3D finite element package. Engineering with Computers 22, 255–276. DOI 10.1007/s00366-006-0036-8, https://doi.org/10.1007/s00366-006-0036-8
17. I. Gomez-Revuelto, L. E. García-Castillo, S. Llorente-Romano, and D. Pardo (2012). A three-dimensional self-adaptive hp finite element method for the characterization of waveguide discontinuities. Computer Methods in Applied Mechanics and Engineering 249–252, 62–74. DOI 10.1016/j.cma.2012.05.013, https://doi.org/10.1016/j.cma.2012.05.013
18. N. Zander, T. Bog, M. Elhaddad, F. Frischmann, S. Kollmannsberger, and E. Rank (2016). The multi-level hp-method for three-dimensional problems: Dynamically changing high-order mesh refinement with arbitrary hanging nodes. Computer Methods in Applied Mechanics and Engineering 310, 252–277. DOI 10.1016/j.cma.2016.07.007, https://doi.org/10.1016/j.cma.2016.07.007
19. V. Darrigrand, D. Pardo, T. Chaumont-Frelet, I. Gómez-Revuelto, and L. E. García-Castillo (2020). A painless automatic hp-adaptive strategy for elliptic problems. Finite Elements in Analysis and Design 178, 103424. DOI 10.1016/j.finel.2020.103424, https://doi.org/10.1016/j.finel.2020.103424
20. F. V. Caro, V. Darrigrand, J. Alvarez-Aramberri, E. Alberdi, and D. Pardo (2022). A painless multi-level automatic goal-oriented hp-adaptive coarsening strategy for elliptic and non-elliptic problems. Computer Methods in Applied Mechanics and Engineering 401, 115641. DOI 10.1016/j.cma.2022.115641, https://doi.org/10.1016/j.cma.2022.115641
21. A. Chakraborty, S. Henneking, and L. Demkowicz (2024). An anisotropic hp-adaptation framework for ultraweak discontinuous Petrov–Galerkin formulations. Computers & Mathematics with Applications 167, 315–327. DOI 10.1016/j.camwa.2024.05.025, https://doi.org/10.1016/j.camwa.2024.05.025
22. M. Bürg (2013). Convergence of an automatic hp-adaptive finite element strategy for Maxwell's equations. Applied Numerical Mathematics 72, 188–204. DOI 10.1016/j.apnum.2013.04.008, https://doi.org/10.1016/j.apnum.2013.04.008
23. P. Daniel and M. Vohralík (2023). Guaranteed contraction of adaptive inexact hp-refinement strategies with realistic stopping criteria. ESAIM: Mathematical Modelling and Numerical Analysis 57, 329–366. DOI 10.1051/m2an/2022082, https://doi.org/10.1051/m2an/2022082