# A quasi-interpolation operator yielding fully computable error bounds

Théophile Chaumont-Frelet, Martin Vohralík

► To cite this version:

Théophile Chaumont-Frelet, Martin Vohralík. A quasi-interpolation operator yielding fully computable error bounds. ESAIM: Mathematical Modelling and Numerical Analysis, In press. hal-05164027

HAL Id: hal-05164027

https://inria.hal.science/hal-05164027v1

Submitted on 16 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A QUASI-INTERPOLATION OPERATOR YIELDING FULLY
COMPUTABLE ERROR BOUNDS

T. CHAUMONT-FRELET\* AND M. VOHRALÍK\textasciitilde, \textasciitilde

ABSTRACT. We design a quasi-interpolation operator from the Sobolev space $H_0^1(\Omega)$ to its finite-
dimensional finite element subspace formed by piecewise polynomials on a simplicial mesh with
a computable approximation constant. The operator 1) is defined on the entire $H_0^1(\Omega)$, no
additional regularity is needed; 2) allows for an arbitrary polynomial degree; 3) works in any
space dimension; 4) is defined locally, in vertex patches of mesh elements; 5) yields optimal
estimates for both the $H^1$ seminorm and the $L^2$ norm error; 6) gives a computable constant
for both the $H^1$ seminorm and the $L^2$ norm error; 7) leads to the equivalence of global-best
and local-best errors; 8) possesses the projection property. Its construction follows the so-called
potential reconstruction from a posteriori error analysis. Numerical experiments illustrate that
our quasi-interpolation operator systematically gives the correct convergence rates in both the
$H^1$ seminorm and the $L^2$ norm and its certified overestimation factor is rather sharp and stable
in all tested situations.

KEYWORDS. finite element method; interpolation operator; stable projection; error estimate; guaranteed bound; minimal regularity

# 1. INTRODUCTION

Due to their ability to operate on general geometries, finite element methods based on un-
structured meshes have become very popular to discretize boundary value problems over the past
decades [9, 14]. The approximation properties of the finite element spaces are one of the key
factors that eventually govern the accuracy of the resulting numerical schemes. Such approxima-
tion properties are often obtained thanks to (quasi-)interpolation operators, i.e., mappings that
associate to a given function $u$ from the infinite-dimensional Sobolev space $H_0^1(\Omega)$ an explicitly
constructed approximation $\mathcal{I}_h u$ in the finite-dimensional finite element space.

Given a polytopal domain $\Omega \subset \mathbb{R}^d$, $d \ge 1$, triangulated by a simplicial mesh $\mathcal{T}_h$, this work focuses on conforming finite elements of polynomial degree $p \ge 1$, $\mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$, forming a finite-dimensional subspace of the Sobolev space $H_0^1(\Omega)$. In this context, a natural way to associate to a function $u \in H_0^1(\Omega)$ an interpolant $\mathcal{I}_h u \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$ is by requiring that $(\mathcal{I}_h u)(\mathbf{x}_\ell) = u(\mathbf{x}_\ell)$ for suitably chosen interpolation nodes $\{\mathbf{x}_\ell\}_\ell$. This process is known as Lagrange interpolation, and presents the advantage that it is fully local: For a mesh element $K \in \mathcal{T}_h$, the interpolation error $(u - \mathcal{I}_h u)|_K$ will only depend on $u|_K$. Besides, fully-explicit error bounds expressed in terms of Sobolev semi-norms of $u$ are available in the literature [1, 6, 18, 20, 21, 22, 24, 25].

The main drawback of Lagrange interpolation, however, is that it requires the nodal values
of $u$ to be well-defined, which is not the case (expect in one space dimension, $d = 1$) if $u$
barely sits in $H_0^1(\Omega)$. When considering, e.g., the Poisson problem in a convex domain with
an $L^2(\Omega)$ right-hand side, this is not an issue, since then the solution $u$ belongs to $H^2(\Omega)$,
enabling the use of Lagrange interpolation. However, the approach fails when considering more
complex geometries and/or coefficients. Additionally, some applications, such as a posteriori
error estimation [4, 27], preconditioning [33], or localized orthogonal decomposition methods [26],
crucially hinge on interpolation operators defined for generic $H_0^1(\Omega)$ target functions (i.e. with
the minimal variational regularity of the PDE under consideration).

<footnote>\*Inria, Univ. Lille, CNRS, UMR 8524 – Laboratoire Paul Painlevé, 40 Av. Halley, 59650 Villeneuve-d'Ascq, France</footnote>

<footnote>†Inria, 48 rue Barault, 75647 Paris, France</footnote>

<footnote>‡CERMICS, Ecole nationale des ponts et chaussées, IP Paris, 77455 Marne-la-Vallée, France</footnote>
There is therefore a need to define interpolation operators for non-smooth functions. The resulting operators are often called “quasi-interpolation” operators, since they do not, strictly-speaking, interpolate the target at any point. The seminal contribution in this direction is the work of Clément [10], where an operator defined over the whole $L^1(\Omega)$ space is introduced with optimal approximation properties. The Clément operator does not preserve boundary conditions, and it is not a projection, i.e., $\mathcal{I}_h u_h \neq u_h$ for $u_h$ already in the finite element space. These two properties have been incorporated later on in the key contribution of Scott and Zhang [34]. For these two operators, the degrees of freedom (dofs) of the interpolant $\mathcal{I}_h u$ are fixed by suitably weighted mean values of the target function $u$. Another family of quasi-interpolation operators is due to Oswald [29]. There, a discontinuous approximation is first built by projecting element-wise $u|_K$ onto a local polynomial space $\mathcal{P}_p(K)$ of the mesh cell $K \in \mathcal{T}_h$. The dofs of this approximation that are shared by several elements are then suitably averaged to produce a conforming interpolant $\mathcal{I}_h u \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$. Since then, these constructions have been improved in several ways, and the design of quasi-interpolation operators is still an active field of research [13, 17, 16, 19, 23, 32, 40].

Because the nodal evaluation of a generic target function $u \in H_0^1(\Omega)$ is not possible, quasi-interpolation operators are not fully local. Indeed, the value of dofs needs to be obtained through some suitable averaging, meaning that the definition $(\mathcal{I}_h u)|_K$ will not only depend on $u|_K$, but rather on $u|_{\omega_K}$, where $\omega_K$ is the domain formed by other mesh elements surrounding $K$. Although this is not a huge drawback in practice, this fact makes the error analysis more complicated than for standard interpolation operators, since it is not possible to use arguments involving a single “reference element”. As a result, a shape-regularity requirement involving $\omega_K$ rather than each element $K$ individually often appears. In fact, the only fully-explicit approximation results we are aware of in this direction are [5, 38, 36, 37] but the scope of these results is limited, as they are specifically established with applications to a posteriori error estimation in mind.

Another important topic related to our work and connected to quasi-interpolation operators is the comparison between local-best and global-best approximations. Namely, we ask whether there exists a constant $C > 0$ such that

$$ (1.1) \qquad \min_{u_h \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)} \|\nabla(u - u_h)\|_{\Omega} \le C \min_{v_h \in \mathcal{P}_p(\mathcal{T}_h)} \|\nabla_h(u - v_h)\|_{\Omega}. $$

In other words, we ask whether, up to a generic constant, the conforming finite element approximation is as good as the element-wise broken polynomial approximation of the target function. It is in fact clear that such a constant exists, because the left-hand side vanishes whenever the right-hand side does. However, the dependence of the constant $C$ on key discretization parameters is a subtle issue. Following [2, 7, 35], recent results in this direction are given in [17, 23, 40]. Notice that once an inequality such as (1.1) is established, approximability estimates for the finite element space $\mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$ follow, since the right-hand side may be easily estimated in terms of (broken) Sobolev norms of $u$ using element-wise Poincaré inequalities (see Section 3.2 below).

In this work, we propose a new quasi-interpolation operator $\mathcal{J}_h^p : H_0^1(\Omega) \to \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$. This operator can be computed locally by solving patch-wise finite element projection problems, works for any polynomial degree $p \ge 1$, preserves boundary conditions, is a projection, and, moreover, satisfies

$$ (1.2) \qquad \|\nabla(u - \mathcal{J}_h^p(u))\|_{\Omega} \le C \min_{v_h \in \mathcal{P}_p(\mathcal{T}_h)} \|\nabla_h(u - v_h)\|_{\Omega} \quad \forall u \in H_0^1(\Omega) $$

for a constant $C$ only depending on the shape-regularity parameter of the mesh $\mathcal{T}_h$, the dimension $d$, and the polynomial degree $p$. In particular, our interpolation operator has the approximation power of discontinuous piecewise polynomials and leads to the local-global equivalence (1.1). We also establish localized versions of (1.2), as well as error estimates in the $L^2(\Omega)$ norm (see Section 4 below). Our construction is inspired by the so-called potential reconstructions used in the context of a posteriori error estimation for nonconforming and mixed finite element methods [15], as well as recent works of the authors on commuting quasi-interpolation operators [8, 12].
The main novelty of our work is that the constant C appearing in (1.2) is fully computable. Specifically, it is calculated through an algorithm that amounts to solving small, uncoupled, matrix eigenvalue problems that stem from patch-wise finite element spaces. We also establish localized versions of (1.2) and $L^2(\Omega)$ error estimates with fully-computable constants. To the best of our knowledge, this work is the first to provide computable constants under minimal regularity, for arbitrary polynomial degree, and yielding the comparison of global-best and local-best errors. We emphasize that the proposed algorithm computes the generic constant C in (1.2) valid for all $u \in H_0^1(\Omega)$ (without needing to know u).

The remainder of this work is organized as follows. Section 2 gives the setting and recalls key preliminary results. In Section 3, we present elementwise Lagrange interpolation, the local-best approximation, and motivate our approach. The construction of our quasi-interpolation operator and statement of its key properties including (1.2) forms the content of Section 4. Finally, Section 5 collects the proofs and Section 6 reports on the actual behavior of our quasi-interpolation operator in several numerical experiments.

## 2. MESHES, SPACES, PATCHES, AND POINCARÉ INEQUALITY

### 2.1. Computational mesh.
We consider an open, bounded, connected, Lipschitz polyhedral domain $\Omega \subset \mathbb{R}^d$, $d \ge 1$, partitioned into a mesh $\mathcal{T}_h$ that consists of (open) simplicial elements $K$. We assume that the mesh is matching in the sense that the intersection $\bar{K}_+ \cap \bar{K}_-$ of two elements $K_\pm \in \mathcal{T}_h$ is either empty or a single common $d'$-dimensional subsimplex of $\bar{K}_+$ and $\bar{K}_-$, $0 \le d' \le d-1$ (a single vertex, edge, or face of $K_-$ and $K_+$ if $d=3$). This assumption is standard (see, e.g., [9, Section 2.2] or [14, Definition 6.11]) We denote by $\mathcal{V}_h$ the set of vertices of $\mathcal{T}_h$.

For a simplex $K \subset \mathbb{R}^d$, $h_K$ denotes the diameter of $K$ and $\rho_K$ the diameter of the largest ball contained in $\bar{K}$. The shape-regularity parameter of $K$ is then defined by

$$ \kappa_K := \frac{h_K}{\rho_K}. $$

