Simple and robust equilibrated flux a posteriori estimates for singularly perturbed reaction-diffusion problems

TL;DR

This paper introduces a simple, robust a posteriori error estimator for finite element solutions to singularly perturbed reaction-diffusion problems, applicable to arbitrary-order methods without special submeshes.

Contribution

It presents a new equilibrated flux estimator that is robust and simple to implement for any finite element order, avoiding complex submesh constructions or estimator combinations.

Findings

Provides guaranteed global error bounds without unknown constants.

Achieves local robustness with respect to mesh size and perturbation parameters.

Includes explicit bounds for inverse inequality constants on simplices.

Abstract

We consider energy norm a posteriori error analysis of conforming finite element approximations of singularly perturbed reaction-diffusion problems on simplicial meshes in arbitrary space dimension. Using an equilibrated flux reconstruction, the proposed estimator gives a guaranteed global upper bound on the error without unknown constants, and local efficiency robust with respect to the mesh size and singular perturbation parameters. Whereas previous works on equilibrated flux estimators only considered lowest-order finite element approximations and achieved robustness through the use of boundary-layer adapted submeshes or via combination with residual-based estimators, the present methodology applies in a simple way to arbitrary-order approximations and does not request any submesh or estimators combination. The equilibrated flux is obtained via local reaction-diffusion problems with suitable weights (cut-off factors), and the guaranteed upper bound features the same weights. We prove that the inclusion of these weights is not only sufficient but also necessary for robustness of any flux equilibration estimate that does not employ submeshes or estimators combination, which shows that some of the flux equilibrations proposed in the past cannot be robust. To achieve the fully computable upper bound, we derive explicit bounds for some inverse inequality constants on a simplex, which may be of independent interest.

Figures

Equations

$$-{\epsilon }^2\Delta u+{\kappa }^{2}u$$

$$u$$

$$a(w,v):={\epsilon }^2(\nabla w,\nabla v)+{\kappa }^2(w,v),w,v\in H_0^1(\Omega ),$$

$$a(u,v)=(f,v)$$

$${|||v|||}^2:=a(v,v),$$

$$a(u_{\mathcal{T}},v_{\mathcal{T}})=(f,v_{\mathcal{T}})\forall v_{\mathcal{T}}\in V_{\mathcal{T}}.$$

$$|||u-u_{\mathcal{T}}|||\le C_{\mathrm{rel}}\eta (u_{\mathcal{T}}),$$

$${\eta }_{\mathrm{res}}(u_{\mathcal{T}}{)}^2:=\sum _{K\in \mathcal{T}}{\alpha }_K^2\Vert r_{\mathcal{T}}{\Vert }_K^2+\sum _{F\in {\mathcal{F}}_{\Omega }}{\epsilon }^{-1}{\alpha }_F\Vert j_{\mathcal{T}}{\Vert }_F^2,$$

$$r_{\mathcal{T}}{|}_K$$

$$j_{\mathcal{T}}{|}_F$$

$${\alpha }_S:=\min \left\{\frac{h_S}{\epsilon },\frac{1}{\kappa }\right\},$$

$$\min \left\{\Vert \epsilon \nabla u_{\mathcal{T}}+{\epsilon }^{-1}{\mathbf{\sigma }}_{\mathcal{T}}{\Vert }_K,C{\epsilon }^{-\frac{1}{2}}{\alpha }_F^{\frac{1}{2}}\Vert j_{\mathcal{T}}{\Vert }_{\partial K\setminus \partial \Omega }\right\}.$$

$$\nabla \cdot {\mathbf{\sigma }}_{\mathcal{T}}+{\kappa }^2{\varphi }_{\mathcal{T}}={\Pi }_{\mathcal{T}}f,$$

$${|||u-u_{\mathcal{T}}|||}^2\le \sum _{K\in \mathcal{T}}{\left[w_K\Vert \epsilon \nabla u_{\mathcal{T}}+{\epsilon }^{-1}{\mathbf{\sigma }}_{\mathcal{T}}{\Vert }_K+\Vert \kappa \left(u_{\mathcal{T}}-{\varphi }_{\mathcal{T}}\right){\Vert }_K+{\tilde{w}}_K\Vert f-{\Pi }_{\mathcal{T}}f{\Vert }_K\right]}^2,$$

$$w_K=\min \left\{1,C_{\ast }\sqrt{\frac{\epsilon }{\kappa h_K}}\right\},{\tilde{w}}_K=\min \left\{\frac{h_K}{\pi \epsilon },\frac{1}{\kappa }\right\},$$

$$w_K^2\Vert \epsilon \nabla u_{\mathcal{T}}+{\epsilon }^{-1}{\mathbf{\sigma }}_{\mathcal{T}}{\Vert }_K^2+\Vert \kappa \left(u_{\mathcal{T}}-{\varphi }_{\mathcal{T}}\right){\Vert }_K^2$$

$$\lesssim \sum _{K^{\prime }\in {\mathfrak{T}}_K}\left[{|||u-u_{\mathcal{T}}|||}_{K^{\prime }}^2+{\alpha }_{K^{\prime }}^2\Vert f-{\Pi }_{\mathcal{T}}f{\Vert }_{K^{\prime }}^2\right],$$

$${|||v|||}_K^2:={\epsilon }^2\Vert \nabla v{\Vert }_K^2+{\kappa }^2\Vert v{\Vert }_K^{2}K\in \mathcal{T},$$

$${\vartheta }_{\mathcal{T}}:=\underset{K\in \mathcal{T}}{\max }\frac{h_K}{{\rho }_K},$$

$${\mathbb{P}}_p(\mathcal{T}):=\{v\in L^2(\Omega ),v{|}_K\in {\mathbb{P}}_p(K)\forall K\in \mathcal{T}\}$$

$$\begin{array}{rl}{\mathbf{RTN}}_p(\mathcal{T})& :=\{{\mathbf{v}}_{\mathcal{T}}\in {\mathbf{L}}^2(\Omega ),{\mathbf{v}}_{\mathcal{T}}{|}_K\in {\mathbf{RTN}}_p(K)\forall K\in \mathcal{T}\},\\ {\mathbf{RTN}}_p(K)& :={\mathbb{P}}_p(K;{\mathbb{R}}^d)+{\mathbb{P}}_p(K)\mathbf{x}.\end{array}$$

$$\Vert v{\Vert }_{\partial K}\le C_{\mathrm{Tr}}\Vert \nabla v{\Vert }_K^{\frac{1}{2}}\Vert v{\Vert }_K^{\frac{1}{2}},$$

$$\Vert v{\Vert }_F^2\le {\vartheta }_{\mathcal{T}}\left(2\Vert \nabla v{\Vert }_K+d/h_K\Vert v{\Vert }_K\right)\Vert v{\Vert }_K.$$

$$h_K^{1/2}\Vert \mathbf{v}\cdot \mathbf{n}{\Vert }_{\partial K}$$

$$C_{p+1,d,\partial K}$$

$$C_{p+1,d,K}$$

$$C_{\ast }:=\frac{1}{\sqrt{2}}\left(\frac{1}{\sqrt{\pi }}{C}_{p+1,d,K}+C_{\mathrm{Tr}}{C}_{p+1,d,\partial K}\right).$$

$${\mathbf{V}}_{\mathcal{T}}^{\mathbf{a}}$$

$$Q_{\mathcal{T}}^{\mathbf{a}}$$

$$({\mathbf{\sigma }}_{\mathcal{T}}^{\mathbf{a}},{\varphi }_{\mathcal{T}}^{\mathbf{a}}):=\underset{\begin{array}{c}({\mathbf{v}}_{\mathcal{T}},q_{\mathcal{T}})\in {\mathbf{V}}_{\mathcal{T}}^{\mathbf{a}}\times Q_{\mathcal{T}}^{\mathbf{a}}\\ \nabla \cdot {\mathbf{v}}_{\mathcal{T}}+{\kappa }^2{q}_{\mathcal{T}}={\Pi }_{\mathcal{T}}(f{\psi }_{\mathbf{a}})-{\epsilon }^2\nabla u_{\mathcal{T}}\cdot \nabla {\psi }_{\mathbf{a}}\end{array}}{\mathrm{argmin}}{w}_{\mathbf{a}}^2\Vert \epsilon {\psi }_{\mathbf{a}}\nabla u_{\mathcal{T}}+{\epsilon }^{-1}{\mathbf{v}}_{\mathcal{T}}{\Vert }_{{\omega }_{\mathbf{a}}}^2+\Vert \kappa \left[{\Pi }_{\mathcal{T}}({\psi }_{\mathbf{a}}{u}_{\mathcal{T}})-q_{\mathcal{T}}\right]{\Vert }_{{\omega }_{\mathbf{a}}}^2$$

$$w_{\mathbf{a}}:=\min \left\{1,C_{\ast }\sqrt{\frac{\epsilon }{\kappa h_{{\omega }_{\mathbf{a}}}}}\right\}.$$

$${\mathbf{\sigma }}_{\mathcal{T}}:=\sum _{\mathbf{a}\in \mathcal{V}}{\mathbf{\sigma }}_{\mathcal{T}}^{\mathbf{a}},{\varphi }_{\mathcal{T}}:=\sum _{\mathbf{a}\in \mathcal{V}}{\varphi }_{\mathcal{T}}^{\mathbf{a}}.$$