If $\mathcal{T} \subset \mathcal{T}_h$ is a submesh, the shape-regularity parameter of $\mathcal{T}$ is $\kappa_{\mathcal{T}} := \max_{K \in \mathcal{T}} \kappa_K$. Notice that $\kappa_K \ge 1$ and that it is possible to design strongly graded meshes $\mathcal{T}_h$ while maintaining $\kappa_{\mathcal{T}_h}$ bounded, as long as the elements remain isotropic.

If $K \subset \mathbb{R}^d$ is a simplex, we employ the notation $\mathcal{F}_K$ for its $(d-1)$-dimensional faces and $\mathcal{V}_K$ for its vertices. Then, if $a \in \mathcal{V}_K$ is a vertex of $K$, the notation $\tau_K^a$ stands for the distance between $a$ and the (hyper)plane generated by the face of $K$ opposite to $a$. Notice that $\rho_K \le \tau_K^a \le h_K$, so that in particular $h_K/\tau_K^a \le \kappa_K$.

### 2.2. Sobolev spaces.
Throughout this manuscript, if $\omega \subset \Omega$ is an open, bounded, connected, Lipschitz subset of $\Omega$, then $L^2(\omega)$ is the Lebesgue space of scalar-valued square-integrable functions and $L^2(\omega) := [L^2(\omega)]^d$ is the space of vector-valued square-integrable functions. We denote by $(\cdot, \cdot)_\omega$ and $\|\cdot\|_\omega$ the inner products and norms of both spaces. We employ the notation $L^\infty(\omega)$ and $L^\infty(\omega)$ for essentially bounded scalar- and vector-valued functions, and denote by $\|\cdot\|_{L^\infty(\omega)}$ and $\|\cdot\|_{L^\infty(\omega)}$ their usual norms.

The Sobolev spaces

$$ H^1(\omega) := \{v \in L^2(\omega) \mid \nabla v \in L^2(\omega)\}, \quad H(\operatorname{div}, \omega) := \{w \in L^2(\omega) \mid \nabla \cdot w \in L^2(\omega)\}, $$

where $\nabla$ and $\nabla \cdot$ respectively denote the gradient and divergence operators defined in the sense of distributions, will be useful. We will also employ the following notation:

$$ H_0^1(\omega) := \{v \in H^1(\omega) \mid v = 0 \text{ on } \partial\omega\}, \quad H(\operatorname{div}^0, \omega) := \{v \in H(\operatorname{div}, \omega) \mid \nabla \cdot v = 0\}, $$

where the boundary value is understood in the trace sense.

For scalar-valued functions, we will also use higher-order Sobolev spaces. Namely for an integer $s \in \mathbb{N}$, $H^s(\omega)$ stands for the space of functions $v \in L^2(\omega)$ such that

$$ \partial_\alpha v \in L^2(\omega) $$
for all $\alpha \in \mathbb{N}^d$ with $|\alpha|_1 \le s$, whereby $|\cdot|_1$ denotes the $l_1$ norm on $\mathbb{N}^d$. We equip $H^s(\Omega)$ with the seminorm

$$ (2.1) \qquad |v|_{H^s(\omega)}^2 := \sum_{\substack{\alpha \in \mathbb{N}^d \\ |\alpha|_1 = s}} \|\partial_\alpha v\|_{\omega}^2 \quad \forall v \in H^s(\omega). $$

The broken Sobolev spaces

$$ H^s(\mathcal{T}_h) = \{v \in L^2(\Omega) \mid v|_{K} \in H^s(K) \quad \forall K \in \mathcal{T}_h \} $$

will also be useful. The broken gradient is in particular defined element-wise as

$$ (\nabla_h v)|_K := \nabla(v|_K) \quad \forall K \in \mathcal{T}_h $$

for all $v \in H^1(\mathcal{T}_h)$. It maps $H^1(\mathcal{T}_h)$ into $L^2(\Omega)$.

**2.3. Finite element spaces.** If $K \subset \mathbb{R}^d$ is a simplex and $q \ge 0$, we denote by $P_q(K)$ the set of polynomial functions of total degree at most $q$ on $K$ and $P_q(K) := [P_q(K)]^d$. The set of Raviart-Thomas polynomials on $K$ is then $RT_q(K) := xP_q(K) + P_q(K)$, see [28, 31].

If $\mathcal{T}$ is a set of simplices with the corresponding domain $\omega \subset \Omega$, then $P_q(\mathcal{T})$ collects the functions $v : \omega \to \mathbb{R}$ such that $v|_K \in P_q(K)$ for all $K \in \mathcal{T}$. The elementwise (broken) space $RT_q(\mathcal{T})$ is defined analogously.

Below, we fix a polynomial degree $p \ge 1$.

**2.4. Hat functions.** Throughout this work, the hat functions $\{\psi_a\}_{a \in \mathcal{V}_h}$ will play an important role. For each mesh vertex $a \in \mathcal{V}_h$, $\psi_a$ is the only element of $P_1(\mathcal{T}_h) \cap H^1(\Omega)$ such that $\psi_a(b) = \delta_{a,b}$ for all $b \in \mathcal{V}_h$, where $\delta$ is the Kronecker symbol. If $K \in \mathcal{T}_h$ and $a \in \mathcal{V}_K$, then

$$ (2.2) \qquad \|\psi_a\|_{L^\infty(K)} = 1 \qquad \|\nabla \psi_a\|_{L^\infty(K)} = \frac{1}{\tau_K^a}. $$

We also denote by $\omega_a$ the open domain corresponding to the support of $\psi_a$. Crucially, the hat functions form a partition of unity as

$$ (2.3) \qquad \sum_{a \in \mathcal{V}_h} \psi_a = 1. $$

**2.5. Vertex patches.** If $a \in \mathcal{V}_h$ we define the patch

$$ \mathcal{T}_a = \{K \in \mathcal{T}_h \mid a \in \mathcal{V}_K\}, $$

of elements having $a$ as a vertex. The elements $K \in \mathcal{T}_a$ correspond to the support $\omega_a$ of the hat function $\psi_a$. For a vertex patch $\mathcal{T}_a$, we define

$$ (2.4) \qquad h_a := \max_{K \in \mathcal{T}_h | a \in \mathcal{V}_K} h_K. $$

**2.6. Poincaré inequality.** Since simplices are convex, for all $K \in \mathcal{T}_h$ and for all $v \in H^1(K)$ such that $(v, 1)_K = 0$, we have

$$ (2.5) \qquad \|v\|_K \le \frac{h_K}{\pi} \|\nabla v\|_K, $$

see, e.g., [3, 30].

### 3. ELEMENTWISE LAGRANGE INTERPOLATION, LOCAL-BEST APPROXIMATION, AND MAIN IDEAS

In this section, we present the elementwise Lagrange interpolation operator and discuss local-best approximation by $H^1$ seminorm elementwise orthogonal projection. We then motivate the construction of our quasi-interpolation operator.
3.1. Elementwise Lagrange interpolation operator. We start by recalling the standard Lagrange interpolation operator on each mesh element $K \in \mathcal{T}_h$. We will only apply it to functions that are polynomials of degree $p+1$ on $K$, so that we have $\mathcal{I}_K^p : \mathcal{P}_{p+1}(K) \to \mathcal{P}_p(K)$,

$$ (3.1a) \qquad (\mathcal{I}_K^p v_h)(\mathbf{x}_\ell) = v_h(\mathbf{x}_\ell) $$

for all Lagrange interpolation nodes $\{\mathbf{x}_\ell\}_\ell$ on $K$ see, e.g., [14, Section 7.4]. Then, the elementwise Lagrange interpolation operator $\mathcal{I}_h^p : \mathcal{P}_{p+1}(\mathcal{T}_h) \to \mathcal{P}_p(\mathcal{T}_h)$ is defined by

$$ (3.1b) \qquad (\mathcal{I}_h^p v)|_K := \mathcal{I}_K^p(v|_K) \quad \forall K \in \mathcal{T}_h, $$

for $v \in \mathcal{P}_{p+1}(\mathcal{T}_h)$. Let us point out that if $v \in \mathcal{P}_{p+1}(\mathcal{T}_h) \cap H_0^1(\Omega)$, then $\mathcal{I}_h^p v \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$ and that $\mathcal{I}_h^p$ is a projection.

**Remark 3.1 (Other interpolation operators).** We focus on the Lagrange interpolation operator $\mathcal{I}_K^p$ on each mesh element $K$ for the sake of simplicity. In fact, any affine-equivalent interpolation operator $\mathcal{I}_K^p : \mathcal{P}_{p+1}(K) \to \mathcal{P}_p(K)$ with the property that $\mathcal{I}_h^p$ maps $\mathcal{P}_{p+1}(\mathcal{T}_h) \cap H_0^1(\Omega)$ into $\mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$ can be used.

3.2. Local-best approximation. Another key ingredient of our quasi-interpolation operator is the so-called “local-best” approximation. For $K \in \mathcal{T}_h$ and $v \in H^1(K)$, we denote by $\pi_K^p v \in \mathcal{P}_p(K)$ the $H^1$-orthogonal projection of $v$. Specifically, it is defined as the only element of $\mathcal{P}_p(K)$ such that $(\pi_K^p v, 1)_K = (v, 1)_K$ and

$$ (3.2a) \qquad (\nabla(\pi_K^p v), \nabla q_h)_K = (\nabla v, \nabla q_h)_K \quad \forall q_h \in \mathcal{P}_p(K). $$

The “local-best” approximation of $v \in H^1(\mathcal{T}_h)$, $\pi_h^p v \in \mathcal{P}_p(\mathcal{T}_h)$, is then defined by setting

$$ (3.2b) \qquad (\pi_h^p v)|_K := \pi_K^p(v|_K) \quad \forall K \in \mathcal{T}_h. $$

Notice that then

$$ \|\nabla_h(v - \pi_h^p v)\|_K = \min_{v_h \in \mathcal{P}_p(K)} \|\nabla_h(v - v_h)\|_K. $$

and

$$ \|\nabla_h(v - \pi_h^p v)\|_{\Omega} = \min_{v_h \in \mathcal{P}_p(\mathcal{T}_h)} \|\nabla_h(v - v_h)\|_{\Omega}. $$

We stress that, in general, $\pi_h^p v$ is discontinuous, i.e., $\pi_h^p v \notin H_0^1(\Omega)$ even if $v \in H_0^1(\Omega)$.

3.3. Main idea of the construction of the quasi-interpolation operator. Here, we explain the spirit of the definition of our operator. Given $u \in H_0^1(\Omega)$, we start by considering the local-best approximation $\pi_h^p u \in \mathcal{P}_p(\mathcal{T}_h)$ of $u$ given by (3.2). This function is locally defined and exhibits the best-possible approximation properties, but it is unfortunately nonconforming since $\pi_h^p u \notin H_0^1(\Omega)$ in general.

The next stage is to then to somehow “locally project” $\pi_h^p u$ to the conforming Lagrange finite element space $\mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$. The optimal projection

$$ (3.3) \qquad \min_{v_h \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)} \|\nabla_h(\pi_h^p u - v_h)\|_{\Omega} = \min_{v_h \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)} \|\nabla(u - v_h)\|_{\Omega} $$

is not satisfactory for a quasi-interpolation operator since it is global: it requires a global system solve and does not lead to local approximation properties (the minimizer from (3.3) restricted to $K \in \mathcal{T}_h$ depends on the values of $u$ in the whole domain $\Omega$ in general). The key idea is therefore to first localize $\pi_h^p u$ and then project it locally in vertex patches.

We employ the partition of unity by the hat functions (2.3) and try to approximate $\psi_a \pi_h^p u$ into a local contribution $s_h^a \in \mathcal{P}_p(\mathcal{T}_a) \cap H_0^1(\omega_a)$. This choice makes sense, since (a) $\psi_a \pi_h^p u$ is expected to be close to $\psi_a u$ which sits in $H_0^1(\omega_a)$ and (b) the boundary conditions ensure that the object resulting from the summation of the local contributions is globally conforming.

At this point, it is therefore tempting to define

$$ (3.4) \qquad s_h^a := \arg \min_{v_h \in \mathcal{P}_p(\mathcal{T}_a) \cap H_0^1(\omega_a)} \|\nabla_h(\psi_a \pi_h^p u - v_h)\|_{\omega_a}, \quad \mathcal{J}_h^p(u) := \sum_{\mathbf{a} \in \mathcal{V}_h} s_h^a. $$
The construction in (3.4) faithfully conveys the ideas behind the construction of our projector. However, as we shall see later, (3.4) would fail because we are trying to represent a polynomial of degree $p+1$ with a polynomial of degree $p$. Indeed, the presence of the hat function $\psi_a$ increases the polynomial degree by 1, and to bring it back, we will use the elementwise Lagrange interpolation operator (3.1).

### 3.4. Main idea for the computable approximation constant.
The local problem (3.4) is an approximation of the piecewise polynomial but discontinuous (nonconforming) datum $\psi_a \pi_h^p u$ by the conforming finite element method. Our key idea to obtain a computable error bound is to employ duality, i.e., define a local counterpart of problem (3.4) via Raviart-Thomas finite elements. This will allow to measure the nonconformity of the local-best approximation $\pi_h^p u$ vertex patch by vertex patch via a local discrete eigenvalue problem and bring us to the precision of the $H^1$-orthogonal projection of (3.2). This projection is defined separately on each mesh element and its approximation bound comes with a fully explicit, computable, constant, see (4.8) below. Combining these ingredients then gives estimates optimal in the mesh size $h$ for both the $H^1$ seminorm and the $L^2$ norm with a fully computable constant.

## 4. MAIN RESULTS

In this section, we construct our quasi-interpolation operator and present its fully computable error bounds. The proofs are postponed until Section 5.

### 4.1. Local potential reconstruction.
Let a vertex $a \in \mathcal{V}_h$ be given. We develop further (3.4). As it will become apparent below, we need here to start from an arbitrary piecewise polynomial $u_h \in \mathcal{P}_p(\mathcal{T}_a)$. We define $s_h^a(u_h)$ as the continuous piecewise polynomial with vanishing trace on the boundary $\partial\omega_a$ of the vertex patch subdomain $\omega_a$, $s_h^a(u_h) \in \mathcal{P}_p(\mathcal{T}_a) \cap H_0^1(\omega_a)$, by

$$ (4.1a) \qquad s_h^a(u_h) := \arg \min_{v_h \in \mathcal{P}_p(\mathcal{T}_a) \cap H_0^1(\omega_a)} \| \nabla_h(\mathcal{T}_h^p(\psi_a u_h) - v_h) \|_{\omega_a}. $$

Equivalently, $s_h^a$ is given by the Euler-Lagrange conditions of (4.1a), reading as: find $s_h^a(u_h) \in \mathcal{P}_p(\mathcal{T}_a) \cap H_0^1(\omega_a)$ such that

$$ (4.1b) \qquad (\nabla s_h^a(u_h), \nabla v_h)_{\omega_a} = (\nabla_h(\mathcal{T}_h^p(\psi_a u_h)), \nabla v_h)_{\omega_a} \quad \forall v_h \in \mathcal{P}_p(\mathcal{T}_a) \cap H_0^1(\omega_a). $$

Since $\|\nabla \cdot \psi_a\|_{\omega_a}$ is a norm on $H_0^1(\omega_a)$, $s_h^a(u_h)$ is indeed uniquely defined. We observe that the Lagrange interpolation operator $\mathcal{T}_h^p$ of (3.1) is employed element-wise in order to reduce the polynomial degree $p+1$ of $\psi_a u_h$ back to $p$. This is legal, since $\psi_a u_h$ is a polynomial (it would have failed if immediately applied to $\psi_a u$). We implicitly extend $s_h^a(u_h)$ by zero outside of $\omega_a$, so that $s_h^a$ is an operator from $\mathcal{P}_p(\mathcal{T}_a)$ to $\mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega)$.

### 4.2. Local flux reconstruction.
For each vertex $a \in \mathcal{V}_h$, consider the local finite element space

$$ (4.2) \qquad \mathbf{W}_h^a := \mathcal{R}\mathcal{T}_p(\mathcal{T}_a) \cap \mathbf{H}(\mathrm{div}^0, \omega_a) $$

of divergence-free $\mathbf{H}(\mathrm{div}, \omega_a)$-conforming Raviart-Thomas piecewise polynomials defined over the vertex patch subdomain $\omega_a$. For an arbitrary piecewise polynomial $u_h \in \mathcal{P}_p(\mathcal{T}_a)$, we define $\mathbf{r}_h^a(u_h)$ as the only element of $\mathbf{W}_h^a$ such that

$$ (4.3a) \qquad \mathbf{r}_h^a(u_h) := \arg \min_{w_h \in \mathbf{W}_h^a} \| \nabla_h(\psi_a u_h) - \mathbf{w}_h \|_{\omega_a}, $$

or, equivalently, $\mathbf{r}_h^a(u_h) \in \mathbf{W}_h^a$ such that

$$ (4.3b) \qquad (\mathbf{r}_h^a(u_h), \mathbf{w}_h)_{\omega_a} = (\nabla_h(\psi_a u_h), \mathbf{w}_h)_{\omega_a} \quad \forall \mathbf{w}_h \in \mathbf{W}_h^a. $$

Notice that we then have

$$ (4.4) \qquad \| \mathbf{r}_h^a(u_h) \|_{\omega_a} = \max_{\substack{w_h \in \mathbf{W}_h^a \\ \|w_h\|_{\omega_a} = 1}} (\nabla_h(\psi_a u_h), \mathbf{w}_h)_{\omega_a}. $$
From a practical standpoint, given $u_h \in \mathcal{P}_p(\mathcal{T}_a)$, $\mathbf{r}_h^a(u_h)$ can be computed as $\mathbf{r}_h^a(u_h) = \mathbf{r}_h$, where $(\mathbf{r}_h, q_h) \in \mathbf{X}_h^a \times \mathbf{Y}_h^a$ is the unique pair such that