$$({\Pi }_{\mathcal{T}}(f{\psi }_{\mathbf{a}})-{\epsilon }^2\nabla u_{\mathcal{T}}\cdot \nabla {\psi }_{\mathbf{a}},1{)}_{{\omega }_{\mathbf{a}}}=(f,{\psi }_{\mathbf{a}}{)}_{{\omega }_{\mathbf{a}}}-{\epsilon }^2(\nabla u_{\mathcal{T}},\nabla {\psi }_{\mathbf{a}}{)}_{{\omega }_{\mathbf{a}}}={\kappa }^2(u_{\mathcal{T}},{\psi }_{\mathbf{a}}{)}_{{\omega }_{\mathbf{a}}}$$

$$(\nabla \cdot {\mathbf{\sigma }}_{\mathcal{T}}^{\mathbf{a}},1{)}_{{\omega }_{\mathbf{a}}}=({\mathbf{\sigma }}_{\mathcal{T}}^{\mathbf{a}}\cdot {\mathbf{n}}_{{\omega }_{\mathbf{a}}},1{)}_{\partial {\omega }_{\mathbf{a}}}=0$$

Full text

Key words: Singular perturbation, a posteriori error analysis, local efficiency, robustness, equilibrated flux

1 Introduction

Let $\Omega$ be a polygonal/polyhedral/polytopal domain in ${\mathbb{R}}^{d}$, $d\geq 1$, with a Lipschitz-continuous boundary. Let $\varepsilon>0$ and $\kappa\geq 0$ be two fixed real parameters, and let $f\in L^{2}(\Omega)$ be a given source term. Consider the problem: find $u:\Omega\rightarrow{\mathbb{R}}$ such that

[TABLE]

Let $a(\cdot,\cdot)$ be the symmetric bilinear form defined by

[TABLE]

where $({\cdot},{\cdot})$ denotes the $L^{2}$-inner product of scalar- and vector-valued functions on $\Omega$, with associated norm $\lVert\cdot\rVert$. The restriction of the $L^{2}$-inner product to an open subset $\omega\subset\Omega$ is denoted by $({\cdot},{\cdot})_{\omega}$, with associated norm $\lVert{\cdot}\rVert_{\omega}$. The weak formulation of problem (1.1) is to find $u\in H^{1}_{0}(\Omega)$ such that

[TABLE]

The energy norm ${\lvert\kern-1.07639pt\lvert\kern-1.07639pt\lvert{\cdot}\rvert\kern-1.07639pt\rvert\kern-1.07639pt\rvert}$ associated to problem (1.1) is then the norm induced by the form $a({\cdot},{\cdot})$, namely

[TABLE]

In this paper, we shall be primarily interested in the case where $\varepsilon\ll\kappa$, when problem (1.1) is said to be singularly perturbed. Then, the accurate numerical approximation can be challenging due to the typical presence of sharp boundary and/or interior layers in the solution.

In order to present more specifically the focus of this work, let us consider a simplicial mesh $\mathcal{T}$ of $\Omega$ and let $V_{{\mathcal{T}}}\coloneqq{\mathbb{P}}_{p}(\mathcal{T})\cap H^{1}_{0}(\Omega)$ denote the subspace of $H^{1}_{0}(\Omega)$ of piecewise polynomial functions of degree at most $p$, where $p\geq 1$ is a fixed integer. The conforming Galerkin finite element approximation of (1.3) consists of finding $u_{{\mathcal{T}}}\in V_{{\mathcal{T}}}$ such that

[TABLE]

The goal is to find a computable a posteriori error estimator $\eta(u_{{\mathcal{T}}})$ that satisfies

[TABLE]

The first inequality in (1.6) is called reliability, while the second inequality is called (global) efficiency. A localized version of the efficiency bound is actually desirable. The quality of the estimator is determined by the product of the two constants $C_{\mathrm{rel}}$ and $C_{\mathrm{eff}}$. A key requirement for singularly perturbed problems is to obtain estimators that are robust in the sense that both constants $C_{\mathrm{rel}}$ and $C_{\mathrm{eff}}$ are independent of the singular perturbation parameters $\varepsilon$ and $\kappa$. Only such estimates can quantify well the error in the numerical approximation and be reliably used in adaptive algorithms which allow for efficient approximation of the localized features of the solution.

Recently, several methodologies for constructing error estimators that satisfy (1.6) in a robust way have been studied. Verfürth [35] (see also [36] or [38, Section 4.3]) was probably the first to show robust bounds, in the framework of the so-called residual-based estimates. For the problem at hand, these estimators take the form (up to the data oscillation term and possible generic constants)

[TABLE]

where the local element and face residuals are defined respectively by

[TABLE]

and where $\Delta_{\mathcal{T}}$ denotes the element-wise Laplacian, $\llbracket\nabla u_{{\mathcal{T}}}{\cdot}\bm{n}_{F}\rrbracket_{F}$ denotes the jump of the normal component of $\nabla u_{{\mathcal{T}}}$ over the face $F$, $\mathcal{F}_{\Omega}$ stands for the set of internal faces of the mesh $\mathcal{T}$, and the weights (cut-off factors) take the form

[TABLE]

with $h_{S}$ being the diameter of $S$, where $S$ is either a simplex ${K}$ or a face $F$. The resulting estimator $\eta_{\mathrm{res}}(u_{{\mathcal{T}}})$ is thus a straightforward extension from the pure diffusion case $\kappa=0$ and is simple to implement in practice. The proof that $\eta_{\mathrm{res}}$ satisfies the second inequality in (1.6) rests on a bubble function technique, where the face bubble functions are defined with respect to a submesh matching the boundary-layer length scales and are possibly very steeply decaying. Their role is to capture the sharp layers caused by the singular perturbation. Note that these bubble functions, and hence the submeshes on which they are defined, are only employed in the analysis; thus they do not need to be constructed in practice. Shortly after, Ainsworth and Babuška [2] extended the method of equilibrated residuals, cf. [3], to satisfy (1.6) in a robust way for lowest-order approximations, i.e. $p=1$. In contrast to the residual-based estimators, a boundary-layer adapted submesh in each mesh element needs to be constructed in practice in order to evaluate the estimator.

Further progress has been made since, although, to the best of our knowledge, only in the case of lowest-order approximations where the polynomial degree $p=1$. Robust estimates that are guaranteed ($C_{\mathrm{rel}}=1$) and where $\eta(u_{{\mathcal{T}}})$ is fully computable have been obtained in Cheddadi et al. [9]. This remedies that $C_{\mathrm{rel}}$ is unknown for residual-based estimates and that exact solutions of some infinite-dimensional boundary value problems on each element (which cannot be performed exactly in practice) are required in the equilibrated residuals approach. The estimator in [9] is based on an equilibrated flux $\bm{\sigma}_{{\mathcal{T}}}$ belonging to a discrete subspace of $\bm{H}(\operatorname{div})$ that satisfies the equilibration identity $\nabla{\cdot}\bm{\sigma}_{{\mathcal{T}}}+\kappa^{2}u_{{\mathcal{T}}}=f_{{\mathcal{T}}}$, where $f_{{\mathcal{T}}}$ is a piecewise polynomial approximation of $f$. The estimator is then composed of terms of the form

[TABLE]

Thus it can be seen as a combination between an equilibrated flux estimator for diffusion problems and the residual-based estimator of [35] for reaction–diffusion problems. No submesh is needed for the construction of the estimator. Subsequently, Ainsworth and Vejchodský [4, 5] proceed in two stages. First, equilibrated face fluxes are computed as in [2], and then, equilibrated fluxes are obtained by face liftings, so that the final estimate $\eta(u_{{\mathcal{T}}})$ is also fully computable and the first inequality in (1.6) is guaranteed with $C_{\mathrm{rel}}=1$. As in [2], though, boundary-layer adapted submeshes appear in the construction of the estimator.

The use of a submesh complicates the construction and implementation of the equilibrated flux estimators of [4, 5]. Moreover, it is likely to be even more involved when moving beyond lowest-order approximations. In this work, by further developing the idea in [9], we show how to obtain simple, i.e. avoiding any submesh, yet robust equilibrated flux estimators for arbitrary-order approximations. The a posteriori error estimates presented in this paper are based on a locally computable flux $\bm{\sigma}_{{\mathcal{T}}}$ and potential approximation $\phi_{{\mathcal{T}}}$, respectively belonging to discrete subspaces of $\bm{H}(\mathrm{div},\Omega)$ and $L^{2}(\Omega)$ of the current mesh $\mathcal{T}$, that satisfy the key equilibration property

[TABLE]

where $\Pi_{{\mathcal{T}}}\colon L^{2}(\Omega)\rightarrow{\mathbb{P}}_{p}(\mathcal{T})$ denotes the $L^{2}$-orthogonal projection operator. The upper bound on the error then has the simple form

[TABLE]

where $w_{{K}}$ is an elementwise computable weight (cut-off factor) such that

[TABLE]

with a fixed computable constant $C_{*}$ given by (2.7); see Theorem 3.1 below for further details. The equilibrated flux $\bm{\sigma}_{{\mathcal{T}}}$ and approximate potential $\phi_{{\mathcal{T}}}$ in (1.10), (1.11) are obtained by an extension of the patchwise equilibration of [12, 8], see also [7, 19].

Furthermore, we prove robustness and efficiency of the estimator (1.11) by showing that its local contributions are bounded, up to a constant, by the local residual estimators. More precisely, for each ${K}\in\mathcal{T}$, we show that

[TABLE]

where $\mathfrak{T}_{K}$ and $\mathfrak{F}_{K}$ denote the set of elements and faces in a suitable neighbourhood of ${K}$ and

[TABLE]