$$ (4.5) \qquad \left\{ \begin{array}{l} (\mathbf{r}_h, \mathbf{v}_h)_{\omega_a} + (q_h, \nabla \cdot \mathbf{v}_h)_{\omega_a} = (\nabla_h(\psi_a u_h), \mathbf{v}_h)_{\omega_a}, \\ (\nabla \cdot \mathbf{r}_h, p_h)_{\omega_a} = 0, \end{array} \right. $$

for all $(\mathbf{v}_h, p_h) \in \mathbf{X}_h^a \times \mathbf{Y}_h^a$, with

$$ \mathbf{X}_h^a := \mathcal{R}\mathcal{T}_p(\mathcal{T}_a) \cap \mathcal{H}(\mathrm{div}, \omega_a), \quad \mathbf{Y}_h^a := \mathcal{P}_p(\mathcal{T}_a). $$

Indeed, although (4.3) and (4.5) are equivalent, in contrast to $\mathbf{W}_h^a$, both $\mathbf{X}_h^a$ and $\mathbf{Y}_h^a$ can be easily equipped with an explicit basis.

**4.3. A local primal-dual quotient and its computation by a local matrix eigenvalue problem.** Consider a vertex $\mathbf{a} \in \mathcal{V}_h$. For $u_h \in \mathcal{P}_p(\mathcal{T}_a)$, let $s_h^a(u_h)$ be defined by (4.1) and $\mathbf{r}_h^a(u_h)$ by (4.3). We define the primal-dual quotient

$$ (4.6) \qquad \lambda_{\mathbf{a}} := \max_{u_h \in \mathcal{P}_p(\mathcal{T}_a)} \frac{\|\nabla_h(\mathcal{I}_h^p(\psi_a u_h) - s_h^a(u_h))\|_{\omega_a}}{\|\mathbf{r}^a(u_h)\|_{\omega_a}}. $$

Theorem 5.1 below warrants that $\lambda_{\mathbf{a}}$ is finite and only depends on the polynomial degree $p$, the space dimension $d$, and the shape-regularity parameter $\kappa_{\mathcal{T}_a}$ of the patch.

The number $\lambda_{\mathbf{a}}$ can be computed through a small-size matrix eigenvalue problem. Specifically, given a basis of $\mathcal{P}_p(\mathcal{T}_a)$, since the operator

$$ u_h \in \mathcal{P}_p(\mathcal{T}_a) \rightarrow \mathcal{P}_p(\mathcal{T}_a) \ni \mathcal{I}_h^p(\psi_a u_h) - s_h^a(u_h) $$

is linear, its action may be encoded as a matrix $D \in \mathbb{R}^{|\mathcal{P}_p(\mathcal{T}_a)| \times |\mathcal{P}_p(\mathcal{T}_a)|}$. Similarly, the action of the operator $\mathbf{r}_h^a : \mathcal{P}_p(\mathcal{T}_a) \rightarrow \mathbf{W}_h^a \subset \mathbf{X}_h^a$ can be easily encoded as a simple matrix $R \in \mathbb{R}^{|\mathcal{P}_p(\mathcal{T}_a)| \times |\mathbf{X}_h^a|}$ after inverting the matrix representation of the left-hand side of (4.5). Then, if M and K respectively denote the mass and stiffness matrices of $\mathbf{X}_h^a = \mathcal{R}\mathcal{T}_p(\mathcal{T}_a) \cap \mathcal{H}(\mathrm{div}, \omega_a)$ and $\mathcal{P}_p(\mathcal{T}_a)$, $\lambda_{\mathbf{a}}$ can be obtained as the largest eigenvalue of the problem: Find $\mu \in \mathbb{R}$ and $v \in \mathbb{R}^{|\mathcal{P}_p(\mathcal{T}_a)|}$ such that

$$ D^T K D v = \mu R^T M R v. $$

**4.4. A local mesh shape-regularity characterization.** Recall the notation of Section 2.1. For a vertex $\mathbf{a} \in \mathcal{V}_h$, define the constant

$$ (4.7) \qquad \rho_{\mathbf{a}} := 1 + \frac{1}{\pi} \max_{K \in \mathcal{T}_a} \frac{h_K}{\tau_K^a}. $$

It is fully computable and only depends on the shape-regularity parameter $\kappa_{\mathcal{T}_a}$ of $\mathcal{T}_a$.

**4.5. Local-best approximation error estimate.** For the local-best approximation of (3.2), we observe that, applying the Poincaré inequality stated in (2.5) repeatedly, see [14, Corollary 12.13] or [39, Theorem 1.8.1] for details,

$$ (4.8) \qquad \| \nabla_h (v - \pi_h^p v) \|_K^2 \le (s+1)! \left( \frac{h_K}{\pi} \right)^{2s} |v|_{H^{1+s}(K)}^2 \quad \forall v \in H^{1+s}(K), \forall K \in \mathcal{T}_h $$

for all integer $0 \le s \le p$. Crucially, the constant appearing in (4.8) is fully explicit (computable).

**Remark 4.1** (Uniform-in-$s$ constant for a redefined Sobolev seminorm). If one redefines the seminorm $|v|_{H^s(\omega)}$ from (2.1) to repetitively contain the indices, starting with $|v|_{H^2(\omega)} = \sum_{i=1}^d \sum_{j=1}^d \|\partial_{x_j} \partial_{x_i} v\|_K^2$ in place of $|v|_{H^2(\omega)} = \sum_{i=1}^d \sum_{j=1}^d \|\partial_{x_j} \partial_{x_i} v\|_K^2$, then the factor $(s+1)!$ from (4.8) can be removed.
4.6. The quasi-interpolation operator $\mathcal{J}_h^p$. Let $u \in H_0^1(\Omega)$. To define our quasi-interpolation operator, we employ (4.1a) with the argument $u_h$ being the local-best approximation $\pi_h^p u$ of the target function $u$ given by (3.2), restricted to the patch subdomain $\omega_a$:

$$ (4.9) \qquad \mathcal{J}_h^p u := \sum_{a \in \mathcal{V}_h} s_h^a ((\pi_h^p u)|_{\omega_a}). $$

Since the computations of $\pi_h^p u$ and $s_h^a((\pi_h^p u)|_{\omega_a})$ amount to element-wise/patch-wise finite element solves, $\mathcal{J}_h^p u$ is computable. This operator is also local in the sense that for any element $K \in \mathcal{T}_h$, $(\mathcal{J}_h^p u)|_K$ only depends on the values of $u$ in the vertex patch subdomains $\omega_a$ for the vertices $a \in \mathcal{V}_K$.

4.7. $H^1$ seminorm error estimate for $\mathcal{J}_h^p$ with computable constants. We now state our main result on the approximation error of $\mathcal{J}_h^p u$ in the $H^1$ seminorm (energy norm). We present it in a Pythagoras form that we are lead to by the orthogonal projection property (3.2). Computable quasi-interpolation estimates follow from combination with (4.8).

**Theorem 4.2** ($H^1$ seminorm error estimate with computable constants). Let $u \in H_0^1(\Omega)$ and let the quasi-interpolation operator $\mathcal{J}_h^p$ be given by (3.2), (4.1), (4.9). Then

$$ (4.10) \qquad \begin{aligned} \|\nabla(u - \mathcal{J}_h^p(u))\|_K &\le \left\{ \|\nabla_h(u - \pi_h^p u)\|_K^2 + \left( \sum_{a \in \mathcal{V}_K} \rho_a \lambda_a \|\nabla_h(u - \pi_h^p u)\|_{\omega_a} \right)^2 \right\}^{1/2} \\ &\le \left( \frac{1}{(d+1)} + (d+1)c_K^2 \right)^{1/2} \left\{ \sum_{a \in \mathcal{V}_K} \|\nabla_h(u - \pi_h^p u)\|_{\omega_a}^2 \right\}^{1/2} \quad \forall K \in \mathcal{T}_h, \end{aligned} $$

with the computable constant

$$ c_K := \max_{a \in \mathcal{V}_K} (\rho_a \lambda_a). $$

There also holds

$$ (4.11) \qquad \|\nabla(u - \mathcal{J}_h^p u)\|_{\Omega} \le (1 + (d+1)^2 c_{\Omega}^2)^{1/2} \|\nabla_h(u - \pi_h^p u)\|_{\Omega} $$

with the computable constant

$$ c_{\Omega} := \max_{a \in \mathcal{V}_h} (\rho_a \lambda_a). $$

4.8. $L^2$ error estimate for $\mathcal{J}_h^p$ with computable constants. We state here our main result on the approximation error of $\mathcal{J}_h^p u$ in the $L^2$ norm. We present it in a triangle-inequality form, since there is no orthogonal projection here. Computable quasi-interpolation estimates follow from combination with (4.8).

**Theorem 4.3** ($L^2$ error estimate with computable constants). Let $u \in H_0^1(\Omega)$ and let the quasi-interpolation operator $\mathcal{J}_h^p$ be given by (3.2), (4.1), (4.9). Then

$$ (4.12) \qquad \|u - \mathcal{J}_h^p u\|_K \le \frac{h_K}{\pi} \|\nabla_h(u - \pi_h^p u)\|_K + \frac{2}{d} \sum_{a \in \mathcal{V}_K} \rho_a \lambda_a h_a \|\nabla_h(u - \pi_h^p u)\|_{\omega_a} $$

$$ (4.13) \qquad \le \left( \frac{1}{\pi} + \frac{2c_K}{d} \right) \sum_{a \in \mathcal{V}_K} h_a \|\nabla_h(u - \pi_h^p u)\|_{\omega_a} \quad \forall K \in \mathcal{T}_h. $$

There also holds

$$ (4.14) \qquad \|u - \mathcal{J}_h^p u\|_{\Omega} \le \left( \frac{1}{\pi \sqrt{d+1}} + \frac{2}{d} \sqrt{d+1} c_{\Omega} \right) \left\{ \sum_{a \in \mathcal{V}_h} h_a^2 \|\nabla_h(u - \pi_h^p u)\|_{\omega_a}^2 \right\}^{1/2}. $$
4.9. **Projection property.** From (4.11) or (4.14), when $u \in \mathcal{P}_p(\mathcal{T}_h)$, then $u = \pi_h^p u$. Hence, we immediately have:

**Theorem 4.4 (Projection).** The quasi-interpolation operator $\mathcal{J}_h^p$ from (3.2), (4.1), (4.9) is a projector, i.e.,

$$ \mathcal{J}_h^p(u) = u \quad \forall u \in \mathcal{P}_p(\mathcal{T}_h) \cap H_0^1(\Omega). $$

5. PROOFS OF THE MAIN RESULTS

5.1. **Preliminary results and notation.** For each vertex $a \in \mathcal{V}_h$, there exists a finite set $\hat{\mathcal{R}}(\kappa_{\mathcal{T}_a})$ of reference patches only depending on the shape-regularity parameter $\kappa_{\mathcal{T}_a}$ and space dimension $d$ and satisfying the following properties. Each $\hat{\mathcal{T}}_0 \in \hat{\mathcal{R}}(\kappa_{\mathcal{T}_a})$ is a conforming mesh of simplices sharing the vertex $0$, and $\hat{\psi}_0$ denotes the piecewise affine “hat function” taking value 1 at $0$ and $0$ at all the other vertices of $\hat{\mathcal{T}}_0$; $\hat{\omega}_0$ is the open domain corresponding to $\hat{\mathcal{T}}_0$. All the elements $\hat{K} \in \hat{\mathcal{T}}_0$ satisfy

$$ c(\kappa_{\mathcal{T}_a}, d) \le \rho_{\hat{K}}, \quad h_{\hat{K}} \le C(\kappa_{\mathcal{T}_a}, d), $$

for two generic constants only depending on $\kappa_{\mathcal{T}_a}$ and $d$. In addition, there exists a piecewise affine bilipchitz mapping $\phi : \hat{\omega}_0 \to \omega_a$ that transforms the elements of $\hat{\mathcal{T}}_0$ into those of $\mathcal{T}_a$. We denote by $\phi^g$ and $\phi^d$ the associated gradient- and divergence-preserving Piola mappings, $\phi^g : L^2(\hat{\omega}_0) \to L^2(\omega_a)$ and $\phi^d : L^2(\hat{\omega}_0) \to L^2(\omega_a)$ respectively defined by

$$ (5.1) \qquad \phi^g(\hat{v}) := \hat{v} \circ \phi^{-1} \quad \text{and} \quad \phi^d(\hat{w}) := \left( \frac{\mathcal{J}}{|\mathcal{J}|} \hat{w} \right) \circ \phi^{-1} $$

for all $\hat{v} \in L^2(\hat{\omega}_0)$ and $\hat{w} \in L^2(\hat{\omega}_0)$, where $\mathcal{J}$ is the Jacobian matrix of $\phi$ and $|\mathcal{J}|$ is its determinant (see, e.g., [14, Section 7.2]). The following properties of Piola mappings will be useful. First, if $\hat{v} \in H^1(\hat{\mathcal{T}}_0)$ and $\hat{w} \in L^2(\hat{\omega}_0)$, we have

$$ (5.2) \qquad (\nabla_h(\phi^g\hat{v}), \phi^d\hat{w})_{\omega_a} = (\nabla_h\hat{v}, \hat{w})_{\hat{\omega}_0}. $$

The Piola mappings are invertible, and $\hat{\psi}_0 := (\phi^g)^{-1}\psi_a$. We have $\phi^g(\hat{\psi}_0\hat{v}) = \psi_a\phi^g(\hat{v})$ and

$$ (5.3) \qquad \|\nabla_h(\phi^g\hat{v})\|_{\omega_a} \le C(\kappa_{\mathcal{T}_a}, d)h_a^{d/2-1}\|\nabla_h\hat{v}\|_{\hat{\omega}_0} $$

for all $\hat{v} \in H^1(\hat{\mathcal{T}}_0)$, with $h_a$ introduced in (2.4). Besides, $\phi^d$ is an isomorphism between $\mathcal{R}\mathcal{T}_p(\hat{\mathcal{T}}_0) \cap \mathcal{H}(\mathrm{div}^0, \hat{\omega}_0)$ and $\mathcal{R}\mathcal{T}_p(\mathcal{T}_a) \cap \mathcal{H}(\mathrm{div}^0, \omega_a)$, and we have

$$ (5.4) \qquad \|\phi^d \hat{w}\|_{\omega_a} \le C(\kappa_{\mathcal{T}_a}, d)h_a^{1-d/2}\|\hat{w}\|_{\hat{\omega}_0} \quad \forall \hat{w} \in L^2(\hat{\omega}_0). $$

We denote by $\hat{\mathcal{F}}_0$ the faces of a reference patch $\hat{\mathcal{T}}_0$. The subset

$$ \hat{\mathcal{F}}_0^i := \{ F \in \hat{\mathcal{F}}_0 \mid \hat{\psi}_0|_F \neq 0, \} $$

coinciding with interior faces for an interior vertex, will be useful. If $F \in \hat{\mathcal{F}}_0$ is a face, $\mathcal{P}_q(F)$ is the set of functions $v : F \to \mathbb{R}$ that are polynomial of degree less than or equal to $q$ in the $(d-1)$ directions tangential to $F$. For a collection of faces $\mathcal{F} \subset \hat{\mathcal{F}}_0$, we write $\mathcal{P}_q(\mathcal{F}) := \cup_{F \in \mathcal{F}} \mathcal{P}_q(F)$. We also associate with each face $F \in \hat{\mathcal{F}}_0$ a unit normal vector $n_F$. If $F \subset \partial\hat{\omega}_0$, we assume that $n_F$ points outward $\hat{\omega}_0$. For interior faces, the orientation of $n_F$ is arbitrary, but fixed. If $v \in \mathcal{P}_{p+1}(\hat{\mathcal{T}}_0)$, the jump of $v$ of through an interior face $F = \partial\hat{K}_- \cap \partial K_+ \in \hat{\mathcal{F}}_0$ is defined by

$$ \llbracket v \rrbracket_F := v_-|_F(n_- \cdot n_F) + v_+|_F(n_+ \cdot n_F) $$

where $v_\pm$ is the restriction of $v$ to $K_\pm$ and $n_\pm$ is the unit normal vector to $\partial K_\pm$ pointing outward $K_\pm$; note that $n_\pm \cdot n_F = \pm 1$ only determines the sign. For exterior faces $F \subset \partial\hat{\omega}_0$, we simply set

$$ \llbracket v \rrbracket_F = v|_F. $$

Thanks to these definitions, the integration by parts formula

$$ (5.5) \qquad (\nabla_h v, w)_{\hat{\omega}_0} = \sum_{F \in \hat{\mathcal{F}}_0} (\llbracket v \rrbracket_F, w \cdot n_F)_F - (v, \nabla \cdot w)_{\hat{\omega}_0} $$
holds true for all $w \in \mathcal{R}\mathcal{T}_p(\hat{\mathcal{T}}_0) \cap \mathcal{H}(\mathrm{div}, \hat{\omega}_0)$ and $v \in \mathcal{P}_{p+1}(\hat{\mathcal{T}}_0)$.

5.2. **Proof of the properties of $\lambda_a$ from (4.6).** Here, we establish the following important result:

**Theorem 5.1** (Upper bound on $\lambda_a$). The constant $\lambda_a$ from (4.6) is finite and only depends on the shape-regularity parameter of $\mathcal{T}_a$, the space dimension $d$, and the polynomial degree $p$, i.e.,

$$ (5.6) \qquad \lambda_a \le C(\kappa_{\mathcal{T}_a}, d, p). $$

In addition, we have

$$ (5.7) \qquad \| \nabla_h (\mathcal{T}_h^p (\psi_a u_h) - s_h^a (u_h)) \|_{\omega_a} \le \lambda_a \min_{v \in H_0^1(\omega_a)} \| \nabla_h (\psi_a u_h - v) \|_{\omega_a} $$

for all $u_h \in \mathcal{P}_p(\mathcal{T}_a)$.

Before proving Theorem 5.1, we first establish a few intermediate results. If $\hat{\mathcal{T}}_0$ is a reference patch, we use the notation $\hat{\mathcal{W}}_h^0 := \mathcal{R}\mathcal{T}_p(\hat{\mathcal{T}}_0) \cap \mathcal{H}(\mathrm{div}^0, \hat{\omega}_0)$ for the counterpart to $\mathcal{W}_h^a$ of (4.2). We first show that even though functions in $\hat{\mathcal{W}}_h^0$ satisfy a divergence constraint, we are still free to assign their normal traces on relevant faces.

**Lemma 5.2** (Traces of Raviart-Thomas elements). Let $\hat{\mathcal{T}}_0$ be a reference patch. For each $F \in \hat{\mathcal{F}}_0^1$, consider $q_F \in \mathcal{P}_p(F)$. Then, there exists $w_h \in \hat{\mathcal{W}}_h^0$ such that $w_h \cdot n_F = q_F$ for all $F \in \hat{\mathcal{F}}_0^1$.

*Proof.* We will construct a suitable function $w_h \in \hat{\mathcal{W}}_h^0$ element by element. Let $K \in \hat{\mathcal{T}}_0$ be fixed. Note that given $r_F \in \mathcal{P}_p(\mathcal{F}_K)$ and $r_K \in \mathcal{P}_p(K)$, there exists a function $w_K \in \mathcal{R}\mathcal{T}_p(K)$ such that $w_K \cdot n_K = r_F$ on $\partial K$ and $\nabla \cdot w_K = r_K$ in $K$ if and only if the Neumann compatibility condition $(r_K, 1)_K = (r_F, 1)_{\partial K}$ is satisfied. The key observation is that there is always a face $F^* \in \mathcal{F}_K$ such that $F^* \notin \hat{\mathcal{F}}_0^1$, i.e., $\hat{\psi}_0|_F = 0$: this is the face opposite to the vertex $a$. Then, we select a boundary datum $r_F \in \mathcal{P}_p(\mathcal{F}_K)$ such that $r_F|_F = q_F$ for all faces $F$ of $K$ such that $F \in \hat{\mathcal{F}}_0^1$ and, by fixing a suitable Neumann value on $F^*$, such that $(r_F, 1)_{\mathcal{F}_K} = 0$. Finally, we are eligible to impose $w_h|_K \cdot n_K = r_F$ on $\partial K$ and $\nabla \cdot (w_h|_K) = 0$. $\square$

By carefully selecting their normal traces, we can then show that functions in $\hat{\mathcal{W}}_h^0$ (and, after mapping, in $\mathcal{W}_h^a$) can be successfully employed to measure the level of nonconformity of a piecewise polynomial function.

**Corollary 5.3** (Control of the nonconformity). For all $u_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0)$, if

$$ (\nabla_h(\hat{\psi}_0 u_h), v_h)_{\hat{\omega}_0} = 0 \quad \forall v_h \in \hat{\mathcal{W}}_h^0, $$

then $\hat{\psi}_0 u_h \in H_0^1(\hat{\omega}_0)$. In particular, $u_h \in H^1(\hat{\omega}_0)$.

*Proof.* Let $u_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0)$. For all $v_h \in \hat{\mathcal{W}}_h^0$, since $\nabla \cdot v_h = 0$, Stokes' formula (5.5) shows that

$$ (\nabla_h(\hat{\psi}_0 u_h), v_h)_{\hat{\omega}_0} = \sum_{F \in \hat{\mathcal{F}}_0} (\llbracket \hat{\psi}_0 u_h \rrbracket_F, v_h \cdot n_F)_F = \sum_{F \in \hat{\mathcal{F}}_0^1} (\hat{\psi}_0 \llbracket u_h \rrbracket_F, v_h \cdot n_F)_F, $$

since $\hat{\psi}_0$ is continuous and $\hat{\psi}_0|_F = 0$ for all $F \in \hat{\mathcal{F}}_0 \setminus \hat{\mathcal{F}}_0^1$. Following Lemma 5.2, we can now pick $v_h^* \in \hat{\mathcal{W}}_h^0$ such that $v_h^* \cdot n_F = \llbracket u_h \rrbracket_F$ for all $F \in \hat{\mathcal{F}}_0^1$. This gives

$$ \sum_{F \in \hat{\mathcal{F}}_0^1} \| \hat{\psi}_0^{1/2} \llbracket u_h \rrbracket_F \|_F^2 = (\nabla_h(\hat{\psi}_0 u_h), v_h^*)_{\hat{\omega}_0} = 0. $$

Since $\hat{\psi}_0 > 0$ a.e. on each $F \in \hat{\mathcal{F}}_0^1$, this shows that $\llbracket u_h \rrbracket_F = 0$ for all $F \in \hat{\mathcal{F}}_0^1$. Hence $u_h \in H^1(\hat{\omega}_0)$. There also holds $u_h = 0$ on $F \subset \partial \hat{\omega}_0$, $F \in \hat{\mathcal{F}}_0^1$ (such faces only exist when the vertex $0$ lies on the boundary of $\hat{\omega}_0$). Thus, it immediately follows that $\hat{\psi}_0 u_h \in H_0^1(\hat{\omega}_0)$. $\square$

We are now ready to establish our main result concerning $\lambda_a$.
Proof of Theorem 5.1. Fix a vertex $\mathbf{a} \in \mathcal{V}_h$, and consider the associated reference patch $\hat{\mathcal{T}}_0 \in \hat{\mathcal{R}}(\kappa_{\mathcal{T}_a})$.

If $\hat{u}_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0)$, we define $\hat{s}_h^0(\hat{u}_h) \in \mathcal{P}_p(\hat{\mathcal{T}}_0) \cap H_0^1(\hat{\omega}_0)$ by requiring that

$$ (5.8) \qquad (\nabla \hat{s}_h^0(\hat{u}_h), \nabla \hat{v}_h)_{\hat{\omega}_0} = (\nabla_h(\mathcal{I}_h^p(\hat{\psi}_0 \hat{u}_h)), \nabla \hat{v}_h)_{\hat{\omega}_0} \quad \forall \hat{v}_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0) \cap H_0^1(\hat{\omega}_0), $$

and similarly we let $\hat{\mathbf{r}}_h^0(\hat{u}_h)$ be the only element of $\hat{\mathbf{W}}_h^0$ such that

$$ (5.9) \qquad (\hat{\mathbf{r}}_h^0(\hat{u}_h), \hat{\mathbf{w}}_h)_{\hat{\omega}_0} = (\nabla_h(\hat{\psi}_0 \hat{u}_h), \hat{\mathbf{w}}_h)_{\hat{\omega}_0} \quad \forall \hat{\mathbf{w}}_h \in \hat{\mathbf{W}}_h^0. $$

These are the natural counterparts to $s_h^a(u_h)$ from (4.1) and $\mathbf{r}_h^a(u_h)$ from (4.3) on the reference patch $\hat{\mathcal{T}}_0$.

Now consider an arbitrary $u_h \in \mathcal{P}_p(\mathcal{T}_a)$. Since $s_h^a(u_h)$ is the minimizer, see (4.1a), we have

$$ \begin{aligned} \| \nabla_h(\mathcal{I}_h^p(\psi_a u_h) - s_h^a(u_h)) \|_{\omega_a} &\le \| \nabla_h(\mathcal{I}_h^p(\psi_a u_h) - \phi^g(\hat{v}_h)) \|_{\omega_a} \\ &\le C(\kappa_{\mathcal{T}_a}, d) h_a^{d/2-1} \| \nabla_h(\mathcal{I}_h^p(\hat{\psi}_0 \hat{u}_h) - \hat{v}_h) \|_{\hat{\omega}_0} \end{aligned} $$

for any $\hat{v}_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0) \cap H_0^1(\hat{\omega}_0)$, with $\hat{u}_h := (\phi^g)^{-1}(u_h) \in \mathcal{P}_p(\hat{\mathcal{T}}_0)$. Here, we employed the fact that $\phi^g$ preserves piecewise polynomial functions as well as point values in the Lagrange interpolation nodes (affine-equivalent degrees of freedom) and used (5.3) in the last inequality. We can minimize the last norm in the right-hand side, leading to

$$ (5.10) \qquad \| \nabla_h(\mathcal{I}_h^p(\psi_a u_h) - s_h^a(u_h)) \|_{\omega_a} \le C(\kappa_{\mathcal{T}_a}, d) h_a^{d/2-1} \| \nabla_h(\mathcal{I}_h^p(\hat{\psi}_0 \hat{u}_h) - \hat{s}_h^0(\hat{u}_h)) \|_{\hat{\omega}_0}, $$

by definition (5.8) of $\hat{s}_h^0(\hat{u}_h)$.

On the other hand, still for the chosen $u_h \in \mathcal{P}_p(\mathcal{T}_a)$, because of (5.2), for any $\hat{w}_h \in \hat{\mathbf{W}}_h^0$, we also have

$$ (5.11) \qquad (\nabla_h(\psi_a u_h), \phi^d(\hat{w}_h))_{\omega_a} = (\nabla_h(\hat{\psi}_0 \hat{u}_h), \hat{w}_h)_{\hat{\omega}_0}. $$

Moreover, as in (4.4), from (5.9), there holds

$$ (5.12) \qquad \| \hat{\mathbf{r}}_h^0(\hat{u}_h) \|_{\hat{\omega}_0} = \max_{\substack{\hat{w}_h \in \hat{\mathbf{W}}_h^0 \\ \|\hat{w}_h\|_{\hat{\omega}_0} = 1}} (\nabla_h(\hat{\psi}_0 \hat{u}_h), \hat{w}_h)_{\hat{\omega}_0}. $$