see Proposition 4.3 and Theorem 4.4 below for full details. Crucially, the constants hidden in $\lesssim$ in (1.12) are independent of the mesh-sizes $h_{K}$ and problem parameters $\varepsilon$ and $\kappa$, depending only on the shape-regularity of $\mathcal{T}$, the space dimension $d$, and the polynomial degree $p$. Hence, just as for residual-based estimates, equilibrated flux estimates have a straightforward extension from the pure diffusion case $\kappa=0$, based on including appropriate weights (cut-off factors) and not requiring computations of quantities over any submesh or combination with the residual estimators. In light of these results, we believe that the claims in [37, 38] of a “structural defect” of the robustness of the equilibrated fluxes estimators are not generally valid.

As a side result, we also prove in Proposition 5.1 that the weights $w_{{K}}$ in (1.11) are necessary for robustness of any equilibrated flux estimate involving the terms $\lVert\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}{\bm{\sigma}}_{\mathcal{T}}\rVert_{K}$ whenever $\bm{\sigma}_{{\mathcal{T}}}$ is a piecewise polynomial on $\mathcal{T}$ (and thus its construction does not involve any submesh), regardless of the precise details of the construction of $\bm{\sigma}_{{\mathcal{T}}}$. This proves that several flux equilibrations proposed in the past cannot be robust with respect to reaction dominance in general (although in many constellations, no loss of robustness may be numerically observed), including those of Repin and Sauter [29], Ainsworth et al. [1], Eigel and Samrowski [16], Eigel and Merdon [15], and Vejchodský [32, 34, 33].

We only treat isotropic meshes. Results for anisotropic meshes can be found in Kunert [27], Grosman [22], Apel et al. [6], Zhao and Chen [40], or Kopteva [24, 25]. Also, we are solely interested in the energy norm. Robust estimates in the maximum norm are obtained in Demlow and Kopteva [11] and, on possibly anisotropic meshes, in Kopteva [23] for $p=1$ any in Linss [28] for any order $p\geq 1$ in one space dimension. We refer to Stevenson [31] for robust convergence, and we refer to Faustmann and Melenk [21] and the references therein for balanced norms. Finally, extensions to variable coefficients $\varepsilon$ and $\kappa$ can be treated easily as in [5], whereas inhomogeneous Dirichlet and Neumann boundary conditions, mixed parallelepipedal–simplicial meshes, meshes with hanging nodes, and approximations with varying polynomial degree $p$ can be treated as in Dolejší et al. [14].

2 Construction of the equilibrated flux

We present in this section the construction of our equilibrated flux $\bm{\sigma}_{{\mathcal{T}}}$ and of the potential approximation $\phi_{{\mathcal{T}}}$.

2.1 Notation

Let $\mathcal{T}$ be a matching simplicial partition of the domain $\Omega$, i.e., $\bigcup_{{K}\in\mathcal{T}}{K}=\overline{\Omega}$, any element ${K}\in\mathcal{T}$ is a closed simplex (interval when $d=1$, triangle when $d=2$, tetrahedron when $d=3$), and the intersection of two different simplices is either empty, or a vertex, or their common $l$-dimensional face, $1\leq l\leq d-1$. We denote by $\sigma_{\mathcal{T}}>0$ the shape-regularity parameter of the mesh $\mathcal{T}$, i.e.

[TABLE]

where $\rho_{K}$ is the diameter of the largest ball contained in $K$. For each element ${K}\in\mathcal{T}$ and for a fixed integer $p\geq 1$, let ${\mathbb{P}}_{p}({K})$ denote the space of polynomials of total degree at most $p$ on ${K}$. Let

[TABLE]

denote the space of scalar piecewise polynomials of degree at most $p$ over $\mathcal{T}$. Let $\Pi_{{\mathcal{T}}}\colon L^{2}(\Omega)\rightarrow{\mathbb{P}}_{p}(\mathcal{T})$ denote the $L^{2}$-orthogonal projection operator from $L^{2}(\Omega)$ onto ${\mathbb{P}}_{p}(\mathcal{T})$. We additionally consider $\bm{L}^{2}(\Omega)\coloneqq L^{2}(\Omega;\mathbb{R}^{d})$ and $\bm{RTN}_{p}(\mathcal{T})\subset\bm{L}^{2}(\Omega)$ the piecewise Raviart–Thomas–Nédélec space defined by

[TABLE]

For any subset $S$ of $\overline{\Omega}$, let $h_{S}$ denote the diameter of $S$. Thus, for instance, $h_{K}$ denotes the diameter of the element ${K}\in\mathcal{T}$. Let $\mathcal{V}$ denote the set of vertices of the mesh $\mathcal{T}$. It is partitioned into the set of interior vertices $\mathcal{V}^{\mathrm{int}}\coloneqq\{{\bm{a}}\in\mathcal{V},\quad{\bm{a}}\in\Omega\}$, and boundary vertices $\mathcal{V}^{\mathrm{ext}}\coloneqq\mathcal{V}\setminus\mathcal{V}^{\mathrm{int}}$. For each vertex ${\bm{a}}\in\mathcal{V}$, the function $\psi_{{\bm{a}}}$ is the hat function associated with ${\bm{a}}$, i.e., $\psi_{{\bm{a}}}\in{\mathbb{P}}_{1}(\mathcal{T})\cap H^{1}(\Omega)$ taking value $1$ in the vertex ${\bm{a}}$ and [math] in the other vertices. The set ${\omega_{\bm{a}}}$ is the interior of the support of $\psi_{{\bm{a}}}$ with associated diameter $h_{{\omega_{\bm{a}}}}$. Furthermore, let ${\mathcal{T}_{{\bm{a}}}}$ denote the restriction of the mesh $\mathcal{T}$ to ${\omega_{\bm{a}}}$, and let $\mathcal{F}_{{\bm{a}}}$ denote the set of interior faces of ${\mathcal{T}_{{\bm{a}}}}$, i.e. the faces of ${\mathcal{T}_{{\bm{a}}}}$ that contain the vertex ${\bm{a}}$ for ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$, without those on $\partial\Omega$ for ${\bm{a}}\in\mathcal{V}^{\mathrm{ext}}$. For each element ${K}\in\mathcal{T}$, we collect in $\mathcal{V}_{K}$ the set of vertices of $\mathcal{V}$ belonging to ${K}$. We also define $\mathfrak{T}_{K}\coloneqq\bigcup_{{\bm{a}}\in\mathcal{V}_{K}}{\mathcal{T}_{{\bm{a}}}}$ and $\mathfrak{F}_{K}\coloneqq\bigcup_{{\bm{a}}\in\mathcal{V}_{K}}\mathcal{F}_{{\bm{a}}}$.

Throughout this work, the notation $a\lesssim b$ means that $a\leq Cb$ with a constant $C$ that only depends on the shape-regularity parameter $\vartheta_{\mathcal{T}}$ of $\mathcal{T}$, on the space dimension $d$, and on the polynomial degree $p$, so that it is in particular independent of the mesh-sizes $h_{{K}}$ and of the problem parameters $\varepsilon$ and $\kappa$; $a\simeq b$ then stands for $a\lesssim b$ and simultaneously $b\lesssim a$.

2.2 Trace and inverse inequalities

We first recall two inequalities that we will rely on.

Lemma 2.1** (Trace inequality with explicit constant).**

For all ${K}\in\mathcal{T}$ and for all $v\in H^{1}({K})$ that satisfy $(v,1)_{{K}}=0$, i.e., that have vanishing mean-value on ${K}$, there holds

[TABLE]

Proof.

We refer the reader to [13, Lemma 1.49] for the explicit constants of the trace inequality for general functions in $H^{1}({K})$; namely, for each face $F\subset\partial{K}$,

[TABLE]

Then, we additionally apply the Poincaré inequality $\lVert v\rVert_{{K}}\leq h_{K}/\pi\lVert\nabla v\rVert_{{K}}$ for functions with vanishing mean-value on ${K}$, and sum over all the faces $F$ to obtain (2.3). ∎

Lemma 2.2** (Inverse inequalities with explicit constants).**

For any $K\in\mathcal{T}$ and any $\bm{v}\in\bm{RTN}_{p}(K)$, we have

[TABLE]

where the constants $C_{p+1,d,\partial K}$ and $C_{p+1,d,K}$ are given by

[TABLE]

Proof.

See Appendix A. ∎

In practice, possibly sharper constants can be obtained for the inequalities in (2.4) by solving numerically small eigenvalue problems on each mesh element, or on a reference element in combination with bounds for the influence of the affine mapping.

We will need below the following constant composed of the constants of the trace and inverse inequalities (2.3) and (2.4):

[TABLE]

2.3 Equilibrated flux $\bm{\sigma}_{{\mathcal{T}}}$ and postprocessed potential $\phi_{{\mathcal{T}}}$

The construction of the auxiliary variables $\bm{\sigma}_{{\mathcal{T}}}$ and $\phi_{{\mathcal{T}}}$ giving the equilibration (1.10) is based on independent local mixed finite element approximations of residual problems over the patches of elements around mesh vertices.

For each ${\bm{a}}\in\mathcal{V}$, let ${\mathbb{P}}_{p}({\mathcal{T}_{{\bm{a}}}})$, respectively $\bm{RTN}_{p}({\mathcal{T}_{{\bm{a}}}})$, be the restriction of the space ${\mathbb{P}}_{p}(\mathcal{T})$, respectively $\bm{RTN}_{p}(\mathcal{T})$, to the patch ${\mathcal{T}_{{\bm{a}}}}$ around the vertex ${\bm{a}}$. The local mixed finite element spaces $\bm{V}_{\mathcal{T}}^{\bm{a}}$ and $Q_{\mathcal{T}}^{{\bm{a}}}$ are defined by

[TABLE]

see Figure 1.

Recall that $u_{{\mathcal{T}}}\in V_{{\mathcal{T}}}$ with $V_{{\mathcal{T}}}={\mathbb{P}}_{p}(\mathcal{T})\cap H^{1}_{0}(\Omega)$ is the finite element solution given by (1.5). Let $C_{*}$ be the constant composed of the constants of the trace and inverse inequalities and given by (2.7). Our construction is:

Definition 2.3** (Flux $\bm{\sigma}_{{\mathcal{T}}}$ and potential $\phi_{{\mathcal{T}}}$).**

For each vertex ${\bm{a}}\in\mathcal{V}$, let $(\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}},\phi_{{\mathcal{T}}}^{{\bm{a}}})\in\bm{V}_{\mathcal{T}}^{\bm{a}}\times Q_{\mathcal{T}}^{{\bm{a}}}$ be defined by the local constrained minimization problem

[TABLE]

We remark that for an interior vertex ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$, we have

[TABLE]

by Galerkin orthogonality with $\psi_{{\bm{a}}}\in V_{{\mathcal{T}}}$ as a test function in (1.5). Since

[TABLE]

by Green’s theorem and the vanishing normal flux condition imposed in the definition (2.8a) of $\bm{V}_{\mathcal{T}}^{\bm{a}}$, it follows that $\phi_{{\mathcal{T}}}^{{\bm{a}}}$ necessarily satisfies the mean-value property

[TABLE]

whenever $\kappa>0$. If $\kappa=0$ instead, then $\phi_{{\mathcal{T}}}^{{\bm{a}}}$ is undefined by (2.9a) but one remarks that it is no longer needed anywhere in the paper. In this case, Definition 2.3 coincides with [7, equation (9)], [18, Definition 6.9], or [19, Construction 3.4]; in particular, the Neumann compatibility condition of problem (2.9a) for ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$ follows from (2.10).

In practice, the constrained minimization problem (2.9a) is solved through its Euler–Lagrange equations, which can be reduced to solving a linear system of dimension $\operatorname{dim}\bm{V}_{\mathcal{T}}^{\bm{a}}+\operatorname{dim}Q_{\mathcal{T}}^{{\bm{a}}}$ in the present context. This problem reads: find $(\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}},\phi_{{\mathcal{T}}}^{{\bm{a}}})\in\bm{V}_{\mathcal{T}}^{\bm{a}}\times Q_{\mathcal{T}}^{{\bm{a}}}$ with $\phi_{{\mathcal{T}}}^{{\bm{a}}}=\gamma_{\mathcal{T}}^{{\bm{a}}}+\Pi_{{\mathcal{T}}}(\psi_{{\bm{a}}}u_{{\mathcal{T}}})$ and $(\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}},\gamma_{\mathcal{T}}^{{\bm{a}}})\in\bm{V}_{\mathcal{T}}^{\bm{a}}\times Q_{\mathcal{T}}^{{\bm{a}}}$ such that

[TABLE]

2.4 Properties of $\bm{\sigma}_{{\mathcal{T}}}$ and $\phi_{{\mathcal{T}}}$

We have constructed $\bm{\sigma}_{{\mathcal{T}}}$ and $\phi_{{\mathcal{T}}}$ such that the following holds:

Proposition 2.4** ($\bm{H}(\mathrm{div},\Omega)$-conformity of $\bm{\sigma}_{{\mathcal{T}}}$, equilibration).**

Let $\bm{\sigma}_{{\mathcal{T}}}\in\bm{RTN}_{p}(\mathcal{T})$ and $\phi_{{\mathcal{T}}}\in{\mathbb{P}}_{p}(\mathcal{T})$ be given by Definition 2.3. Then $\bm{\sigma}_{{\mathcal{T}}}$ belongs to $\bm{H}(\mathrm{div},\Omega)$, and $\bm{\sigma}_{{\mathcal{T}}}$ and $\phi_{{\mathcal{T}}}$ satisfy the equilibration property (1.10).

Proof.

First, the $\bm{H}(\mathrm{div},\Omega)$-conformity of $\bm{\sigma}_{{\mathcal{T}}}$ follows from the fact that, for any vertex ${\bm{a}}\in\mathcal{V}$, the zero extension of $\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}}$ belongs to $\bm{H}(\mathrm{div},\Omega)$ as a result of the vanishing normal flux boundary conditions in the space $\bm{V}_{\mathcal{T}}^{\bm{a}}$. Then, to show (1.10), we employ the constraint in (2.9a) together with (2.9c):

[TABLE]

where we have used the fact that the hat functions $\{\psi_{{\bm{a}}}\}_{{\bm{a}}\in\mathcal{V}}$ form a partition of unity over $\Omega$, i.e. $\sum_{{\bm{a}}\in\mathcal{V}}\psi_{{\bm{a}}}=1$. ∎

3 A computable guaranteed a posteriori error estimate

This section presents our guaranteed and fully computable a posteriori error estimate. The following upper bound on the energy norm of the error builds on [9, Theorems 3.1 and 4.4] and [5, Lemma 2]. It employs additionally the concept of a potential reconstruction $\phi_{{\mathcal{T}}}$ that will turn out crucial for a simple and robust flux equilibration. Moreover, it relies on the trace and inverse inequalities of Section 2.2 to make appear the crucial weighs (cut-off factors), with the constant $C_{*}$ given by (2.7).

Theorem 3.1** (Guaranteed a posteriori error estimate).**

Let $u$ be the weak solution of problem (1.1) given by (1.3) and let $u_{{\mathcal{T}}}\in V_{{\mathcal{T}}}$ be its finite element approximation given by (1.5). Let $\bm{\sigma}_{{\mathcal{T}}}\in\bm{RTN}_{p}(\mathcal{T})\cap\bm{H}(\mathrm{div},\Omega)$ and $\phi_{{\mathcal{T}}}\in{\mathbb{P}}_{p}(\mathcal{T})$ be given by Definition 2.3. Then the following upper bound for the energy norm of the error holds:

[TABLE]

where the weights $w_{{K}}$ and $\widetilde{w}_{{K}}$ are respectively defined by

[TABLE]

Proof.

First, we note that the energy norm of the error ${\lvert\kern-1.07639pt\lvert\kern-1.07639pt\lvert u-u_{{\mathcal{T}}}\rvert\kern-1.07639pt\rvert\kern-1.07639pt\rvert}$ is related to the residual $\mathcal{R}(u_{{\mathcal{T}}})\in H^{-1}(\Omega)$, defined by

[TABLE]

through the identity

[TABLE]

cf., e.g., [35, equation (4.1)]. Consider now $\langle\mathcal{R}(u_{{\mathcal{T}}}),v\rangle$ for a fixed function $v\in H^{1}_{0}(\Omega)$. Since $\bm{\sigma}_{{\mathcal{T}}}\in\bm{H}(\mathrm{div},\Omega)$ and $v\in H^{1}_{0}(\Omega)$, Green’s theorem gives $(\bm{\sigma}_{{\mathcal{T}}},\nabla v)+(\nabla{\cdot}\bm{\sigma}_{{\mathcal{T}}},v)=0$, so

[TABLE]

where we have also used the equilibration identity (1.10). We now proceed by estimating each term in (3.4) elementwise.

For each element ${K}\in\mathcal{T}$, we use the identity $(f-\Pi_{{\mathcal{T}}}f,v)_{K}=(f-\Pi_{{\mathcal{T}}}f,v-\Pi_{{\mathcal{T}}}v)_{K}$ and the Poincaré–Friedrichs inequality on the convex element ${K}$, i.e. $\lVert v-\Pi_{{\mathcal{T}}}v\rVert_{K}\leq\frac{h_{K}}{\pi}\lVert\nabla v\rVert_{K}$ for any $v\in H^{1}({K})$, together with the energy error definition (1.13), to obtain the following bound

[TABLE]

Here, actually, a little sharper bound is possible by a convex combination of the two possibilities, but we prefer to use the simple form (3.5) with $\widetilde{w}_{{K}}$ in the form of minimum given by (3.2).

Next, it is clear that

[TABLE]

for each ${K}\in\mathcal{T}$. However, this is not necessarily the sharpest possible estimate in the singularly perturbed regime $\kappa\gg\varepsilon$. Therefore, following the idea of [9, Proof of Theorem 4.4], we use Green’s theorem elementwise together with the fact that $\nabla v=\nabla(v-\overline{v}_{{K}})$, where $\overline{v}_{{K}}$ denotes the mean-value of $v$ on ${K}$. This gives

[TABLE]

The $L^{2}({K})$-stability of the mean-value, $\lVert v-\overline{v}_{{K}}\rVert_{K}\leq\lVert v\rVert_{K}$, Young’s inequality

[TABLE]

and the multiplicative trace inequality (2.3) altogether lead to

[TABLE]

Combined with the inverse inequality (2.4), we find that

[TABLE]

The $L^{2}({K})$-stability of the mean-value, the Poincaré–Friedrichs inequality in the form $\lVert v-\overline{v}_{{K}}\rVert_{K}\leq\frac{h_{K}}{\pi}\lVert\nabla v\rVert_{K}$, and (3.7) yield

[TABLE]

Thus, combined with the inverse inequality (2.4), we find that