Thus, using the maximizer from (5.12), $\hat{w}_h^* \in \hat{\mathbf{W}}_h^0$ with $\|\hat{w}_h^*\|_{\hat{\omega}_0} = 1$, and employing $\phi^d(\hat{w}_h^*)/\|\phi^d(\hat{w}_h^*)\|_{\omega_a}$ in (4.4), we have

$$ \begin{aligned} \| \mathbf{r}_h^a(u_h) \|_{\omega_a} &\stackrel{(4.4)}{\ge} \frac{1}{\| \phi^d(\hat{\mathbf{w}}_h^*) \|_{\omega_a}} (\nabla_h(\psi_a u_h), \phi^d(\hat{\mathbf{w}}_h^*))_{\omega_a} \\ &\stackrel{(5.11)}{=} \frac{1}{\| \phi^d(\hat{\mathbf{w}}_h^*) \|_{\omega_a}} (\nabla_h(\hat{\psi}_0 \hat{u}_h), \hat{\mathbf{w}}_h^*)_{\hat{\omega}_0} \stackrel{(5.12)}{=} \frac{1}{\| \phi^d(\hat{\mathbf{w}}_h^*) \|_{\omega_a}} \| \hat{\mathbf{r}}_h^0(\hat{u}_h) \|_{\hat{\omega}_0}. \end{aligned} $$

Finally, we use (5.4) to establish that

$$ \| \phi^d(\hat{\mathbf{w}}_h^*) \|_{\omega_a} \le C(\kappa_{\mathcal{T}_a}, d) h_a^{1-d/2} \| \hat{\mathbf{w}}_h^* \|_{\hat{\omega}_0} = C(\kappa_{\mathcal{T}_a}, d) h_a^{1-d/2}, $$

which leads to

$$ (5.13) \qquad \frac{1}{\| \mathbf{r}_h^a(u_h) \|_{\omega_a}} \le C(\kappa_{\mathcal{T}_a}, d) \frac{h_a^{1-d/2}}{\| \hat{\mathbf{r}}_h^0(\hat{u}_h) \|_{\hat{\omega}_0}}. $$

Combining (5.10) and (5.13) gives

$$ \frac{\| \nabla_h(\mathcal{I}_h^p(\psi_a u_h) - s_h^a(u_h)) \|_{\omega_a}}{\| \mathbf{r}_h^a(u_h) \|_{\omega_a}} \le C(\kappa_{\mathcal{T}_a}, d) \frac{\| \nabla_h(\mathcal{I}_h^p(\hat{\psi}_0 \hat{u}_h) - \hat{s}_h^0(\hat{u}_h)) \|_{\hat{\omega}_0}}{\| \hat{\mathbf{r}}_h^0(\hat{u}_h) \|_{\hat{\omega}_0}} \le C(\kappa_{\mathcal{T}_a}, d) \hat{\lambda}_0, $$

where

$$ \hat{\lambda}_0 := \max_{\hat{u}_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0)} \frac{\| \nabla_h(\mathcal{I}_h^p(\hat{\psi}_0 \hat{u}_h) - \hat{s}_h^0(\hat{u}_h)) \|_{\hat{\omega}_0}}{\| \hat{\mathbf{r}}_h^0(\hat{u}_h) \|_{\hat{\omega}_0}}. $$

In other words, since the last inequality is valid for all $u_h \in \mathcal{P}_p(\mathcal{T}_a)$, we have actually shown that $\lambda_a \le C(\kappa_{\mathcal{T}_a}, d) \hat{\lambda}_0$.
We next show that $\hat{\lambda}_0$ is finite, only depending on the shape regularity $\kappa_{\mathcal{T}_a}$, the space dimension $d$, and the polynomial degree $p$. Since $\mathcal{P}_p(\hat{\mathcal{T}}_0)$ is a finite-dimensional vector space (with dimension only depending on $\kappa_{\mathcal{T}_a}$, $d$, and $p$), it suffices to show that $\|\hat{\mathbf{r}}_h^0(\hat{u}_h)\|_{\hat{\omega}_0} = 0$ implies that $\|\nabla_h(\mathcal{I}_h^p(\hat{\psi}_0\hat{u}_h) - \hat{s}_h^0(\hat{u}_h))\|_{\hat{\omega}_0} = 0$ for all $\hat{u}_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0)$. Thus, we consider $\hat{u}_h \in \mathcal{P}_p(\hat{\mathcal{T}}_0)$ such that $\|\hat{\mathbf{r}}_h^0(\hat{u}_h)\|_{\hat{\omega}_0} = 0$. Due to (5.12) and Corollary 5.3, $\hat{\psi}_0\hat{u}_h \in \mathcal{P}_{p+1}(\hat{\mathcal{T}}_0) \cap H_0^1(\hat{\omega}_0)$. Then, using the minimization property of $\hat{s}_h^0(\hat{u}_h)$ that follows from (5.8), cf. (4.1a), we have

$$ \|\nabla_h(\mathcal{I}_h^p(\hat{\psi}_0\hat{u}_h) - \hat{s}_h^0(\hat{u}_h))\|_{\hat{\omega}_0} \le \|\nabla_h(\mathcal{I}_h^p(\hat{\psi}_0\hat{u}_h) - \mathcal{I}_h^p(\hat{\psi}_0\hat{u}_h))\|_{\hat{\omega}_0} = 0, $$

since $\mathcal{I}_h^p(\hat{\psi}_0\hat{u}_h) \in \mathcal{P}_p(\hat{\mathcal{T}}_0) \cap H_0^1(\hat{\omega}_0)$. Thus, we have established (5.6).

It remains to establish the bound (5.7). We first note that

$$ \|\nabla_h(\mathcal{I}_h^p(\psi_a u_h) - s_h^a(u_h))\|_{\omega_a} \le \lambda_a \|\mathbf{r}^a(u_h)\|_{\omega_a} \quad \forall u_h \in \mathcal{P}_p(\mathcal{T}_a) $$

by the definition of $\lambda_a$. On the other hand, using (4.4), we have

$$ \begin{aligned} \|\mathbf{r}^a(u_h)\|_{\omega_a} &= \max_{\substack{w_h \in \mathbf{W}_h^a \\ \|w_h\|_{\omega_a} = 1}} (\nabla_h(\psi_a u_h), w_h)_{\omega_a} = \max_{\substack{w_h \in \mathbf{W}_h^a \\ \|w_h\|_{\omega_a} = 1}} (\nabla_h(\psi_a u_h - v), w_h)_{\omega_a} \\ &\le \|\nabla_h(\psi_a u_h - v)\|_{\omega_a} \end{aligned} $$

for all $v \in H_0^1(\omega_a)$, since we always have

$$ (\nabla_h v, w_h)_{\omega_a} = (\nabla v, w_h)_{\omega_a} = -(v, \nabla \cdot w_h)_{\omega_a} = 0 $$

for all $w_h \in \mathbf{W}_h^a$.

□

**5.3. Proof of $H^1$ seminorm error estimates of Theorem 4.2.** In this section, we now derive the local and global error estimates in the $H^1$ norm stated in Theorem 4.2. We start with a fundamental result concerning the local contributions $s_h^a$.

**Lemma 5.4 (Patch-wise estimates).** *The estimate*

$$ (5.14) \qquad \|\nabla_h(\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a((\pi_h^p u)|_{\omega_a}))\|_{\omega_a} \le \rho_a \lambda_a \|\nabla_h(u - \pi_h^p u)\|_{\omega_a} $$

holds true for each vertex $a \in \mathcal{V}_h$.

*Proof.* We start with (5.7). Next, since $\psi_a u \in H_0^1(\omega_a)$, we have

$$ \min_{v \in H_0^1(\omega_a)} \|\nabla_h(\psi_a \pi_h^p u - v)\|_{\omega_a} \le \|\nabla_h(\psi_a(\pi_h^p u - u))\|_{\omega_a}. $$

Recalling that $(\pi_h^p u, 1)_K = (u, 1)_K$ for all $K \in \mathcal{T}_a$ from (3.2), we then employ the product rule element-wise to show that

$$ \begin{aligned} \|\nabla(\psi_a(u - \pi_h^p u))\|_K &\le \|\nabla\psi_a\|_{L^\infty(K)} \|u - \pi_h^p u\|_K + \|\psi_a\|_{L^\infty(K)} \|\nabla(u - \pi_h^p u)\|_K \\ &\le \left(1 + \frac{h_K}{\pi \tau_K^a}\right) \|\nabla(u - \pi_h^p u)\|_K, \end{aligned} $$

where we employed (2.2) and (2.5). The desired result (5.14) follows by the definition of $\rho_a$ in (4.7). □

*Proof of Theorem 4.2.* From (3.2a), the Pythagoras equality yields

$$ \|\nabla(u - \mathcal{J}_h^p(u))\|_K^2 = \|\nabla(u - \pi_h^p u)\|_K^2 + \|\nabla(\pi_h^p u - \mathcal{J}_h^p(u))\|_K^2. $$

The partition of unity (2.3) and the projection and linearity properties of $\mathcal{I}_h^p$ give

$$ (5.15) \qquad \sum_{a \in \mathcal{V}_h} \mathcal{I}_h^p(\psi_a \pi_h^p u) = \mathcal{I}_h^p \left( \sum_{a \in \mathcal{V}_h} (\psi_a \pi_h^p u) \right) = \mathcal{I}_h^p(\pi_h^p u) = \pi_h^p u. $$
Thus, together with definition (4.9) and Lemma 5.4, we estimate the second term as

$$
\begin{align*}
\| \nabla (\pi_h^p u - \mathcal{J}_h^p(u)) \|_K &= \left\| \sum_{a \in \mathcal{V}_K} \nabla (\mathcal{I}_h^p (\psi_a \pi_h^p u) - s_h^a ((\pi_h^p u)|_{\omega_a})) \right\|_K \\
&\le \sum_{a \in \mathcal{V}_K} \| \nabla (\mathcal{I}_h^p (\psi_a \pi_h^p u) - s_h^a ((\pi_h^p u)|_{\omega_a})) \|_K \\
&\le \sum_{a \in \mathcal{V}_K} \| \nabla_h (\mathcal{I}_h^p (\psi_a \pi_h^p u) - s_h^a ((\pi_h^p u)|_{\omega_a})) \|_{\omega_a} \\
&\le \sum_{a \in \mathcal{V}_K} \rho_a \lambda_a \| \nabla_h (u - \pi_h^p u) \|_{\omega_a},
\end{align*}
$$

and the first inequality in (4.10) immediately follows.

We also have

$$
\left( \sum_{a \in \mathcal{V}_K} \rho_a \lambda_a \| \nabla_h (u - \pi_h^p u) \|_{\omega_a} \right)^2
\le (d+1) \sum_{a \in \mathcal{V}_K} (\rho_a \lambda_a)^2 \| \nabla_h (u - \pi_h^p u) \|_{\omega_a}^2
\le (d+1) \max_{a \in \mathcal{V}_K} (\rho_a \lambda_a)^2 \sum_{a \in \mathcal{V}_K} \| \nabla_h (u - \pi_h^p u) \|_{\omega_a}^2,
$$

which proves the second inequality in (4.10).

Finally, to establish (4.11), we treat sharply