[TABLE]

Therefore, combining inequalities (3.6), (3.8), and (3.9), we get

[TABLE]

with $w_{{K}}$ given by (3.2) and $C_{*}$ given in (2.7). As a side remark, it is possible to obtain a slightly sharper, at the expense of making the weight $w_{{K}}$ more complicated than the simple form given by (3.2).

Finally, we can apply the Cauchy–Schwarz inequality to see that $|(\kappa(\phi_{{\mathcal{T}}}-u_{{\mathcal{T}}}),\kappa v)_{{K}}|\leq\lVert\kappa\left(u_{{\mathcal{T}}}-\phi_{{\mathcal{T}}}\right)\rVert_{{K}}{\lvert\kern-1.07639pt\lvert\kern-1.07639pt\lvert v\rvert\kern-1.07639pt\rvert\kern-1.07639pt\rvert}_{{K}}$. Therefore, we deduce from (3.4) and the above inequalities that

[TABLE]

which implies the upper bound on the error (3.1) after another Cauchy–Schwarz inequality, using (3.3) and $\sum_{{K}\in\mathcal{T}}{\lvert\kern-1.07639pt\lvert\kern-1.07639pt\lvert v\rvert\kern-1.07639pt\rvert\kern-1.07639pt\rvert}_{K}^{2}={\lvert\kern-1.07639pt\lvert\kern-1.07639pt\lvert v\rvert\kern-1.07639pt\rvert\kern-1.07639pt\rvert}^{2}$. ∎

4 Efficiency and robustness of the estimate

This section establishes the local (and consequently global) efficiency and robustness of our a posteriori error estimate.

4.1 A basic stability result

The main tool in the analysis of efficiency is the following stability result, where, we recall, the broken and the patchwise ${\bm{H}}(\mathrm{div})$-conforming Raviart–Thomas–Nédélec spaces $\bm{RTN}_{p}({\mathcal{T}_{{\bm{a}}}})$ and $\bm{V}_{\mathcal{T}}^{\bm{a}}$ are respectively given by (2.2) and (2.8).

Lemma 4.1** (Stability of patchwise flux equilibration).**

Let a vertex ${\bm{a}}\in\mathcal{V}$ be fixed, and let $g_{{\mathcal{T}}}\in{\mathbb{P}}_{p}({\mathcal{T}_{{\bm{a}}}})$ and $\bm{\tau}_{{\mathcal{T}}}\in\bm{RTN}_{p}({\mathcal{T}_{{\bm{a}}}})$ be given discontinuous piecewise polynomial functions, with the Neumann compatibility condition $(g_{{\mathcal{T}}},1)_{{\omega_{\bm{a}}}}=0$ satisfied if ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$. Then, there holds

[TABLE]

where $H^{1}_{*}({\omega_{\bm{a}}})$ is the subspace of functions in $H^{1}({\omega_{\bm{a}}})$ that have mean-value zero on the patch subdomain ${\omega_{\bm{a}}}$ if ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$ is an interior vertex, or that vanish on $\partial{\omega_{\bm{a}}}\cap\partial\Omega$ if ${\bm{a}}\in\mathcal{V}^{\mathrm{ext}}$ is a boundary vertex.

The above result holds for any dimension $d\geq 1$, although some additional properties are known for $d\leq 3$. Indeed, in the case where $d=2$, it is shown in [7, Theorem 7] that the constant in (4.1) is in fact independent of the polynomial degree $p$, i.e. $p$-robust. The extension of the $p$-robustness of the bound to the case of $d=3$ was shown in [20, Corollaries 3.3 and 3.6]. It is also possible to extend similar results of this kind to situations with hanging nodes and locally refined submeshes, as shown in [17].

4.2 Stability with respect to residual estimators

The next lemma shows that the local contributions of the equilibrated flux a posteriori estimators of Definition 2.3 lie below the local residual estimators as defined in (1.7), with the element residuals $r_{\mathcal{T}}$ and face residuals $j_{\mathcal{T}}$ are defined by (1.8) and the weights $\alpha_{{K}}$ and $\alpha_{F}$ defined by (1.9).

Lemma 4.2** (Stability of patchwise flux equilibration with respect to residual estimators).**

For each ${\bm{a}}\in\mathcal{V}$, let $\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}}$ and $\phi_{{\mathcal{T}}}^{{\bm{a}}}$ be defined by (2.9a). Then

[TABLE]

Proof.

Let a vertex ${\bm{a}}\in\mathcal{V}$ be fixed. Since $\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}}$ and $\phi_{{\mathcal{T}}}^{{\bm{a}}}$ are defined as minimizers of the functional in the right-hand side of (2.9a), it is enough to prove that there always exists some $\bm{v}_{{\mathcal{T}}}^{*}\in\bm{V}_{\mathcal{T}}^{\bm{a}}$ and $q_{\mathcal{T}}^{*}\in Q_{\mathcal{T}}^{{\bm{a}}}$ that satisfy the constraint $\nabla{\cdot}\bm{v}_{{\mathcal{T}}}^{*}+\kappa^{2}q_{\mathcal{T}}^{*}=\Pi_{{\mathcal{T}}}(f\psi_{{\bm{a}}})-\varepsilon^{2}\nabla u_{{\mathcal{T}}}{\cdot}\nabla\psi_{{\bm{a}}}$ and that satisfy the bound (4.2) with $\bm{v}_{{\mathcal{T}}}^{*}$ in place of $\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}}$ and $q_{\mathcal{T}}^{*}$ in place of $\phi_{{\mathcal{T}}}^{{\bm{a}}}$. The specific construction depends on the mesh size and the problem parameters $\varepsilon$ and $\kappa$, as we now show.

Case 1, $\varepsilon/h_{{\omega_{\bm{a}}}}\leq\kappa$ (reaction dominance). Up to a constant, we have $\kappa^{-1}\lesssim h_{{K}}/\varepsilon$ and $\kappa^{-1}\lesssim h_{F}/\varepsilon$ for all elements ${K}\in{\mathcal{T}_{{\bm{a}}}}$ and all interior faces $F\in\mathcal{F}_{{\bm{a}}}$. In this case, we adopt the following construction. Let

[TABLE]

and $\rho_{{\bm{a}}}\coloneqq 0$ otherwise. Next, we define

[TABLE]

where

[TABLE]

It is easy to check that if ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$, then $(g_{{\mathcal{T}}}^{*},1)_{{\omega_{\bm{a}}}}=0$, since the Galerkin orthogonality (take $v_{{\mathcal{T}}}=\psi_{{\bm{a}}}$ in (1.5)) implies that

[TABLE]

Therefore, it follows that $q_{\mathcal{T}}^{*}\in Q_{\mathcal{T}}^{{\bm{a}}}$ and $\bm{v}_{{\mathcal{T}}}^{*}\in\bm{V}_{\mathcal{T}}^{\bm{a}}$ are well-defined and that they satisfy the constraint $\nabla{\cdot}\bm{v}_{{\mathcal{T}}}^{*}+\kappa^{2}q_{\mathcal{T}}^{*}=\Pi_{{\mathcal{T}}}(f\psi_{{\bm{a}}})-\varepsilon^{2}\nabla u_{{\mathcal{T}}}{\cdot}\nabla\psi_{{\bm{a}}}$.

We now bound $w_{{\bm{a}}}^{2}\lVert\varepsilon^{2}\psi_{{\bm{a}}}\nabla u_{{\mathcal{T}}}+\bm{v}_{{\mathcal{T}}}^{*}\rVert_{{\omega_{\bm{a}}}}^{2}$ and $\lVert\kappa\left[\Pi_{{\mathcal{T}}}(\psi_{{\bm{a}}}u_{{\mathcal{T}}})-q_{\mathcal{T}}^{*}\right]\rVert_{{\omega_{\bm{a}}}}^{2}$. First, we obtain

[TABLE]

where we have used the stability of the $L^{2}$-projection (note that $\rho_{{\bm{a}}}$ is also the mean value of $\Pi_{{\mathcal{T}}}(\psi_{{\bm{a}}}r_{\mathcal{T}})$ on ${\omega_{\bm{a}}}$ for ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$) and the fact that $\lVert\psi_{{\bm{a}}}\rVert_{\infty,{\omega_{\bm{a}}}}=1$ to bound $\lVert\Pi_{{\mathcal{T}}}(\psi_{{\bm{a}}}r_{\mathcal{T}})-\rho_{{\bm{a}}}\rVert_{{\omega_{\bm{a}}}}$. Next, we apply Lemma 4.1 to bound $w_{{\bm{a}}}^{2}\lVert\varepsilon^{2}\psi_{{\bm{a}}}\nabla u_{{\mathcal{T}}}+\bm{v}_{{\mathcal{T}}}^{*}\rVert_{{\omega_{\bm{a}}}}^{2}$. Note first that for an interior vertex ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$, $(\rho_{{\bm{a}}},v)_{{\omega_{\bm{a}}}}=0$ since $v\in H^{1}_{*}({\omega_{\bm{a}}})$ implies that $v$ is orthogonal to constant functions on ${\omega_{\bm{a}}}$. We find that

[TABLE]

where the last line follows by elementwise integration by parts. It is then straightforward to deduce from the trace inequalities $\lVert v\rVert_{F}\lesssim h_{K}^{-\frac{1}{2}}\lVert v\rVert_{{K}}+\lVert\nabla v\rVert_{{K}}^{\frac{1}{2}}\lVert v\rVert_{{K}}^{\frac{1}{2}}$ and the Poincaré–Friedrichs inequality for functions in $H^{1}_{*}({\omega_{\bm{a}}})$ $\lVert v\rVert_{\omega_{\bm{a}}}\lesssim h_{\omega_{\bm{a}}}\lVert\nabla v\rVert_{\omega_{\bm{a}}}$ that

[TABLE]

Consequently, using definition (2.9b) of the weight $w_{{\bm{a}}}$

[TABLE]

Therefore, if $\varepsilon/h_{{\omega_{\bm{a}}}}\leq\kappa$, we have shown that there exists $\bm{v}_{{\mathcal{T}}}^{*}$ and $q_{\mathcal{T}}^{*}$ satisfying the constraint $\nabla{\cdot}\bm{v}_{{\mathcal{T}}}^{*}+\kappa^{2}q_{\mathcal{T}}^{*}=\Pi_{{\mathcal{T}}}(f\psi_{{\bm{a}}})-\varepsilon^{2}\nabla u_{{\mathcal{T}}}{\cdot}\nabla\psi_{{\bm{a}}}$ and such that

[TABLE]

As explained above, this implies (4.2) in the case $\varepsilon/h_{{\omega_{\bm{a}}}}\leq\kappa$.

Case 2, $\varepsilon/h_{{\omega_{\bm{a}}}}>\kappa$ (diffusion dominance). We select

[TABLE]

where

[TABLE]

Notice that Galerkin orthogonality implies that $(g_{{\mathcal{T}}}^{*},1)_{{\omega_{\bm{a}}}}=0$ if ${\bm{a}}\in\mathcal{V}^{\mathrm{int}}$ as in (4.4), and also $\nabla{\cdot}\bm{v}_{{\mathcal{T}}}^{*}+\kappa^{2}q_{\mathcal{T}}^{*}=\Pi_{{\mathcal{T}}}(f\psi_{{\bm{a}}})-\varepsilon^{2}\nabla u_{{\mathcal{T}}}{\cdot}\nabla\psi_{{\bm{a}}}$, so the requested constraint is satisfied. It then follows directly from Lemma 4.1 that

[TABLE]

where we use the fact that elementwise integration by parts shows that, as in (4.5),

[TABLE]

Thus, proceeding as in (4.5)–(4.6) for the face residuals term and using the stability of the $L^{2}$-projection, $\lVert\psi_{{\bm{a}}}\rVert_{\infty,{\omega_{\bm{a}}}}=1$, and the Poincaré–Friedrichs inequality for functions in $H^{1}_{*}({\omega_{\bm{a}}})$, $\lVert v\rVert_{\omega_{\bm{a}}}\lesssim h_{\omega_{\bm{a}}}\lVert\nabla v\rVert_{\omega_{\bm{a}}}$, for the element residuals term, we get

[TABLE]

Consequently,

[TABLE]

Hence, on noting that $w_{{\bm{a}}}\leq 1$ and that $\lVert\kappa\left[\Pi_{{\mathcal{T}}}(\psi_{{\bm{a}}}u_{{\mathcal{T}}})-q_{\mathcal{T}}^{*}\right])\rVert_{{\omega_{\bm{a}}}}=0$, we see that (4.2) also holds for the case $\varepsilon/h_{{\omega_{\bm{a}}}}>\kappa$. ∎

Recall that $\mathfrak{T}_{K}\coloneqq\bigcup_{{\bm{a}}\in\mathcal{V}_{K}}{\mathcal{T}_{{\bm{a}}}}$ and $\mathfrak{F}_{K}\coloneqq\bigcup_{{\bm{a}}\in\mathcal{V}_{K}}\mathcal{F}_{{\bm{a}}}$.

Proposition 4.3** (Bound on flux estimators by the residual estimators).**

Let $\bm{\sigma}_{{\mathcal{T}}}$ and $\phi_{{\mathcal{T}}}$ be given by Definition 2.3. Additionally, let the volume and face residual functions $r_{\mathcal{T}}$ and $j_{\mathcal{T}}$ be defined by (1.8). Then, for each element ${K}\in\mathcal{T}$, we have the bound

[TABLE]

Proof.

For each mesh element ${K}\in\mathcal{T}$, we have $\bm{\sigma}_{{\mathcal{T}}}|_{{K}}=\sum_{{\bm{a}}\in\mathcal{V}_{K}}\bm{\sigma}_{{\mathcal{T}}}^{{\bm{a}}}|_{{K}}$ and $\phi_{{\mathcal{T}}}|_{{K}}=\sum_{{\bm{a}}\in\mathcal{V}_{K}}\phi_{{\mathcal{T}}}^{{\bm{a}}}|_{{K}}$. Furthermore, since $\{\psi_{{\bm{a}}}\}_{{\bm{a}}\in\mathcal{V}_{K}}$ form a partition of unity over ${K}$ and since $\Pi_{{\mathcal{T}}}$ is the elementwise $L^{2}$ projection of degree $p$, it follows that $u_{{\mathcal{T}}}|_{{K}}=\Pi_{{\mathcal{T}}}u_{{\mathcal{T}}}|_{{K}}=\sum_{{\bm{a}}\in\mathcal{V}_{K}}\Pi_{{\mathcal{T}}}\left(\psi_{{\bm{a}}}u_{{\mathcal{T}}}\right)|_{{K}}$. Furthermore, (3.2) and (2.9b) together with the mesh shape regularity imply that $w_{{K}}\lesssim w_{{\bm{a}}}$ for each ${\bm{a}}\in\mathcal{V}_{K}$, where the constant depends only on $\vartheta_{\mathcal{T}}$. Therefore, we obtain

[TABLE]

Therefore, we can use (4.2) for each ${\bm{a}}\in\mathcal{V}_{K}$ to get (4.7). ∎

4.3 Local efficiency and robustness of the estimate

We now recall the well-known efficiency and robustness results for residual estimators, see [35, Proposition 4.1] and [38] for details. For each ${K}\in\mathcal{T}$ and $F\in\mathcal{F}_{\Omega}$, there holds

[TABLE]

Therefore, the combination of Proposition 4.3 with (4.8) shows that the equilibrated flux estimator of Theorem 3.1 is locally efficient and robust.

Theorem 4.4** (Local efficiency and robustness).**

Let $u$ be the weak solution of problem (1.1) given by (1.3) and let $u_{{\mathcal{T}}}\in V_{{\mathcal{T}}}$ be its finite element approximation given by (1.5). Let $\bm{\sigma}_{{\mathcal{T}}}\in\bm{RTN}_{p}(\mathcal{T})\cap\bm{H}(\mathrm{div},\Omega)$ and $\phi_{{\mathcal{T}}}\in{\mathbb{P}}_{p}(\mathcal{T})$ be given by Definition 2.3. Then, for each mesh element ${K}\in\mathcal{T}$, there holds

[TABLE]

where the constant in $\lesssim$ depends only on the dimension $d$, the shape-regularity constant $\vartheta_{\mathcal{T}}$ of $\mathcal{T}$, and on the polynomial degree $p$, so that it is independent of the parameters $\varepsilon$ and $\kappa$ and the mesh-sizes $h_{K}$.

5 Necessity of the weights $w_{{K}}$ in the upper bound

Theorems 3.1 and 4.4 show that the estimator $w_{{K}}\lVert\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}\bm{\sigma}_{{\mathcal{T}}}\rVert_{{K}}+\lVert\kappa\left(u_{{\mathcal{T}}}-\phi_{{\mathcal{T}}}\right)\rVert_{{K}}$ obtained from the flux equilibration of Definition 2.3 is a reliable, locally efficient, and robust energy error estimator for singularly perturbed reaction–diffusion problems. Here we show the necessity of the weight $w_{{K}}$ for robustness of equilibrated flux estimators that involve only piecewise polynomial vector fields on $\mathcal{T}$. We also recall that an alternative option, related to the approach in [2, 4, 5, 25], is to perform an equilibrations on a submesh.

5.1 Necessity of the weights $w_{{K}}$

The following proposition applies to any flux equilibration on $\mathcal{T}$:

Proposition 5.1** (Best-possible bound by piecewise polynomials of the mesh $\mathcal{T}$).**

Let $u_{{\mathcal{T}}}\in{\mathbb{P}}_{p}(\mathcal{T})\cap H^{1}_{0}(\Omega)$ be an arbitrary piecewise $p$-degree polynomial, $p\geq 1$, and let the face residual term $j_{\mathcal{T}}$ be defined by (1.8b). Let ${\mathbb{P}}_{p^{\prime}}(\mathcal{T};\mathbb{R}^{d})$ denote the space of $\mathbb{R}^{d}$-valued piecewise polynomials of degree at most $p^{\prime}$ over $\mathcal{T}$, where $p^{\prime}\geq 0$ is an arbitrary nonnegative integer. Then,

[TABLE]

where $\underline{h}\coloneqq\min_{{K}\in\mathcal{T}}h_{K}$, and where the constant depends only on the polynomial degrees $p$ and $p^{\prime}$, the dimension $d$, and the shape-regularity $\vartheta_{\mathcal{T}}$ of $\mathcal{T}$.

Proof.