$$
\begin{align}
\| \nabla_h (\pi_h^p u - \mathcal{J}_h^p(u)) \|_{\Omega}^2 &= \sum_{K \in \mathcal{T}_h} \| \nabla (\pi_h^p u - \mathcal{J}_h^p(u)) \|_K^2 \nonumber \\
&\le \sum_{K \in \mathcal{T}_h} (d+1) \sum_{a \in \mathcal{V}_K} \| \nabla (\mathcal{I}_h^p (\psi_a \pi_h^p u) - s_h^a ((\pi_h^p u)|_{\omega_a})) \|_K^2 \nonumber \\
&= (d+1) \sum_{a \in \mathcal{V}_h} \| \nabla_h (\mathcal{I}_h^p (\psi_a \pi_h^p u) - s_h^a ((\pi_h^p u)|_{\omega_a})) \|_{\omega_a}^2 \nonumber \\
&\le (d+1) \sum_{a \in \mathcal{V}_h} (\rho_a \lambda_a)^2 \| \nabla_h (u - \pi_h^p u) \|_{\omega_a}^2 \nonumber \\
&\le \max_{a \in \mathcal{V}_h} (\rho_a \lambda_a)^2 (d+1)^2 \| \nabla_h (u - \pi_h^p u) \|_{\Omega}^2. \tag{5.16}
\end{align}
$$

□

5.4. **Proof of L² error estimates of Theorem 4.3.** We finally establish the local and global error estimates presented in Theorem 4.3. We start by establishing a Poincaré inequality.

**Lemma 5.5** (Element-wise Poincaré). *Consider an element* $K \in \mathcal{T}_h$, *one of its vertices* $a \in \mathcal{V}_K$
*and the face* $F \in \mathcal{F}_K$ *opposite* $a$. *Then, if* $v \in H^1(K)$ *satisfies* $v = 0$ *on F, we have*

$$
(5.17) \qquad \|v\|_K \le \frac{2h_K}{d} \|\nabla v\|_K.
$$

*Proof.* Let $\mathbf{y}: \mathbf{x} \to \mathbf{x} - \mathbf{a}$. Since $\nabla \cdot \mathbf{y} = \mathbf{d}$, we have

$$
d \int_K |v|^2 = \int_K \nabla \cdot \mathbf{y} |v|^2 = \int_{\partial K} \mathbf{y} \cdot \mathbf{n} |v|^2 - 2 \int_K \mathbf{y} \cdot \nabla v v.
$$

We then observe that the boundary term vanishes as $\mathbf{y} \cdot \mathbf{n} = 0$ on all the faces $F' \in \mathcal{F}_K$ sharing
the vertex $\mathbf{a}$, and $v = 0$ on the remaining face $F$. We then use that $|\mathbf{y}| \le h_K$, and therefore

$$
d \|v\|_K^2 \le 2h_K \|\nabla v\|_K \|v\|_K
$$

using the Cauchy–Schwarz inequality. The result follows.

□
**Corollary 5.6** (Closeness to the local-best approximation in $L^2$). For all $\mathbf{a} \in \mathcal{V}_h$, we have

$$ (5.18) \qquad \|\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a(\pi_h^p u)|_{\omega_a}\|_{\omega_a} \le \frac{2}{d} h_a \|\nabla_h(\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a(\pi_h^p u)|_{\omega_a})\|_{\omega_a}. $$

In addition, the estimates

$$ (5.19) \qquad \|\pi_h^p u - \mathcal{J}_h^p u\|_K \le \frac{2}{d} \sum_{\mathbf{a} \in \mathcal{V}_K} \rho_a \lambda_a h_a \|\nabla_h(u - \pi_h^p u)\|_{\omega_a} $$

and

$$ (5.20) \qquad \|\pi_h^p u - \mathcal{J}_h^p u\|_{\Omega} \le \frac{2}{d} \sqrt{d+1} c_{\Omega} \left\{ \sum_{\mathbf{a} \in \mathcal{V}_h} h_{\mathbf{a}}^2 \|\nabla_h(u - \pi_h^p u)\|_{\omega_a}^2 \right\}^{1/2} $$

hold true.

*Proof.* The estimate in (5.18) simply follows by applying Lemma 5.5 in each element of the patch, noting that, for each $K \in \mathcal{T}_a$, both $\mathcal{I}_h^p(\psi_a \pi_h^p u)$ and $s_h^a(\pi_h^p u)|_{\omega_a}$ vanish on the face of $K$ opposite to $\mathbf{a}$, see (3.1) and (4.1).

For the second estimate (5.19), we observe from (5.15) and (4.9) that

$$ (\pi_h^p u - \mathcal{J}_h^p u)|_K = \sum_{\mathbf{a} \in \mathcal{V}_K} \{\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a(\pi_h^p u)|_{\omega_a}\}|_K $$

and therefore, using the triangle inequality together with (2.4) and (5.14),

$$ \begin{aligned} \|\pi_h^p u - \mathcal{J}_h^p u\|_K &\le \sum_{\mathbf{a} \in \mathcal{V}_K} \|\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a(\pi_h^p u)|_{\omega_a}\|_K \\ &\le \frac{2}{d} \sum_{\mathbf{a} \in \mathcal{V}_K} h_K \|\nabla_h(\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a(\pi_h^p u)|_{\omega_a})\|_K, \\ &\le \frac{2}{d} \sum_{\mathbf{a} \in \mathcal{V}_K} h_a \|\nabla_h(\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a(\pi_h^p u)|_{\omega_a})\|_{\omega_a} \\ &\le \frac{2}{d} \sum_{\mathbf{a} \in \mathcal{V}_K} h_a \rho_a \lambda_a \|\nabla_h(u - \pi_h^p u)\|_{\omega_a}. \end{aligned} $$

To establish (5.20), we estimate sharply, as in (5.16),

$$ \begin{aligned} \|\pi_h^p u - \mathcal{J}_h^p u\|_{\Omega}^2 &\le (d+1) \sum_{\mathbf{a} \in \mathcal{V}_h} \|\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a((\pi_h^p u)|_{\omega_a})\|_{\omega_a}^2 \\ &\le \frac{4}{d^2} (d+1) \sum_{\mathbf{a} \in \mathcal{V}_h} h_{\mathbf{a}}^2 \|\nabla_h(\mathcal{I}_h^p(\psi_a \pi_h^p u) - s_h^a((\pi_h^p u)|_{\omega_a}))\|_{\omega_a}^2 \\ &\le \frac{4}{d^2} (d+1) \max_{\mathbf{a} \in \mathcal{V}_h} (\rho_a \lambda_a)^2 \sum_{\mathbf{a} \in \mathcal{V}_h} h_{\mathbf{a}}^2 \|\nabla_h(u - \pi_h^p u)\|_{\omega_a}^2. \end{aligned} \tag{5.18} $$

□

*Proof of Theorem 4.3.* The estimates follow from the ones in Corollary 5.6, triangle inequalities, and the fact that

$$ \|u - \pi_h^p u\|_K \le \frac{h_K}{\pi} \|\nabla_h(u - \pi_h^p u)\|_K $$
for all $K \in \mathcal{T}_h$. For (4.14), we namely use $\|u - \mathcal{J}_h^p u\|_{\Omega} \le \|u - \pi_h^p u\|_{\Omega} + \|\pi_h^p u - \mathcal{J}_h^p u\|_{\Omega}$ together with (5.20) and

$$ \begin{aligned} \|u - \pi_h^p u\|_{\Omega}^2 &\le \sum_{K \in \mathcal{T}_h} \frac{h_K^2}{\pi^2} \|\nabla(u - \pi_h^p u)\|_K^2 = \frac{1}{d+1} \sum_{a \in \mathcal{V}_h} \sum_{K \in \mathcal{T}_h | a \in \mathcal{V}_K} \frac{h_K^2}{\pi^2} \|\nabla(u - \pi_h^p u)\|_K^2 \\ &\le \frac{1}{\pi^2(d+1)} \sum_{a \in \mathcal{V}_h} h_a^2 \|\nabla_h(u - \pi_h^p u)\|_{\omega_a}^2. \end{aligned} $$

□

## 6. NUMERICAL EXAMPLES

We present here the results of several numerical experiments illustrating the actual practical behavior of our quasi-interpolate.

### 6.1. Setting.
For different two-dimensional domains $\Omega \subset \mathbb{R}^2$, functions $u \in H_0^1(\Omega) \cap C^0(\bar{\Omega})$ with various additional regularity, and triangular meshes $\mathcal{T}_h$, we compute the (broken) local-best approximation $\pi_h^1 u \in \mathcal{P}_1(\mathcal{T}_h)$,

$$ \pi_h^1 u := \arg \min_{v_h \in \mathcal{P}_1(\mathcal{T}_h)} \|\nabla_h(u - v_h)\|_{\Omega}, $$

the global-best approximation

$$ u_h^1 := \arg \min_{v_h \in \mathcal{P}_1(\mathcal{T}_h) \cap H_0^1(\Omega)} \|\nabla(u - u_h^1)\|_{\Omega}, $$

the Lagrange interpolant $\mathcal{I}_h^1 u$ of Section 3.1, and our projection $\mathcal{J}_h^1 u$ of Section 4.6, both for the polynomial degree $p = 1$. Except from the local-best approximation, all these approximations of $u$ are conforming and sit in the finite element space $\mathcal{P}_1(\mathcal{T}_h) \cap H_0^1(\Omega)$.

We employ both structured and unstructured meshes $\mathcal{T}_h$. Our structured meshes are simply constructed by first building a Cartesian grid of squares, and then breaking each square into four triangles by connected its vertices to its barycenter. The unstructured meshes are generated by the `mmg` software package [11]. We use both uniform and adapted unstructured meshes. In the uniform case, the meshes are simply generated using the `-hmax` flag of `mmg`, whereas the local mesh size is specified through a “metic file” via the `-sol` flag for adapted meshes.

For a given mesh $\mathcal{T}_h$, $N := \dim \mathcal{P}_1(\mathcal{T}_h) \cap H_0^1(\Omega)$ denotes the number of degrees of freedom of the associated finite element space. For a uniform mesh with maximal size $h$, we have $N \sim (\text{diam}(\Omega)/h)^2$.

For each configuration of domain and function, our results are summarized in a figure with 4 panels. Panels (a) and (b) aim at comparing the quality of approximation of different projectors. Specifically, the projection errors measured in the $H^1$ seminorm and in the $L^2$ norm are respectively represented on panels (a) and (b). On the other hand, in panels (c) and (d), we measure the quality of our guaranteed upper bounds, using that the local-best approximation error $\|\nabla_h(u - \pi_h^1 u)\|_{\Omega}$ is known. On panel (c), we plot $\|\nabla(u - \mathcal{J}_h^1 u)\|_{\Omega}$ and $\|u - \mathcal{J}_h^1 u\|_{\Omega}$ together with the right-hand sides in the guaranteed estimations in (4.11) and (4.14), written as

$$ \|\nabla(u - \mathcal{J}_h^1 u)\|_{\Omega} \le (1 + (d+1)^2 c_{\Omega}^2)^{1/2} \|\nabla_h(u - \pi_h^1 u)\|_{\Omega} := \eta_{H^1}(\pi_h^1 u) $$

and

$$ \|u - \mathcal{J}_h^1 u\|_{\Omega} \le \left( \frac{1}{\pi \sqrt{d+1}} + \frac{2}{d} \sqrt{d+1} c_{\Omega} \right) \left\{ \sum_{a \in \mathcal{V}_h} h_a^2 \|\nabla_h(u - \pi_h^1 u)\|_{\omega_a}^2 \right\}^{1/2} =: \eta_{L^2}(\pi_h^1 u). $$