Let $\bm{v}_{{\mathcal{T}}}\in\bm{H}(\mathrm{div},\Omega)\cap{\mathbb{P}}_{p^{\prime}}(\mathcal{T};\mathbb{R}^{d})$ be arbitrary. Then, for each interior face $F\in\mathcal{F}_{\Omega}$, the $\bm{H}(\mathrm{div},\Omega)$-conformity of $\bm{v}_{{\mathcal{T}}}$ implies that $\llbracket\bm{v}_{{\mathcal{T}}}{\cdot}\bm{n}_{F}\rrbracket_{F}=0$, and hence $j_{\mathcal{T}}|_{F}=-\varepsilon\llbracket(\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}\bm{v}_{{\mathcal{T}}}){\cdot}\bm{n}_{F}\rrbracket_{F}$. Since $(\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}\bm{v}_{{\mathcal{T}}})|_{K}\in{\mathbb{P}}_{\max(p^{\prime},p-1)}({K};\mathbb{R}^{d})$ for each element $K\in\mathcal{T}$, we can apply the triangle inequality and the inverse inequality (analogous to (2.4)) to find that, for any $F\in\mathcal{F}_{\Omega}$,

[TABLE]

Therefore, we get (5.1) by summing (5.2) over all faces $F\in\mathcal{F}_{\Omega}$, and recalling that $\bm{v}_{{\mathcal{T}}}$ was arbitrary. ∎

The upshot of Proposition 5.1 is that for any problem where the jump estimators are sufficiently dominant, i.e. when

[TABLE]

then any error estimator involving a term of the form $\lVert\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}\bm{v}_{{\mathcal{T}}}\rVert$ without any weight will necessarily be non-robust when $\kappa\underline{h}/\varepsilon$ takes large values, since (5.1) and (5.3) then imply

[TABLE]

In other words, the effectivity index can be become arbitrarily large in the singularly-perturbed regime when the weight $w_{{K}}$ is not included. It is then seen that the inclusion of the weight term $w_{{K}}$ in Theorem 3.1 is necessary when considering flux equilibrations from vector-valued piecewise polynomial subspaces of $\bm{H}(\mathrm{div},\Omega)$ on the mesh $\mathcal{T}$, regardless of the precise details of the construction of the flux. Examples of flux equilibrations proposed in the past that cannot be robust in general include Repin and Sauter [29], Ainsworth et al. [1], Eigel and Samrowski [16], Eigel and Merdon [15], and Vejchodský [32, 33, 34].

We now present an example of a situation where (5.3) holds and where $\kappa\underline{h}/\varepsilon$ can be arbitrarily large. In fact the example is similar to the one in [2, Section 2.3], albeit with some suitable adjustments.

Example 5.2** (Dominant jump estimators).**

Let $\Omega\coloneqq(-1/2,1/2)$ and let $m$ be an odd integer that will later on be chosen sufficiently large. Consider a uniform mesh $\mathcal{T}$ of $\Omega$ with $2N=(m+1)^{2}$ intervals, $N\coloneqq(m+1)^{2}/2$, and mesh size $h\coloneqq 1/(2N)=1/(m+1)^{2}$. Hence, the interior nodes are $x_{i}=ih$, where $i\in\{-N+1,\dots,N-1\}$. Let

[TABLE]

denote the piecewise affine Lagrange interpolant (preserving the point values) of the function $x\mapsto\cos(m\pi x)$; it follows from the fact that $m$ is odd that $f\in H^{1}_{0}(\Omega)$. Note that in the example of [2], the function $f$ was chosen as $\cos(\pi x)$ instead.

Consider now problem (1.1) along with its finite element approximation (1.5) in the space $V_{{\mathcal{T}}}={\mathbb{P}}_{1}(\mathcal{T})\cap H^{1}_{0}(\Omega)$. It is easy to show that

[TABLE]

is the discrete solution, where

[TABLE]

as a result of the identity

[TABLE]

Then, noting that interior vertices and faces coincide for problems in one space dimension, it is found that

[TABLE]

Now, since $\lim_{m\rightarrow\infty}\frac{1-\cos(m\pi h)}{h}=\frac{\pi^{2}}{2}$ when $h=h(m)=1/(m+1)^{2}$, we can pick $m$ sufficiently large such that $\mu_{h}\simeq h^{-1}$. Suppose also henceforth that $\kappa h/\varepsilon\geq 1$, so that $\alpha_{K}$ given by (1.9) takes the value $1/\kappa$. Then, we find that

[TABLE]

We also obtain

[TABLE]

where we have used the trigonometric identity

[TABLE]

Since $\varepsilon\kappa h\leq\kappa^{2}h^{2}$, we see that

[TABLE]

where we note that there is no data oscillation since $f\in{\mathbb{P}}_{1}(\mathcal{T})$. Hence this provides an example where (5.3) holds, and the factor $\kappa h/\varepsilon$ can be made arbitrarily large.

5.2 Flux equilibration on a submesh

We finish with the following remark:

Remark 5.3** (Flux equilibration on boundary-layer adapted submeshes).**

The approach in [4, 5, 25], following [2], can be seen as defining a flux $\bm{\sigma}_{\widetilde{\mathcal{T}}}\in\bm{H}(\mathrm{div},\Omega)$ that satisfies an equilibration property similar to (1.10), yet with the key difference that $\bm{\sigma}_{\widetilde{\mathcal{T}}}$ is defined with respect to a submesh $\widetilde{\mathcal{T}}$ of $\mathcal{T}$ with thin elements that are adapted to the parameters $\varepsilon$ and $\kappa$ and local mesh-size (see e.g. [4, Fig. 3]). In this case, the argument in the proof of Proposition 5.1 does not apply, because the inverse inequality $\lVert(\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}\bm{\sigma}_{\widetilde{\mathcal{T}}}){\cdot}{\bm{n}}_{F}\rVert_{F}\lesssim h_{K}^{-\frac{1}{2}}\lVert\varepsilon\nabla u_{{\mathcal{T}}}+\varepsilon^{-1}\bm{\sigma}_{\widetilde{\mathcal{T}}}\rVert_{{K}}$, $F\subset\partial{K}$, ${K}\in\mathcal{T}$, is not applicable when $\bm{\sigma}_{\widetilde{\mathcal{T}}}\in{\mathbb{P}}_{p^{\prime}}(\widetilde{\mathcal{T}};\mathbb{R}^{d})$ but $\bm{\sigma}_{\widetilde{\mathcal{T}}}\notin{\mathbb{P}}_{p^{\prime}}(\mathcal{T};\mathbb{R}^{d})$. This essentially shows how there are now two different approaches to constructing robust equilibrated flux estimators. Either the flux is computed as a piecewise polynomial vector field with respect to the original mesh, in which case the inclusion of a weight of the form of $w_{{K}}$ from (3.2) in the upper bound is necessary, or one constructs the flux with respect to some other sufficiently rich subspace of $\bm{H}(\mathrm{div},\Omega)$, such as a piecewise polynomial subspace with respect to an adapted submesh $\widetilde{\mathcal{T}}$ of $\mathcal{T}$, in which case the weights are not necessary.

Appendix A Explicit constants for the inverse inequality

For each polynomial degree $p\geq 0$, let $C_{p,1}$ denote the best constant of the inverse inequality for the unit interval $(0,1)$, i.e.

[TABLE]

where ${\mathbb{P}}_{p}(0,1)$ denotes the space of univariate polynomials of degree at most $p$ on $(0,1)$. It was shown in [26] that, for all $p\geq 0$,

[TABLE]

where we have taken into account the fact that we consider $C_{p,1}$ on the unit interval $(0,1)$ rather than the interval $(-1,1)$ as in [26]. This improves on earlier bounds, e.g. in [30].

We will show here explicit bounds for the constants of the inverse inequality for hypercubes and simplices in terms of $C_{p,1}$.

A.1 Unit hypercube

For an integer $d\geq 1$, let $\{1{:}d\}$ be a shorthand notation for $\{1,\dots,d\}$. Let $Q_{d}\coloneqq\{x\in\mathbb{R}^{d},\;\lvert x\rvert_{\infty}\leq 1,\;x_{i}\geq 0\;\;\forall i\in\{1{:}d\}\}$ denote the unit hypercube in $\mathbb{R}^{d}$, where $\lvert x\rvert_{\infty}\coloneqq\max_{i\in\{1{:}d\}}\lvert x_{i}\rvert$. Let ${\mathbb{P}}_{p}(Q_{d})$ denote the space of polynomials of total degree at most $p$ on $Q_{d}$.

Lemma A.1**.**

For all $d\geq 1$ and all $p\geq 0$, we have

[TABLE]

Proof.

After a possible re-labelling of the indices, it is enough to show that (A.3) holds for the case $i=1$. Then, writing $x=(x_{1},x^{\prime})$ with $x^{\prime}\in\mathbb{R}^{d-1}$, we see that

[TABLE]

where we use the fact that $x_{1}\mapsto v(x_{1},x^{\prime})$ is in ${\mathbb{P}}_{p}(0,1)$ for all $x^{\prime}$. ∎

A.2 Unit simplex

For a parameter $t>0$, let $K^{d}_{t}\coloneqq\{x\in\mathbb{R}^{d},\;\lvert x\rvert_{1}\leq t,\;x_{i}\geq 0\;\;\forall i\in\{1{:}d\}\}$, where $\lvert x\rvert_{1}\coloneqq\sum_{i=1}^{d}\lvert x_{i}\rvert$, denote the simplex in $\mathbb{R}^{d}$ with side-length $t$. If $t=1$, we adopt the simpler notation $K^{d}\coloneqq K^{d}_{1}$. Let $C_{p,d}$ denote the best constant such that

[TABLE]

We shall obtain here an explicit bound for the constant $C_{p,d}$ in terms of the space dimension $d$ and the constant $C_{p,1}$ of (A.1).

Theorem A.2**.**

For all $d\geq 1$ and for all $p\geq 0$, the best constant $C_{p,d}$ in (A.4) satisfies

[TABLE]

Proof.

The proof is based on an induction on the dimension, where we seek to bound $C_{p,d}$ in terms of $C_{p,d-1}$, $C_{p,1}$, and $d$. Without loss of generality, it is enough to consider only the case $i=1$ in (A.4), after a possible re-labelling of the indices. Then, writing $x=(x^{\prime},x_{d})$ with $x^{\prime}\in\mathbb{R}^{d-1}$, we have $\int_{K^{d}}\lvert v_{x_{1}}\rvert^{2}\mathrm{d}x=\int_{0}^{1}\int_{K^{d-1}_{1-x_{d}}}\lvert v_{x_{1}}\rvert^{2}\mathrm{d}x^{\prime}\mathrm{d}x_{d}$. Since, for fixed $x_{d}\in(0,1)$, $x^{\prime}\mapsto v(x^{\prime},x_{d})$ is a polynomial of degree at most $p$ on $K^{d-1}_{1-x_{d}}$, it would be natural to apply the inverse inequality for simplices of dimension $d-1$ after a suitable scaling. However, a difficulty arises for $x_{d}$ close to $1$ due to the appearance of a negative power of $1-x_{d}$ inside the resulting integral. We can overcome this obstacle using an appropriate subdivision of the unit simplex and a change of variables.

The proof proceeds in two steps. We first treat the case $d=2$ and show that (A.5) holds (we actually consider $d\geq 2$ below for the sake of generality), and then the induction is carried out on $d$ with a different argument, leading to a sharper bound than that would result from step 1 only.

Step 1. Let $d\geq 2$ and consider the partition of $K$ into $K_{*}\coloneqq\{x\in K,x_{d}<1-1/d\}$ and $K_{\dagger}\coloneqq K\setminus K_{*}$. Then, $\lVert v_{x_{1}}\rVert_{L^{2}(K^{d})}^{2}=\lVert v_{x_{1}}\rVert_{L^{2}(K_{*})}^{2}+\lVert v_{x_{1}}\rVert_{L^{2}(K_{\dagger})}^{2}$, and the first term can be bounded as follows:

[TABLE]

where crucially we use the fact that $(1-x_{d})^{-2}\leq d^{2}$ for $x_{d}\leq 1-1/d$. In order to bound the second term $\lVert v_{x_{1}}\rVert^{2}_{L^{2}(K_{\dagger})}$, we introduce a change of coordinates in terms of the affine map $F$ defined by

[TABLE]

where $e_{0}=0$, and $e_{i}$ is the $i$-th unit vector for $1\leq i\leq d$. Letting $x=F(\xi)$, we have $x_{j}=\xi_{j+1}$ for $j\leq d-1$, and $x_{d}=1-\sum_{i=1}^{d}\xi_{i}$. The inverse is then given by $\xi_{1}=1-\sum_{i=1}^{d}x_{i}$, and $\xi_{j}=x_{j-1}$ for $2\leq j\leq d$. It is thus easily seen that $F$ is a bijection from $K$ onto itself, and that $F(0)=x_{d}$. Thus $F$ corresponds to a change of coordinates on the unit simplex. Additionally, it can be shown that the Jacobian $\lvert\mathrm{det}DF\rvert=1$.

Let $Q^{d}_{1/d}\coloneqq\{\xi\in Q^{d},\,\lvert\xi\rvert_{\infty}\leq 1/d\}$ be a hypercube with side length $1/d$, and let $Q_{\dagger}$ be the parallelepiped obtained as the image of $Q^{d}_{1/d}$ under the mapping $F$, i.e. $Q_{\dagger}=F(Q^{d}_{1/d})$. It is then easy, but tedious, to show that

[TABLE]

Figure 2 illustrates the sets $K_{\dagger}$, $Q_{\dagger}$, and $K$ for the case $d=3$. Now, let $\widetilde{v}(\xi)=v(F(\xi))$ be the pullback of $v$ under $F$. Since $F$ is affine, $\widetilde{v}\in{\mathbb{P}}_{p}(K^{d})$. It is also easy to check that $v_{x_{1}}=\widetilde{v}_{\xi_{2}}-\widetilde{v}_{\xi_{1}}$. Using the change of variables and the fact that $\lvert\mathrm{det}DF\rvert=1$, it follows from (A.7) that

[TABLE]

Applying the inverse inequality for hypercubes, namely $\lVert\widetilde{v}_{\xi_{i}}\rVert_{L^{2}(Q^{d}_{1/d})}^{2}\leq d^{2}C_{p,1}^{2}\lVert\widetilde{v}\rVert_{L^{2}(Q^{d}_{1/d})}^{2}$, and changing back to the original variables, we then obtain from the second inclusion in (A.7) that

[TABLE]

Therefore, combining (A.6) and (A.8), we arrive at $\lVert v_{x_{1}}\rVert^{2}_{L^{2}(K^{d})}\leq d^{2}(C_{p,d-1}^{2}+4C_{p,1}^{2})\lVert v\rVert_{L^{2}(K^{d})}^{2}$, for any $v\in{\mathbb{P}}_{p}(K^{d})$. This implies $C_{p,d}^{2}\leq d^{2}(C_{p,d-1}^{2}+4C_{p,1}^{2})$, and thus (A.1) and an induction argument show that

[TABLE]

for any $d\geq 2$. This shows (A.4), but with a worse constant than that of (A.5) for $d\geq 3$. For this reason, we proceed in a second step in a different way.

Step 2. Let $d\geq 3$. We again subdivide the simplex $K$, this time as

[TABLE]

Furthermore, for any fixed $x_{d-1}$, $x_{d-1}^{\prime}=(x_{1},\dots,x_{d-2},x_{d})^{\prime}\mapsto v(x)$ is a polynomial of degree at most $p$ on a simplex that is isometric to $K^{d-1}_{1-x_{d-1}}$. Let also $x_{d}^{\prime}=(x_{1},\dots,x_{d-1})$. Crucially, since $d\geq 3$ and we subdivide above into two subsets, we can avoid the critical subset $K_{\dagger}$ of Step 1 as

[TABLE]

It then follows by induction that $C_{p,d}\leq(2\sqrt{2})^{d-2}C_{p,2}$ for all $d\geq 3$. Since $C_{p,2}\leq 2\sqrt{5}\,C_{p,1}$ by (A.9), we get (A.5). ∎

Applying Theorem A.2 to the cases $d=2$ and $d=3$ gives the following explicit bounds

[TABLE]

A.3 General simplex

Let $K$ be a simplex in $\mathbb{R}^{d}$, $d\geq 2$, and let $\hat{K}$ denote the unit simplex. Let $J_{K}$ denote the differential of the affine transformation mapping $T_{K}\colon\hat{K}\rightarrow K$. For $\bm{v}\in\bm{RTN}_{p}(K)$, we define the Piola transformation $\bm{\hat{v}}\in\bm{RTN}_{p}(\hat{K})$ by

[TABLE]

Lemma A.3** (Ciarlet [10] Thm 3.1.2 and [13]).**

There holds

[TABLE]

Note that in Lemma A.3, we have used the fact that the diameter of the unit simplex is $\sqrt{2}$ for all $d\geq 2$.

Lemma A.4** (Hesthaven & Warburton [39]).**

Let $v\in{\mathbb{P}}_{p}(K)$. Then

[TABLE]

Therefore, for $\bm{v}\in\bm{RTN}_{p}(K)$, we have

[TABLE]

Lemma A.5**.**

Let $K$ be a simplex in $\mathbb{R}^{d}$ and $\bm{v}\in\bm{RTN}_{p}(K)$. Then,

[TABLE]

where $C_{p,d}$ is characterized in Theorem A.2.

Proof.

Using the Piola Transformation, we have $\nabla\cdot\bm{v}=\nabla_{\hat{x}}\cdot\bm{\hat{v}}/\lvert\det J_{K}\rvert$, therefore, $\lVert\nabla\cdot\bm{v}\rVert_{K}^{2}\leq\lvert\det J_{K}\rvert^{-1}\lVert\nabla\cdot\bm{\hat{v}}\rVert_{\hat{K}}^{2}$. Then, since $\bm{\hat{v}}_{i}\in{\mathbb{P}}_{p+1}(\hat{K})$ for $i\in\{1{:}d\}$, we apply Theorem A.2 to obtain $\lVert\nabla\cdot\bm{\hat{v}}\rVert_{\hat{K}}^{2}\leq d\sum_{i=1}^{d}\lVert\bm{\hat{v}}_{i,x_{i}}\rVert_{\hat{K}}^{2}\leq dC_{p+1,d}^{2}\lVert\bm{\hat{v}}\rVert_{\hat{K}}^{2}$. Then, using the definition of the Piola transformation, it is seen that $\lVert\bm{\hat{v}}\rVert_{\hat{K}}^{2}\leq\lVert J_{K}^{-1}\rVert_{2}^{2}\lvert\det J_{K}\rvert\lVert\bm{v}\rVert_{K}^{2}$. We then use the bound $\lVert J_{K}^{-1}\rVert_{2}\leq\sqrt{2}\vartheta_{\mathcal{T}}h_{K}^{-1}$ from (A.13) to find that $h_{K}\lVert\nabla\cdot\bm{v}\rVert_{K}\leq\sqrt{2d}\vartheta_{\mathcal{T}}C_{p+1,d}\lVert\bm{v}\rVert_{K}$, which finishes the proof. ∎