Finally, we plot the ratios on panel (d), thereby measuring the overestimation factor in our guaranteed bounds.

We obtain similar results in all the test cases below, so that we give a discussion in Section 6.5.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Smooth solution example</figcaption></figure>

6.2. Smooth function. We start with the square $\Omega := (-1, 1)^2$ and the smooth function

$$u(x) = \sin(\pi x_1) \sin(\pi x_2).$$

We employ uniform unstructured meshes generated by mmg. The result are presented in Figure 1.

6.3. Circular irregularity. We keep the domain $\Omega := (-1, 1)^2$ and consider the function

$$u(x) = (1 - |x|)\mathbf{1}_{|x|<1}.$$

This function is $C^0$, but not $C^1$. In fact, the gradient has a line of discontinuity across the unit circle, so that $h^{1/2}$ and $h^{3/2}$ (i.e. $N^{-1/4}$ and $N^{-3/4}$) rates are expected for the $H^1$ and $L^2$ errors. We employ uniform structured meshes. Figure 2 summarizes the results.

6.4. Corner singularity. We now consider the L-shape domain $\Omega := (-1, 1)^2 \setminus [0, 1]^2$, and the standard corner singularity function

$$u(x) = \chi(x)|x|^{\alpha} \sin(\alpha\theta(x))$$

where $\alpha := 2/3$ and $\theta(x)$ is an angular variable defined such that $\theta = 0$ for $x \in (0, 1) \times \{0\}$ and $3\pi/2$ when $x \in \{0\} \times (0, 1)$ and

$$\chi(x) := (1 - |x|^2)\mathbf{1}_{|x|<1}$$

is a cutoff function that enforces boundary conditions. We first use uniform meshes generated by mmg and present the corresponding results in Figure 3. In this case, we expect and observe convergence rates of order $N^{-1/3}$ and $N^{-5/6}$ in the $H^1$ and $L^2$ norm. We then employed adapted
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Circular interface solution example</figcaption></figure>

meshes (still generated by mmg) with corresponding results presented in Figure 4. Optimal convergence rates are expected and observed in this case. The adapted meshes are generated by requiring that $h_K \lesssim \max(h_{\max}^{3/2}, |x_K|^{1/3}h_{\max})$ for all $K \in \mathcal{T}_h$, where $x_K$ is the barycenter of K.

**6.5. Discussion.** In all the benchmarks, we obtain the expected rates of convergence for our quasi-interpolation operator $\mathcal{J}_h^1$. In the $H^1$ seminorm, the errors of the global-best approximation $u_h^1$, the Lagrange interpolant $\mathcal{I}_h^1 u$, and our quasi-interpolant $\mathcal{J}_h^1 u$ are always similar. In the $L^2$ norm, our quasi-interpolant $\mathcal{J}_h^1 u$ is always more accurate than the Lagrange interpolant $\mathcal{I}_h^1 u$ and the global-best approximation $u_h^1$, and sometimes significantly so. In both the $H^1$ seminorm and the $L^2$ norm, the local-best approximation provides better errors than the conforming approximation, but the improvement is typically marginal. Our guaranteed upper bounds never underestimate the errors as predicted by the theory of Theorems 4.2 and 4.3, and behave with the correct rates. They allow for error certification. The overestimation factor is remarkably stable in all tested situations, with values around 10 in the $H^1$ seminorm and 100 in the $L^2$ norm.

## REFERENCES

1. R. Arcangeli and J.L. Gout, *Sur l'évaluation de l'erreur d'interpolation de Lagrange dans un ouvert de $\mathbb{R}^n$*, R.A.I.R.O. Analyse numérique **10** (1976), no. 3, 5-27.

2. M. Aurada, M. Feischl, J. Kemetmüller, M. Page, and D. Praetorius, *Each $H^{1/2}$-stable projection yields convergence and quasi-optimality of adaptive FEM with inhomogeneous Dirichlet data in $\mathbb{R}^d$, ESAIM Math. Model. Numer. Anal.* **47** (2013), no. 4, 1207-1235.
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Singular function example</figcaption></figure>

3. M. Bebendorf, *A note on the Poincaré inequality for convex domains*, Z. Anal. Anwendungen **22** (2003), 751–756.

4. C. Carstensen, *Quasi-interpolation and a posteriori error analysis in finite element methods*, ESAIM Math. Model. Numer. Anal. **33** (1999), no. 6, 1187–1202.

5. C. Carstensen and S.A. Funken, *Constants in Clément-interpolation error and residual based a posteriori error estimates in finite element methods*, East-West J. Numer. Math. **8** (2000), no. 3, 153–175.

6. C. Carstensen, J. Gedicke, and D. Rim, *Explicit error estimates for Courant, Crouzeix-Raviart and Raviart-Thomas finite element methods*, J. Comput. Math. **30** (2012), no. 4, 337–353.

7. C. Carstensen, D. Peterseim, and M. Schedensack, *Comparison results of finite element methods for the Poisson model problem*, SIAM J. Numer. Anal. **50** (2012), no. 6, 2803–2823.

8. T. Chaumont-Frelet and M. Vohralík, *A stable local commuting projector and optimal hp approximation estimates in H(curl)*, Numer. Math. **156** (2024), no. 6, 2293–2342.

9. P.G. Ciarlet, *The finite element method for elliptic problems*, SIAM, 2002.

10. Ph. Clément, *Approximation by finite element functions using local regularization*, ESAIM Math. Model. Numer. Anal. **9** (1975), no. R2, 77–84.

11. C. Dobrzynski, *MMG3D: User guide*, Tech. Report 422, Inria, 2012.

12. A. Ern, T. Gudi, I. Smears, and M. Vohralík, *Equivalence of local- and global-best approximations, a simple stable local commuting projector, and optimal hp approximation estimates in H(div)*, IMA J. Numer. Anal. **42** (2022), no. 2, 1023–1049.

13. A. Ern and J.-L. Guermond, *Finite element quasi-interpolation and best approximation*, ESAIM Math. Model. Numer. Anal. **51** (2017), 1367–1385.

14. A. Ern and J.-L. Guermond, *Finite Elements I. Approximation and Interpolation*, Texts in Applied Mathematics, vol. 72, Springer International Publishing, Springer Nature Switzerland AG, 2021.
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Singular function example with adapted meshes</figcaption></figure>

<footnote n="15">15. A. Ern and M. Vohralík, *Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations*, SIAM J. Numer. Anal. **53** (2015), no. 2, 1058–1081.</footnote>

<footnote n="16">16. R.S. Falk and R. Winther, *Local bounded cochain projections*, Math. Comp. **83** (2014), no. 290, 2631–2656.</footnote>

<footnote n="17">17. E. Gawlik, M.J. Holst, and M.W. Licht, *Local finite element approximation of Sobolev differential forms*, ESAIM Math. Model. Numer. Anal. **55** (2021), no. 5, 2075–2099.</footnote>

<footnote n="18">18. Tingting Hao, Xiaofei Guan, Shipeng Mao, and Shaochun Chen, *Computable interpolation error constants for the geometric simplex finite elements*, J. Sci. Comput. **87** (2021), no. 1, Paper No. 35, 21.</footnote>

<footnote n="19">19. R. Hiptmair and C. Pechstein, *Discrete regular decompositions of tetrahedral discrete 1-forms*, in collection *Maxwell's Equations*, De Gruyter, 2019, pp. 199–258.</footnote>

<footnote n="20">20. F. Kikuchi and X. Liu, *Estimation of interpolation error constants for the $P_0$ and $P_1$ triangular finite elements*, Comput. Methods Appl. Mech. Engrg. **196** (2007), no. 37-40, 3750–3758.</footnote>

<footnote n="21">21. K. Kobayashi and T. Tsuchiya, *A Babuška-Aziz proof of the circumradius condition*, Japan J. Indust. Appl. Math. **31** (2014), 193–210.</footnote>

<footnote n="22">22. Kenta Kobayashi, *On the interpolation constants over triangular elements*, Applications of mathematics 2015, Czech. Acad. Sci., Prague, 2015, pp. 110–124.</footnote>

<footnote n="23">23. M.W. Licht, *Averaging-based local projections in finite element exterior calculus*, arXiv preprint arXiv:2301.03007, 2023.</footnote>

<footnote n="24">24. X. Liu and F. Kikuchi, *Analysis and estimation of error constants for $P_0$ and $P_1$ interpolations over triangular finite elements*, J. Math. Sci. Univ. Tokyo **17** (2010), 27–78.</footnote>

<footnote n="25">25. X. Liu and C. You, *Explicit bound for quadratic Lagrange interpolation constant on triangular finite elements*, Appl. Math. Comput. **319** (2018), 693–701.</footnote>

<footnote n="26">26. A. Målqvist and D. Peterseim, *Localization of elliptic multiscale problems*, Math. Comp. **83** (2014), no. 290, 2583–2603.</footnote>
27. J.M. Melenk, hp-interpolation of nonsmooth functions and an aplication to hp-a posteriori error estimation, SIAM J. Numer. Anal. **43** (2005), no. 1, 127–155.

28. J.C. Nédélec, *Mixed finite elements in R³*, Numer. Math. **35** (1980), 315–341.

29. P. Oswald, *On a BPX-preconditioner for P1 elements*, Computing **51** (1993), no. 2, 125–133.

30. L.E. Payne and H.F. Weinberger, *An optimal Poincaré inequality for convex domains*, Arch. Ration. Mech. Anal. **5** (1960), 1005–1036.

31. P.A. Raviart and J.M. Thomas, *A mixed finite element method for 2nd order elliptic problems*, Mathematical Aspects of Finite Element Methods, Springer-Verlag, 1977.

32. J. Schöberl, *Commuting quasi-interpolation operators for mixed finite elements*, Tech. Report ISC-01-10-MATH, Texas A&M University, 2001.

33. J. Schöberl, J.M. Melenk, C. Pechstein, and S. Zaglmayr, *Additive Schwarz preconditioning for p-version triangular and tetrahedral finite elements*, IMA J. Numer. Anal. **28** (2008), no. 1, 1–24.

34. L.R. Scott and S. Zhang, *Finite element interpolation of nonsmooth functions satisfying boundary conditions*, Math. Comp. **54** (1990), no. 190, 483–493.

35. A. Veeser, *Approximating gradients with continuous piecewise polynomial functions*, Found. Comput. Math. **16** (2016), 723–750.

36. A. Veeser and R. Verfürth, *Explicit upper bounds for dual norms of residuals*, SIAM J. Numer. Anal. **47** (2009), no. 3, 2387–2405.

37. ——, *Poincaré constants for finite element stars*, IMA J. Numer. Anal. **32** (2012), no. 1, 30–47.

38. R. Verfürth, *Error estimates for some quasi-interpolation operators*, ESAIM Math. Model. Numer. Anal. **33** (1999), no. 4, 695–713.

39. M. Vohralík, *Analysis and approximation of partial differential equations by finite elements*, ENSTA, Institut Polytechnique de Paris, 2025, Lecture notes.

40. ——, *p-robust equivalence of global continuous and local discontinuous approximation, a p-stable local projector, and optimal elementwise hp approximation estimates in H¹*, HAL Preprint 04436063, submitted for publication, 2025.