Stabilization-free HHO a posteriori error control

Hybrid high-order methods (HHO) were introduced in [26, 27] and are examined in the textbooks [25, 29] as a promising class of flexible nonconforming discretization methods for partial differential equations that involve a parameter-free stabilization term for the link between the volume and skeletal variables.

1.1 Known a posteriori error estimator

The a priori error analysis of HHO involves the stability terms in extended norms as part of the methodology and motivated a first explicit residual-based a posteriori error estimator in [25] with a reformulation of the stabilization in the upper bound. Let \(s_h(u_h,u_h)\) denote the stabilization at the discrete solution \(u_h\in V_h\) and let the (elliptic) reconstruction \(Ru_h\) of \(u_h\) denote a piecewise polynomial of degree at most \(k+1\) that approximates \(u\in H^1(\Omega )\), cf. (2) and Sect. 3 below for further details. Then a possible error term reads

It is disputable if \(s_h(u_h,u_h)\ge 0\) is an error contribution, but if the total error includes \(s_h(u_h,u_h) \) (or an equivalent form), then the error estimator may also include this term (or a computable equivalent) for a reliable and efficient a posteriori error control. Amongst the many skeletal schemes like (nonconforming) virtual elements, hybridized (weak) discontinuous Galerkin schemes et al., the HHO methodology has a clear and efficacious stabilization

with the abbreviation \(S_{TF}v_h :=\Pi _{F,k} \left( v_{\mathcal {T}} + (1-\Pi _{T,k})R v_{h} \right) |_{T}-v_{\mathcal {F}}|_{F} \) for \(v_h=(v_\mathcal {T}, v_\mathcal {F})\in V_h\) in terms of the \(L^2\) projections \(\Pi _{K,k}\) onto polynomials of degree at most k on a facet or simplex \(K\in \mathcal {F}\cup \mathcal {T}\) of diameter \(h_K = \textrm{diam}(K)\); cf. Sect. 1.4 for further details. The original residual-based estimator \(\eta _{\textrm{HHO}}\) from the textbook [25] for the Poisson model problem \(-\Delta u = f\) includes (2) and an interpolation \({\mathcal {A}} Ru_h\in V\) of \(Ru_h\) by nodal averaging in

(Multiplicative constants are undisplayed in this introduction for simplicity.) The results from Theorem 4.3 and 4.7 in [25] show reliability and efficiency for the total error (1) and piecewise polynomial source terms \(f\in P_{k+1}(\mathcal {T})\),

1.2 Stabilization-free a posteriori error control

There are objections against the double role of \(s_h(u_h,u_h)\) on both sides of the efficiency and reliability estimate. First, the term \(s_h(u,u_h) \) may dominate both sides of the error estimate. In other words, the total error might be equivalent to \(s_h(u_h,u_h)\), but the quantity of interest may exclusively be

Second, since the stabilization (2) incorporates a negative power of the mesh-size, a reduction property for local refinements remains unclear but is inevitable in the proofs of optimal convergence of an adaptive algorithm [5, 16]. This paper, therefore, asks a different question about the control of the error without the stabilization term (2) in the upper bound and introduces two stabilization-free error estimators (multiplicative constants are undisplayed)

for some parameter \(p \in {\mathbb {N}}_0\). The explicit residual-based a posteriori error estimator \(\eta _{\textrm{res}}\) follows from the a posteriori methodology in the spirit of [14, 18, 20, 21] with a piecewise volume residual \(f+\Delta _\textrm{pw}Ru_h\) and the jumps \([\nabla _{\textrm{pw}}Ru_h]_F\) across a facet F (on the boundary this is only the tangential component of \(\nabla _{\textrm{pw}}Ru_h\)). The equilibrated error estimator \(\eta _{\textrm{eq},p}\) includes the post-processed quantity \(Q_p \in RT_{k+p}(\mathcal {T})\) in the space \( RT_{k+p}(\mathcal {T})\) of Raviart-Thomas functions of degree \(k+p\) for \(p \in {\mathbb {N}}_0\) and the nodal average \({\mathcal {A}} R u_h \in S^{k+1}_0({\mathcal {T}})\) of \(R u_h\). The main results establish reliability and efficiency

for any \(p, q\in {\mathbb {N}}_0\) up to data-oscillations \(\textrm{osc}_{q}^2(f, \mathcal {T}):=\Vert h_\mathcal {T}(1-\Pi _{q})f\Vert ^2_{L^2(\Omega )}\) of the source term \(f\in L^2(\Omega )\) and without any stabilization terms. Computational benchmarks with adaptive mesh-refinement driven by any of these estimators provide numerical evidence for optimal convergence rates.

1.3 Further contributions and outline

The higher-order Crouzeix-Raviart finite element schemes are complicated at least in 3D [23] and then the HHO methodology is an attractive alternative even for simplicial triangulations with partly unexpected advantages like the computation of higher-order guaranteed eigenvalue bounds [15]. Higher convergence rates rely on an appropriate adaptive mesh-refining algorithm and hence stabilization-free a posteriori error estimators are of particular interest. The recent paper [36] establishes the latter for virtual elements with an over-penalization strategy as an extension of [8] for the discontinuous Galerkin schemes. A disadvantage is the quantification of the restriction on the stabilization parameter in practise and poor condition for larger parameters. The stabilization-free a posteriori error control in this paper is based on two observations for the HHO schemes on simplicial triangulations. First, the \(P_1\)-conforming finite element functions let the stabilization vanish and, second, the divergence-free lowest-order Raviart-Thomas functions are \(L^2\) perpendicular to the piecewise gradients \(\nabla _{\textrm{pw}}R u_h\). In fact, those two fairly general properties lead in Sect. 2 to a reliable explicit residual-based a posteriori error estimator. In contrast to the simplified introduction above, the paper also focuses on multiplicative constants that lead to the GUB

cf. Table 1 for explicit quantities and Theorem 1 and Theorem 3 for further details.

Numerical comparisons of \(\eta _{\textrm{HHO}}\) with \(\eta _{\textrm{res}}\) and \(\eta _{\textrm{eq},p}\) favour the latter. Section 2 identifies general building blocks of the a posteriori error analysis for discontinuous schemes with emphasis on explicit constants. An application to HHO leads to the new stabilization-free residual-based estimator \(\eta _{\textrm{res}}\) in Sect. 3. The alternative stabilization-free error estimator \(\eta _{\textrm{eq},p}\) follows from an equilibration strategy plus post-processing in Sect. 4. This paper also contributes to the HHO literature a local equivalence of two stabilizations and the efficiency of the stabilization terms up to data-oscillations in extension of [32]. Numerical comparisons of the different error estimators and an error estimator competition for guaranteed error control of the piecewise energy norm in 2D conclude this paper in Sect. 5. Three computational benchmarks provide striking numerical evidences for the optimality of the associated adaptive algorithms. The appendix provides algorithmic details on the computation of the post-processed contribution \(\Vert Q_p-\nabla _{\textrm{pw}}Ru_h\Vert _{L^2(\Omega )}\) in \(\eta _{\textrm{eq}, p}\).

1.4 Overall notation

Standard notation for Sobolev and Lebesgue spaces and norms apply with \(\Vert \bullet \Vert :=\Vert \bullet \Vert _{L^2(\Omega )}\) and \(|\!|\!|\bullet |\!|\!|:=\Vert \nabla \bullet \Vert _{L^2(\Omega )}\). In particular, \(H(\mathop {\textrm{div}}\limits , \Omega )\) is the space of Sobolev functions with weak divergence in \(L^2(\Omega )\) and \(H(\mathop {\textrm{div}}\limits =0, \Omega )\) contains only divergence-free functions in \(H(\mathop {\textrm{div}}\limits , \Omega )\). Throughout this paper, \(\mathcal {T}\) denotes a shape-regular triangulation of the polyhedral bounded Lipschitz domain \(\Omega \subset \mathbb R^n\) into n-simplices with facets \(\mathcal {F}\) (edges for \(n=2\) and faces for \(n=3\)) and vertices \({\mathcal {V}}\). Let \(\mathcal {F}(\Omega )\) (resp. \({\mathcal {V}}(\Omega )\)) denote the set of interior facets (resp. vertices) and (resp. ). Given \(v \in \Omega \rightarrow {\mathbb {R}}^n\) and \(w: \Omega \rightarrow {\mathbb {R}}^{2n-3}\), let if \(n = 2\) and and \(\mathop {\textrm{Curl}}\limits w :=\mathop {\textrm{curl}}\limits w\) if \(n = 3\). For \(s\in \mathbb R\), let \(H^s(\mathcal {T}), H(\mathop {\textrm{div}}\limits ,\mathcal {T})\), and \(H(\mathop {\textrm{curl}}\limits ,\mathcal {T})\) denote the space of piecewise Sobolev functions with restriction to \(T \in \mathcal {T}\) in \(H^s(T), H(\mathop {\textrm{div}}\limits ,T)\), and \(H(\mathop {\textrm{curl}}\limits ,T)\). To simplify notation, \(H^s(K)\) abbreviates \(H^s(\textrm{int}(K))\) for the open interior \(\textrm{int}(K)\) of a compact set K. The \(L^2\)-scalar product reads \((\bullet , \bullet )_{L^{2}(\omega )}\) for volumes \(\omega \subseteq \Omega \) and \(\langle \bullet , \bullet \rangle _{L^2(\gamma )}\) for surfaces \(\gamma \subset {\overline{\Omega }}\) of co-dimension one; the same symbol applies to scalars and to vectors. For and , let \(\langle \bullet , \bullet \rangle \) denote the duality-brackets in \(V^* \times V\) for the dual space \(V^*\) of V equipped with the operator norm .

Define the energy scalar product \(a(v,w) :=(\nabla v, \nabla w)_{L^2(\Omega )}\) for \(v, w \in H^{1}(\Omega )\) and its piecewise version \(a_{\textrm{pw}}(v, w)=(\nabla _{\textrm{pw}}v,\nabla _{\textrm{pw}}w)_{L^2(\Omega )}\) for \(v, w\in H^1(\mathcal {T})\). The latter induces the seminorm \(|\!|\!|\bullet |\!|\!|_\textrm{pw}:=a_\textrm{pw}(\bullet ,\bullet )^{1/2}\) in \(H^1({\mathcal {T}})\). Here and throughout the paper, \(\nabla _{\textrm{pw}}, \textrm{div}_{\textrm{pw}}, \textrm{curl}_{\textrm{pw}}, {\Delta }_{\textrm{pw}}\), denote the piecewise evaluation of the differential operators \(\nabla \), div, \(\textrm{curl}_{\textrm{pw}}, \Delta \) without explicit reference to the underlying shape-regular triangulation \(\mathcal {T}\).

The vector space \(P_k(K)\) of polynomials of degrees at most \(k\in {\mathbb N}_0\) over a facet or simplex \(K\in \mathcal {F}\cup \mathcal {T}\) defines the piecewise polynomial spaces

and the space of piecewise Raviart-Thomas functions

The associated \(L^2\) projections read \(\Pi _{K,k}:L^2(\Omega )\rightarrow P_k(K), \Pi _k: L^2(\Omega ) \rightarrow P_k(\mathcal {T}),\) and \(\Pi _{\mathcal F, k}: L^2(\Omega ) \rightarrow P_k(\mathcal {F})\) with the convention \(\Pi _{-1}:=0\). Abbreviate \(S^{k+1}_0(\mathcal {T}) :=P_{k+1}(\mathcal {T}) \cap V\) and \(RT_k(\mathcal {T}) :=RT^\textrm{pw}_k(\mathcal {T}) \cap H(\mathop {\textrm{div}}\limits ,\Omega )\) for all \(k \in {\mathbb {N}}_0\). The piecewise constant mesh-size function \(h_\mathcal {T}\in P_0(\mathcal {T})\) satisfies \(h_{\mathcal {T}|T}:=h_T\) for \(T\in \mathcal {T}\) with the diameter \(h_K :=\textrm{diam}(K)\in P_0(K)\) of \(K\in \mathcal {F}\cup \mathcal {T}\).

If not explicitly stated otherwise, constants are independent of the mesh-size in the triangulation but may depend on the shape-regularity and on the polynomial degree k. The abbreviation \(A\lesssim B\) hides a generic constant C (independent of the mesh-size) in \(A\le C\; B; A\approx B\) abbreviates \(A\lesssim B\lesssim A\).

This section investigates general building blocks of the a posteriori error analysis and revisits arguments from [14, 18, 20, 21] with emphasis on multiply connected domains \(\Omega \subset \mathbb R^n\) for \(n=2,3\). The general setting of this section results in reliability for an error estimator that is applicable beyond the HHO methodology. Consider the weak solution \(u\in V = H^1_0(\Omega )\) to the Poisson model problem \(-\Delta u=f\) a.e. in \(\Omega \) and \(u = 0\) on \(\partial \Omega \) for a given source \(f\in L^2(\Omega )\); i.e., \(u\in V\) satisfies

An approximation of the gradient gives rise to the residual seen as a linear functional on V, i.e.,

Let \(\nu _T\) denote the unit outer normal along the boundary \(\partial T\) of each simplex \(T\in \mathcal {T}\) and fix the orientation of a unit normal \(\nu _F=\pm \nu _T\) for each facet \(F\in \mathcal {F}(T)\) of T such that it matches the outer unit normal \(\nu \) of \(\partial \Omega \) at the boundary. The jump \([{G}]_F\) of a piecewise function in \(m\in {\mathbb N}\) components \({G}\in H^1(\mathcal {T};\mathbb R^m)\) reads \([{G}]_F:={G}_{|T_+}-{G}_{|T_-}\) on interior facets \(F=T_+\cap T_-\in \mathcal {F}(\Omega )\) (with \(T_\pm \) labelled such that \(\nu _{T_+|F} = \nu _F = -\nu _{T_-|F})\) and \([{G}]_F:={G}\) on the boundary \(F\in \mathcal {F}(\partial \Omega )\). The main result of this section establishes the residual-based error estimator

as a GUB \(\Vert \nabla u - {G}\Vert \le \eta (\mathcal {T}, {G})\) under minimal assumptions on the approximation \({G}\in H^1(\mathcal {T};\mathbb R^n)\subset L^2(\Omega ; \mathbb R^n)\). The constants \(C_1, C_2\), and \(C_H\) (or upper bounds thereof) are computable; cf. Table 1 for an example in 2D with details in Example 1 at the end of Sect. 2. The first assumption is a weakened discrete solution property

The second assumption is the orthogonality to the lowest-order divergence-free Raviart-Thomas functions

Theorem 1

(residual-based GUB) Suppose that \({G}\in H^1(\mathcal T;\mathbb R^n)\) and \(f\in L^2(\Omega )\) satisfy (5)–(6). Then the error estimator \(\eta (\mathcal {T}, {G})\) from (4) is a GUB

of the error \(\Vert \nabla u-{G}\Vert \) for the solution \(u\in V\) to (3). The constants \(C_\textrm{1}, C_{2}, C_{\textrm{H}}\) exclusively depend on \(\Omega \) and the shape-regularity of \(\mathcal T\).

The remaining parts of this section are devoted to the proof of Theorem 1 and the computation of (upper bounds of) the constants \(C_\textrm{1}, C_\textrm{2}\), and \(C_{\textrm{H}}\) in (4). The point of departure is the subsequent decomposition that appears necessary in the nonconforming and mixed finite element a posteriori error analysis. It leads to a split of the error \(\Vert \nabla u-{G}\Vert \) into some divergence part and some consistency part.

Lemma 1

(decomposition) Any \(v\in V\) and \({G}\in L^2(\Omega ; \mathbb R^n)\) satisfy the decomposition

(7)

with the (unique) minimizer \(w\in V\) of the distance

of \({G}\) to the gradients \(\nabla V\) of Sobolev functions. The solution \(u\in V\) to (3) satisfies

Proof

The minimizer \(w \in V\) of \(\Vert {G}- \nabla \varphi \Vert \) among \(\varphi \in V\) satisfies the variational formulation \(a(w, \varphi )_{L^2(\Omega )} = ({G},\nabla \varphi )_{L^2(\Omega )}\) for all \(\varphi \in V\). (Notice that w is the unique weak solution to the Poisson model problem \(-\Delta w = -\mathop {\textrm{div}}\limits {G}\in V^*\).) In particular, \({G}- \nabla w\) is \(L^2\) orthogonal onto \( \nabla V\) and the Pythagoras theorem proves (7). Given \(\varphi \in V\) with \( |\!|\!|\varphi |\!|\!| =1\), the orthogonality of \({G}-\nabla w\) to \(\nabla \varphi \) and (3) show

with the duality brackets \(\langle \bullet , \bullet \rangle \) in \(V^* \times V\). Since the supremum of (9) over all \(\varphi \in V\) with \( |\!|\!|\varphi |\!|\!| = 1\) is equal to \(|\!|\!|u-w|\!|\!|= |\!|\!| f+\mathop {\textrm{div}}\limits {G}|\!|\!|_{*} \), this and (7) conclude the proof of (8). \(\square \)

The split (7) of the error \(\Vert \nabla u-{G}\Vert \) allows for and enforces a separate estimation of the equilibrium and consistency contribution in residual-based a posteriori error estimators.

In order to derive explicit constants, two lemmas are recalled. The first has a long tradition in the a posteriori error control in form of a Helmholtz decomposition on simply connected domains [4, 13] and introduces the constant \(C_{\textrm{H}}\) from Theorem 1. The following version includes the general case of multiply connected domains as in [33] for \(n=2\) or \(n=3\) dimensions and weak assumptions on a divergence-free function \(\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )\).

Lemma 2

(Helmholtz-decomposition) Suppose the divergence-free function \(\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )\) is \(L^2\) orthogonal onto \(RT_0(\mathcal T)\cap H(\mathop {\textrm{div}}\limits =0, \Omega )\). Then there exists \(\beta \in H^1(\Omega ;\mathbb R^N), N=2n-3\), such that any \(\beta _C\in S^1(\mathcal T)^N\) satisfies

The constant \(C_{\textrm{H}}>0\) exclusively depends on \(\Omega \).

Proof

The compact polyhedral boundary \(\partial \Omega \) of the bounded Lipschitz domain \(\Omega \) has \(J+1\) connectivity components \(\Gamma _0,\dots ,\Gamma _J\) for some finite \(J\in \mathbb N_0\). Those connectivity components have a positive surface measure \(|\Gamma _j|\) and a positive distance of each other. So the integral mean

is well defined and depends continuously on \(\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )\) in the sense that \(|\gamma _j|\le c_1\Vert \varrho \Vert \) (recall \(\mathop {\textrm{div}}\limits \varrho =0\)) for each \(j=0,\dots , J\) and \(c_1>0\). This constant \(c_1\) and the constants \(c_2,c_3,c_4\) below exclusively depend on the domain \(\Omega \). The finite real numbers \(\gamma _0,\dots ,\gamma _J\) define the Neumann data for the harmonic function \(z\in H^1(\Omega )/\mathbb R\) with

The elliptic regularity theory for polyhedral domains lead to \(z\in H^{1+\alpha }(\Omega )\) for some \(\alpha >1/2\) and \(c_2>0\) with \(\Vert z \Vert _{H^{1+\alpha }(\Omega )}\le c_2\, (|\gamma _0|+\dots +|\gamma _J|)\). The Raviart-Thomas interpolation operator defines a bounded linear operator on \(H(\mathop {\textrm{div}}\limits ,\Omega )\cap L^p(\Omega ;\mathbb R^n)\) for \(p>2\). It is generally accepted that, for \(\alpha >0\) and \( \nabla z\in H(\mathop {\textrm{div}}\limits =0,\Omega )\cap H^\alpha (\Omega ;\mathbb R^n)\), the Fortin interpolation \(I_{\textrm{F}} \nabla z \in RT_0(\mathcal T)\cap H(\mathop {\textrm{div}}\limits =0,\Omega )\) is well defined and \(\Vert I_{\textrm{F}} \nabla z\Vert \le c_3 \Vert \nabla z\Vert _{H^\alpha (\Omega )}\) follows for some \(c_3>0\). The additional property \( \nabla z\in L^p(\Omega ;\mathbb R^n)\) for some \(p>2\) allows the definition of \(\int _F \nabla z\cdot \nu _F \ \textrm{d} s\) as a Lebesgue integral over a facet \(F\in \mathcal {F}\). One consequence for the boundary facets is the vanishing integral

Since \(\varrho -I_{\textrm{F}} \nabla z\in H(\mathop {\textrm{div}}\limits =0,\Omega )\) is divergence-free, Theorems 3.1 and 3.4 in [33] prove the existence of \(c_4>0\) and \(\beta \in H^1(\Omega ;\mathbb R^N)\) with

Recall that \(\varrho \perp I_{\textrm{F}}\nabla z\) and \(\varrho \perp \mathop {\textrm{curl}}\limits \beta _C\in RT_0(\mathcal {T})\cap H(\mathop {\textrm{div}}\limits =0, \Omega )\). This concludes the proof of (10) with \(C_{\textrm{H}}:=\sqrt{1 + (c_1c_2c_3(1+J))^2}c_4\). \(\square \)

The subsequent version of the trace inequality on the facets \(\mathcal {F}\) leads to the piecewise constant \(\ell \in P_0(\mathcal {F})\) defined by

Lemma 3

(trace inequality) Any \(f\in H^1(\Omega )\) satisfies

with the constant \(C_{\textrm{tr}}:=\max _{T\in \mathcal {T}}\max _{x\in T}|x-\textrm{mid}(T)|/h_T<n/(n+1)\).

Proof

The center of inertia \(\textrm{mid}(T)=\sum _{j=0}^n P_j/(n+1)\) of the n-simplex \(T=\textrm{conv}\{P_0,..., P_n\}\in \mathcal {T}\) and the \(n+1\) faces \(F_j=\textrm{conv}\{P_0,..., P_{j-1},P_{j+1},..., P_n\}\in {\mathcal {F}}(T)\) for \(j=0,..., n\) give rise to the decomposition of T into \(n+1\) sub-simplices \(T_j' = \textrm{conv}(F_j, \textrm{mid}(T))\) with volume

\(|T_j'|=|T|/(n+1)\).

Standard arguments like the trace identity on \(T_j'\subset T\) [17, Lemma 2.1] for \(|f|^2\in W^{1,1}(T')\) and a Cauchy inequality show

The distance \(\max _{x\in T_j'}|x-\textrm{mid}(T)|=|P_k - \textrm{mid}(T)|\le C_{\textrm{tr}}h_T\) is attained at a vertex \(P_k\) for \(k\in \{0,..., j-1, j+1,..., n\}\). Since the centroid \(\textrm{mid}(T)\) divides each median of T in the ratio n to 1 and the length of each median is strictly bounded by \(h_T\),

the bound \(C_{\textrm{tr}}<n/(n+1)\) follows and cannot be improved in the absence of further assumptions on the shape of the simplex T. Since \(|T_j'| = |T|/(n+1)\), the previously displayed estimate leads to

Let \(\mathcal {T}'\) be the refinement of \(\mathcal {T}\), obtained by replacing \(T\in \mathcal {T}\) with \(T_0',..., T_d'\) from above. The triangulation \(\mathcal {T}'\) allows for the facet based decomposition \(\{\omega '(F)\}_{F\in \mathcal {F}}\) of \(\Omega \), where \(\omega '(F)\) is either the patch \(\omega '(F)=\textrm{int}(T'_+\cup T'_-)\) for an interior facet \(F=T'_+\cap T'_-\) or \(\omega '(F)=\textrm{int}(T')\) for \(F\in \mathcal {F}(\partial \Omega )\cap \mathcal {F}(T')\). This establishes, for any \(F\in \mathcal {F}\), the estimate

Since the family \(\{\omega '(F): F\in \mathcal {F}\}\) has no overlap, the sum of the last displayed inequality over all \(F \in \mathcal {F}\) and a Cauchy inequality conclude the proof of Lemma 3. \(\square \)

The next lemma utilizes a quasi-interpolation operator \(J: H^1(\Omega ) \rightarrow S^1(\mathcal {T})\) with the restriction \(J(V) \subset S^1_0(\mathcal {T})\), e.g., \(J = J_1 \circ I_{NC}\) from [19, Section 5] with explicit constants for \(n=2\), and the approximation and stability properties

for constants \(C_{1}\) and \(C_\textrm{st}\) exclusively depending on the shape-regularity of \(\mathcal {T}\). For the precise definition of \(J_1\) and \(I_{NC}\), we refer to [19, eq. (47) and Section 5]. Recall the constant \(C_{\textrm{tr}}\) from Lemma 3 and set \(C_{2}:=(C_{1}(C_{1}+2C_{\textrm{tr}}C_\textrm{st}\,/n))^{1/2}\).

Lemma 4

(equilibrium) Suppose that \({G}\in H(\mathop {\textrm{div}}\limits ,\mathcal T)\) and \(f\in L^2(\Omega )\) satisfy (5) and suppose \(({G}|_T)|_F\cdot \nu _F\in L^2(F)\) for all \(F\in {\mathcal {F}}(T)\) and \(T\in \mathcal T\). Then

Proof

Given \(\varphi \in V\) with \( |\!|\!|\varphi |\!|\!| =1\), set \(\psi :=\varphi -\varphi _C\) for some quasi-interpolation \(\varphi _C:=J\varphi \in S^1_0(\mathcal T)\) with (11). Since (5) implies \(\langle f+\mathop {\textrm{div}}\limits {G},\varphi \rangle =\langle f+\mathop {\textrm{div}}\limits {G},\psi \rangle \), a piecewise integration by parts and the collection of jump contributions show

The first bound follows from a Cauchy inequality and (11),

The second bound additionally exploits the trace inequality of Lemma 3,

with \({\left| \hspace{-1.0625pt}\left| \hspace{-1.0625pt}\left| \varphi \right| \hspace{-1.0625pt}\right| \hspace{-1.0625pt}\right| }=1\) in the last step. Since (13)–(14) hold for all \(\varphi \in V\) with \(|\!|\!|\varphi |\!|\!|=1\), the supremum in (12) over all such \(\varphi \) concludes the proof.\(\square \)

The final ingredient for the proof of Theorem 1 controls the second term \(\delta \) in the decomposition of Lemma 1 for \(\varrho :={G}-\nabla w\). Recall \(C_{\textrm{H}}\) from Lemma 2 and \({C_{1}}\) from (11), and \(C_{2}\) from page 9.

Lemma 5

(conformity) Suppose the divergence-free function \(\varrho \in H(\mathop {\textrm{div}}\limits =0,\Omega )\cap H(\mathop {\textrm{curl}}\limits ,\mathcal T)\) is \(L^2\) orthogonal onto \(RT_0(\mathcal T) \cap H(\mathop {\textrm{div}}\limits =0,\Omega )\) and satisfies \((\varrho |_T)|_F\times \nu _F\in L^2(F)\) for all \(F\in {\mathcal {F}}(T)\) and \(T\in T\). Then

Proof

Lemma 2 provides \(\beta \in H^1(\Omega ; \mathbb R^N)\) with (10) for a (component-wise) quasi-interpolation \(\beta _C\in S^1(\mathcal {T})^N\) with (11) as in the proof of Lemma 4; set \(\psi :=\beta -\beta _C\). A piecewise integration by parts and the collection of jump contributions shows

Stability and approximation properties of the quasi-interpolation (11) and the trace inequality of Lemma 3 eventually lead to

In fact, the routine estimation with element and jump terms is completely analogous to the proof of Lemma 4 and leads to the same constants \({C_{1}}, C_{2}\). This and \(|\!|\!|\beta |\!|\!|\le C_{\textrm{H}}\Vert \varrho \Vert \) conclude the proof. \(\square \)

Proof

(Theorem 1) The trace of \({G}|_T\in H^1(T; \mathbb R^n)\) is well defined on any facet \(F\in \mathcal {F}(T)\) of the simplex \(T\in \mathcal {T}\). Lemma 1 provides \(w\in V=H^1_0(\Omega )\) with \(\Vert \nabla u-{G}\Vert ^2 = |\!|\!| f + \mathop {\textrm{div}}\limits {G}|\!|\!|_{*}^2+\Vert {G}- \nabla w\Vert ^2\). Since \({G}\) satisfies (5), Lemma 4 establishes

The assumption (6) on \({G}\) and an integration by parts prove that Lemma 5 is applicable to \(\varrho :={G}-\nabla w\in H(\mathop {\textrm{div}}\limits =0, \Omega )\cap H(\mathop {\textrm{curl}}\limits , \mathcal {T}; \mathbb R^n)\). Since \(\mathop {\textrm{curl}}\limits \nabla w = 0\) and \(\nabla w\times \nu = 0\), this reveals

The above estimates together with the decomposition of Lemma 1 establish \(\eta (\mathcal {T}, {G})\) as a GUB for the error \(\Vert \nabla u - {G}\Vert \). \(\square \)

Example 1

(constants for right-isosceles triangles) In two space dimensions, \(\Vert \mathop {\textrm{Curl}}\limits \bullet \Vert = |\!|\!|\bullet |\!|\!|\) and so \(C_{\textrm{H}}\le 1\) for a simply connected domain \(\Omega \) in Lemma 2. The choice \(J:=J_1\circ I_{\textrm{NC}}\) from [19, Section 5] of the quasi-interpolation operator J in the proof of Lemma 4 allows for the explicit estimates

where \(j_{1,1} = 3.8317\) denotes the first positive root of the first Bessel function. For triangulations into right-isosceles triangles, the constant \(c_{\textrm{apx}}\le \sqrt{3}/(2 - 2\cos (\pi /\max \{4,M_{\textrm{bd}}\}))\) from [19, Lemma 4.8] depends on the domain by the maximal number \(M_{\textrm{bd}}\le 4\max \{\pi ,\omega _{\textrm{max}}\}/\pi \le 8\) of triangles sharing a boundary vertex. Given the maximal interior angle \(\omega _{\textrm{max}}\) of \(\Omega \), Table 1 displays those constants for the maximal possible value \(M_{\textrm{bd}}= 4\max \{\pi ,\omega _{\textrm{max}}\}/\pi \). The geometric quantity \(\max _{x\in T}|x-\textrm{mid}(T)|\) equals two-thirds of the maximum median of T. Thus, \(C_{\textrm{tr}}=\sqrt{5}/(3\sqrt{2})\le 0.5271\) and \(\ell (F) = 6h_F\) for interior edges \(F\in \mathcal {F}(\Omega )\) and \(\ell (F)=12h_F\) for boundary edges \(F\in \mathcal {F}(\partial \Omega )\) of triangulations into right-isosceles triangles. Consequently,

The arguments from Sect. 2 apply to the HHO method and result in a stabilization-free reliable a posteriori error control. In combination with the efficiency estimate from this section, this leads to a new explicit residual-based a posteriori error estimator for the HHO method that is equivalent to the error up to data oscillations.

3.1 Hybrid high-order methodology

The HHO ansatz space reads \( V_h:=P_k(\mathcal {T}) \times P_k(\mathcal {F}(\Omega ))\) for \(k\in {\mathbb N}_0\) with the subspace \(P_k(\mathcal {F}(\Omega )) \subset P_k(\mathcal {F})\) of piecewise polynomials \(p\in P_k(\mathcal {F})\) under the convention \(p_{|\partial \Omega }=0\). The interpolation \(\textrm{I}: V \rightarrow V_h\) maps \(v \in V\) onto \(\textrm{I} v :=(\Pi _k v, \Pi _{{\mathcal {F}},k} v) \in V_h\). Given any \(v_h = (v_{\mathcal {T}},v_{\mathcal {F}})\in V_h\), the reconstruction operator \(R: V_h\rightarrow P_{k+1}(\mathcal {T})\) defines the unique piecewise polynomial \(Rv_h \in P_{k+1}(\mathcal {T})\) with \(\Pi _{0}(Rv_h - v_{\mathcal {T}}) = 0\) such that, for all \(w_{k+1} \in P_{k+1}(\mathcal {T})\),

Let \(u_h\in V_h\) solve the HHO discrete formulation of (3) with

for the HHO bilinear form

and the stabilization term \(s_h(u_h, v_h)\) from (2). Given any \(w_C \in S^1_0(\mathcal {T})= P_1(\mathcal {T})\cap H^1_0(\Omega )\), the definition of the reconstruction operator R in (17) verifies \(R \textrm{I} w_C=w_C\) with the interpolation \(\textrm{I}\) onto \(V_h\). Hence, \(S_{TF}\textrm{I} w_C=0\) vanishes for all \(F \in \mathcal {F}(T)\) and \(T \in \mathcal {T}\). This and (18) show, for all \(w_C\in S^1_0(\mathcal {T})\), that

3.2 Explicit a posteriori error estimator

As a result of (20), \(\nabla _{\textrm{pw}}Ru_h\) satisfies the solution property (5) if \(k\ge 1\) and, in the lowest order case \(k=0\), (5) holds with f replaced by \(\Pi _0 f\). This allows the application of the theory from Sect. 2 to the HHO method with minor modifications for the case \(k=0\). Define the error estimator contributions

Since \(\nabla _{\textrm{pw}}Ru_h\) is a piecewise gradient, its piecewise \(\mathop {\textrm{curl}}\limits \) vanishes. This leads to the explicit residual-based a posteriori error estimator

(Recall \({C_{1}}, C_{2} \) from Lemma 4 and \(C_{\textrm{H}}\) from Lemma 2 as well as the Poincaré constant \(C_P\le \pi ^{-1}\).) The main result of this section verifies the assumptions in Theorem 1 and proves reliability and efficiency of \(\eta _{\textrm{res}}(\mathcal {T})\).

Theorem 2

(residual-based GUB for HHO) Let \(u\in V\) solve the Poisson equation (3) and let \(u_h \in V_h\) solve the discrete formulation (18). Then

and \(\textrm{osc}_{k-1}(f, \mathcal T)\le C_{4} \eta _{\textrm{res}}(\mathcal {T})\) hold for any \(q\in {\mathbb N}_0\). The constants \(C_{3}\) and \(C_{4}\) exclusively depend on k, q and on the shape-regularity of the triangulation \(\mathcal T\).

3.3 Proof of Theorem 2

The orthogonality of \(\nabla _\textrm{pw}R u_h\) to the divergence-free Raviart-Thomas space of lowest degree is an assumption in Theorem 1 and verified below.

Lemma 6

(orthogonality) The piecewise gradients \(\nabla _{\textrm{pw}}R V_h\) are \(L^2\) orthogonal to the space  \(RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )\), i.e., any \(v_h\in V_h\) and \(q_{RT}\in RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )\) satisfy

Proof

Given any \(q_{RT}\in RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )\), \(\mathop {\textrm{div}}\limits q_{RT}=0\) shows \(q_{RT}\in P_0(\mathcal {T};\mathbb R^n)\) [28, Lemma 14.9]. Since \(P_0(\mathcal {T};\mathbb R^n) = \nabla _{\textrm{pw}}P_1(\mathcal {T})\), there exists a piecewise affine function \(\phi _1 \in P_1(\mathcal {T})\) with \(q_{RT}= \nabla _\textrm{pw}\phi _1\) a.e. in \(\Omega \). This and the definition of \( Rv_h\) from (17) imply, for any \(v_h = (v_{\mathcal {T}},v_{\mathcal {F}})\in V_h\), that

This, a piecewise integration by parts, and \(\Delta _\textrm{pw}\phi _1 \equiv 0\) lead to

Since \(q_{RT}\in RT_0(\mathcal {T})\) has continuous normal components, the jump term \(\langle v_{\mathcal {F}}, [q_{RT}]_F\cdot \nu _F\rangle _{L^2(F)} =0\) vanishes for all \(F\in \mathcal F(\Omega )\). This, \(v_{\mathcal {F}}\equiv 0\) on \(\partial \Omega \), and (24) conclude \((\nabla _{\textrm{pw}} Rv_h,q_{RT})_{L^2(\Omega )} = 0\). \(\square \)

The following lemma concerns the efficiency of the jump contributions. Each facet \(F\in \mathcal {F}\) has at most two adjacent simplices that define a triangulation \(\mathcal {T}(F):=\{T\in \mathcal {T}: F\in \mathcal {F}(T)\}\) of the facet-patch \(\omega (F):=\textrm{int}(\bigcup _{T\in \mathcal {T}(F)}T)\).

Lemma 7

(efficiency of jumps) The solution \(u \in V\) to (3) and the discrete solution \(u_h \in V_h\) to (18) satisfy (a) for all \(F \in \mathcal {F}\) and (b) for all \(F \in \mathcal {F}(\Omega )\).

1. (a)

   \(h_F^{1/2}\Vert [\nabla _{\textrm{pw}} Ru_h]_F\times \nu _F \Vert _{L^2(F)} \lesssim \min _{v \in V} \Vert \nabla v - \nabla _{\textrm{pw}} Ru_h\Vert _{L^2(\omega (F))}\),

2. (b)

   \(h_F^{1/2}\Vert [\nabla _{\textrm{pw}} Ru_h]_F\cdot \nu _F \Vert _{L^2(F)} \lesssim \Vert \nabla u - \nabla _{\textrm{pw}} Ru_h\Vert _{L^2(\omega (F))} + \textrm{osc}_k(f, \mathcal {T}(F))\).

Proof

The proof is based on the following extension argument. Given a polynomial \(p\in P_k(F)\) of degree at most k along the side \(F\in \mathcal {F}\), the coefficients determine a polynomial (also denoted by p) along the hyperplane H that enlarges F. The intersection \({\widehat{F}}:= H\cap \textrm{conv}({ \omega (F)}) \) of the hyperplane H with the convex hull of the facet-patch \(\omega (F)\) may be strictly larger than F. The shape-regularity of \(\mathcal {T}\) bounds the size of \({\widehat{F}}\) in terms of F and an inverse estimate leads to a bound \(\Vert p \Vert _{L^\infty ({\widehat{F}})} \le C(k)\Vert p \Vert _{L^\infty (F)}\) with a constant C(k) that depends on the shape-regularity of \(\mathcal {T}\) and on k. The extension of p from H to \(\mathbb R^n\) by constant values along the side normal \(\nu _F\) leads to a polynomial \({\widehat{p}}\in P_k(\mathbb R^ n)\) with

Proof of (a). The tangential jump \(\varrho _F:= [\nabla _{\textrm{pw}} Rv_h]_F\times \nu _F \in P_{k}(F;{\mathbb {R}}^N)\) is a polynomial in \(N=2n-3\) components on \(F\in \mathcal {F}\) for \(n=2,3\). Let \(p=\varrho _F(j)\) be one of the components of \(\varrho _F\in P_k(F)^N\), for \(j=1,...,N\), and extend it as explained above to \({\hat{p}}\in P_k(\mathbb R^n)\) and call this \({\widehat{\varrho }}(j)\) in the vector \(\widehat{\varrho _F}\in P_k(\mathbb R^n;{\mathbb {R}}^N)\). The proof involves the piecewise polynomial facet-bubble function \(b_F:=n^n \Pi _{j=1}^n \varphi _j \) for the n nodal basis function \(\varphi _1, \dots , \varphi _n \in S^1(\mathcal T)\) associated with the vertices of F. An inverse estimate [38, Proposition 3.37] shows

Since \(\varrho := b_F\widehat{\varrho _F} \in S^{k+n}_0(\mathcal {T}(F);{\mathbb {R}}^N)\) vanishes on \(\partial \omega (F){\setminus }\textrm{int}(F)\), (26) and a piecewise integration by parts show

This and \((\mathop {\textrm{curl}}\limits \varrho , \nabla v)_{L^2(\Omega )} = 0\) for any \(v\in V\) imply

An inverse estimate, \(\Vert b_F\Vert _{L^\infty (\omega (F))} = 1\), and (25) imply

In combination with (28), this concludes the proof of (a).

Proof of (b). The efficiency of normal jumps (b) follows from the arguments for conforming FEMs, cf. [38, Section 1.4.5]; further details are omitted. \(\square \)

The following lemma reveals that the order \(k \ge {\mathbb {N}}_0\) of the oscillations \(\textrm{osc}_{k}(f, T)\) in Lemma 7 (b) can be any natural number. It is certainly known to the experts but hard to find in the literature. Recall the convention \(\Pi _{-1}:=0\).

Lemma 8

(efficiency of lower-order oscillations) Given any simplex \(T \in {\mathcal {T}}\) and parameters \(k, q \in {\mathbb {N}}_0\), the solution \(u \in V\) to (3) satisfies

The constant \(C_5\) exclusively depends on q and the shape of T.

Proof

The assertion (30) is trivial for \(q \le k-1\), so suppose \(k \le q\). Any \(v_{k+1}\in P_{k+1}(T)\) and \(\varrho _T:=\Pi _q f + \Delta v_{k+1}\in P_q(T)\) satisfy

Let \(b_T\in S^{n+1}_0(T)\) with \(0\le b_T\le 1=\max b_T\) denote the volume bubble-function on \(T\in \mathcal {T}\). The equivalence of norms in the finite-dimensional space \(P_q(T)\) provides

A more detailed analysis of the mass matrices reveals that the constant \(C_{6}\) exclusively depends on the polynomial degree q. An integration by parts with \(b_T\varrho _T\in S^{q+n+1}_0(T)\subset V\) and the weak formulation (3) result in

A Cauchy inequality, the inverse estimate \(h_T\Vert \nabla (b_T\varrho _T)\Vert _{L^2(T)}\le C_7\Vert \varrho _T\Vert _{L^2(T)}\) with a constant \(C_7\) that exclusively depends on \(q+n+1\) and the shape of T, and (32) lead to

The combination of (31) with (33) and a Cauchy inequality conclude the proof of (30), e.g., with \(C_{5}= 1 + {C_{6}}^4(1 + {C_{7}}^2).\) \(\square \)

Proof

(of Theorem 2) Recall the definition of \(\eta _{\textrm{res}}(\mathcal {T})\) for \(k\ge 1\) and \(k=0\) in (21). Since \( \textrm{osc}_{k-1}^2(f, \mathcal {T})\le \eta _{\textrm{res},1}^2(\mathcal {T}) + \eta _{\textrm{res},2}^2(\mathcal {T})\lesssim \eta _{\textrm{res}}^2(\mathcal {T})\), the remaining parts of this proof discuss the reliability and efficiency of \(\eta _{\textrm{res}}(\mathcal {T})\).

Lemma 6 provides the orthogonality of \(\nabla _{\textrm{pw}}Ru_h\in H^1(\mathcal {T}; \mathbb R^n)\) to the divergence-free Raviart-Thomas function \(RT_0(\mathcal {T})\cap H(\textrm{div} = 0,\Omega )\). This and (20) show that the assumptions in Theorem 1 hold for \({G}:=\nabla _{\textrm{pw}}Ru_h\) and \(k\ge 1\), whence the reliability of \(\eta _{\textrm{res}}(\mathcal {T})\) follows with a reliability constant 1.

Minor modifications to the proof of Theorem 1 lead to reliability in the case \(k=0\). In fact, the only modifications required concern the upper bound of \(|\!|\!|f + \textrm{div}\,\nabla _{\textrm{pw}}R u_h|\!|\!|_* \le |\!|\!|(1 - \Pi _0) f|\!|\!|_* + |\!|\!|\Pi _0 f + \textrm{div}\,\nabla _{\textrm{pw}}R u_h|\!|\!|_*\). A piecewise Poincaré inequality shows \(|\!|\!|(1 - \Pi _0) f|\!|\!|_* \le C_P\textrm{osc}_{0}(\mathcal {T},f)\) with the Poincaré constant \(C_P\)(\(\le 1/\pi \) for simplices). Lemma 4 proves \(|\!|\!|\Pi _0 f + \textrm{div}\,\nabla _{\textrm{pw}}R u_h|\!|\!|_* \le C_\textrm{1}\eta _{\textrm{res},1}(\mathcal {T}) + C_\textrm{2}\eta _{\textrm{res},3}(\mathcal {T})\). Hence, the decomposition of Lemma 1 and Lemma 5 result in \(|\!|\!|u-Ru_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{res}}(\mathcal {T})\).

This provides the reliability and it remains to verify the efficiency \(\eta _{\textrm{res}}(\mathcal {T})\lesssim |\!|\!|u-Ru_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,\mathcal {T})\) for any \(q\in {\mathbb {N}}_0\). The Pythagoras theorem and (33) with \(\varrho _T:=\Pi _kf + \Delta R u_h\in P_k(T)\) and \(v_{k+1}:=Ru_h\in P_{k+1}(T)\) lead to the local efficiency of the volume contributions

Lemma 7 considers the remaining terms in the error estimator and establishes their efficiency namely,

with the modified jump \([\nabla _\textrm{pw}R u_h]_F = \nabla _\textrm{pw}R u_h \times \nu _F\) on boundary facets \(F \in {\mathcal {F}}(\partial \Omega )\). This and Lemma 8 establish the existence of some mesh-independent constant \(C_3>0\) with \(C_3^{-1}\eta _{\textrm{res}}(\mathcal {T})\le |\!|\!|u-Ru_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,\mathcal {T})\) for arbitrary \(q\in {\mathbb {N}}_0\). This concludes the proof.\(\square \)

While the focus of this paper is on the HHO methodology, the framework of Sect. 2 also applies to other skeletal methods as well. The following example covers a hybridized discontinuous Galerkin (HDG) FEM from [35] with the Lehrenfeld-Schöberl stabilization.

Example 2

(HDG FEM) Let \(V_h :=P_{k+1}({\mathcal {T}}) \times P_k({\mathcal {F}}(\Omega ))\) for \(k \in {\mathbb {N}}_0\). An equivalent formulation to the HDG FEM from [35] seeks \(u_h \in V_h\) with

Here, \(R: V_h \rightarrow P_{k+1}({\mathcal {T}})\) is defined as in (17) and

for any \(v_h = (v_{\mathcal {T}}, v_{\mathcal {F}}), w_h = (w_{\mathcal {T}}, w_{\mathcal {F}}) \in V_h\). This method is also known under the label of weak Galerkin FEM [39]. It is straightforward to verify that \(\nabla _\textrm{pw}R u_h\) satisfies (5)–(6). Notice that (5) also holds for \(k = 0\) without any modification. Therefore, Theorem 1 leads to the reliable a posteriori estimate

The efficiency \(\eta _{\textrm{res}}({\mathcal {T}}) \lesssim |\!|\!|u-R u_h|\!|\!|_\textrm{pw}\) follows from the arguments in the proofs of Lemma 7–8.

The residual-based guaranteed upper bound (GUB) of the error \(|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\) from Sect. 3.2 employs explicit constants that may lead to overestimation in higher dimensions and for different triangular shapes. This section utilizes equilibrated flux reconstructions [1, 2, 5, 6, 30] to establish, up to the well-known Poincaré constant \(C_P \le 1/\pi \), a constant-free guaranteed upper bounds for a tight error control.

4.1 Guaranteed error control

The guaranteed upper bounds of this section involves two post-processings of the potential reconstruction \(R u_h \in P_{k+1}(\mathcal {T})\) of the discrete solution \(u_h\) to (18). First, the patch-wise design of a flux reconstruction \(Q_p \in RT_{k+p}(\mathcal {T})\) with \(p \in {\mathbb {N}}_0\) from [3, 10, 31] provides an \(H(\textrm{div},\Omega )\)-conforming approximation to \(\nabla _{\textrm{pw}}R u_h\) with the equilibrium \(\Pi _{r} f + {\text { div }}Q_p = 0\) in \(\Omega \) and r from (35) below. Second, the nodal average \({\mathcal {A}} R u_h \in S^{k+1}_0(\mathcal {T})\subset V\) results in an V-conforming approximation of \(R u_h\) by averaging all values of the discontinuous function \(R u_h\) at each Lagrange point of \(S^{k+1}_0(\mathcal {T})\). This, the split (7), and the solution property (20) give rise to the guaranteed upper bound (GUB)

with \(r \in {\mathbb {N}}_0\) defined by

The main result of this section states the reliability and efficiency (up to data oscillations) of \(\eta _{\textrm{eq},p}\) for all parameters \(p \in {\mathbb {N}}_0\).

Theorem 3

(equilibrium-based GUB for HHO) Let \(u \in V\) resp. \(u_h \in V_h\) solve (3) resp. (18). Given a parameter \(p \in {\mathbb {N}}_0\), there exists \(Q_p\in RT_{k+p}({\mathcal {T}})\) such that the error estimator \(\eta _{\textrm{eq},p}(\mathcal {T})\) from (34) is an efficient GUB

for any \(q\in {\mathbb {N}}_0\) and \(\textrm{osc}_{k-1}(f, {\mathcal {T}}) \le {C_{9}}\eta _{\textrm{eq},p}(\mathcal {T})\).

The constants \(C_8\) and \(C_9\) exclusively depend on the polynomial degree \(k \in {\mathbb {N}}_0\), the parameter \(q \in {\mathbb {N}}_0\), and the shape-regularity of \(\mathcal {T}\).

At least two technical contributions for the proof of Theorem 3 are of broader interest. A first contribution to the HHO literature is the local equivalence of the original HHO stabilization \(s_h\) from (2) and the alternative stabilization \({\tilde{s}}_h(v_h,v_h) :=\sum _{T \in \mathcal {T}} {\tilde{s}}_T(v_h,v_h)\) from [25] defined, for \(v_h=(v_\mathcal {T}, v_\mathcal {F})\in V_h\), by

A second result of separate interest in the HHO literature (cf. [25] where the efficiency in (38) is left open) is the efficiency of the stabilizations from Theorems 4–5 below,

The subsequent subsection continues with some explanations on the flux reconstruction \(Q_p\in RT_{k+p}({\mathcal {T}})\) that is defined by local minimization problems on each vertex patch. Appendix A complements the discussion with an algorithmic two-step procedure for the computation of \(Q_{p}-\nabla _{\textrm{pw}}R u_h\) in 2D. The efficiency of the averaging \(|\!|\!|(1 - {\mathcal {A}})R u_h|\!|\!|_\textrm{pw}\) up to data oscillations follows in Sects. 4.3–4.4. Section 4.6 concludes with the proof of Theorem 3.

4.2 Construction of equilibrated flux

This subsection defines the post-processed \(H(\textrm{div},\Omega )\)-conforming equilibrated flux \(Q_p\in RT_{k+p}({\mathcal {T}})\) that enters the GUB \(\eta _{\textrm{eq},p}\) from (34) based on local patch-wise minimization problems in the spirit of [10, 30, 31].

Consider the shape-regular vertex-patch \(\omega (z):=\textrm{int}(\bigcup \{T\in \mathcal {T}(z)\})\) covered by the neighbouring simplices \(\mathcal {T}(z):=\{T\in \mathcal {T}: z\in T\}\) sharing a given vertex \(z \in {\mathcal {V}}\) with the facet spider \(\mathcal {F}(z):=\{F\in \mathcal {F}\,\ z\in F\}\). Recall the space of piecewise Raviart-Thomas functions \(RT_{k}^{\textrm{pw}}(\mathcal {T})\) from Sect. 1.4 and define

Throughout the remaining parts of this section, abbreviate \(G_h:=\nabla _\textrm{pw}R u_h \in P_k({\mathcal {T}};\mathbb R^n)\). Given a vertex \(z\in {\mathcal {V}}\) with the \(P_1\)-conforming nodal basis function \(\varphi _z\in S^1(\mathcal {T})\), the property (20) provides compatible data

such that the discrete affine space

is not empty. Consequently,

is well defined as the \(L^2\) projection \(\Pi _{{\mathcal {Q}}_h(z)} {\mathcal {I}}_{RT} (\varphi _z G_h)\) of \({\mathcal {I}}_{RT} (\varphi _z G_h)\) onto \({\mathcal {Q}}_h(z)\) with the piecewise Raviart-Thomas interpolation \({\mathcal {I}}_{RT}: H^1(\mathcal {T};\mathbb R^n)\rightarrow RT^\textrm{pw}_{k+p}(\mathcal {T})\) [7, Section III.3.1]. In the case \(p \ge 1\), \(\varphi _z G_h\in P_{k+1}({\mathcal {T}}(z)) \subset RT_{k+p}^\textrm{pw}({\mathcal {T}}(z))\) is a piecewise Raviart-Thomas function of degree \(k+p\). Hence \({\mathcal {I}}_{RT}(\varphi _z G_h) = \varphi _z G_h\) and \({\mathcal {I}}_{RT}\) could be omitted in the formula (41). The partition of unity \(\sum _{z\in {\mathcal {V}}}\varphi _z\equiv 1\) and \(G_h= {\mathcal {I}}_{RT} G_h= \sum _{z \in {\mathcal {V}}} {\mathcal {I}}_{RT} (\varphi _z G_h)\) show that the sum \(Q_p = \sum _{z \in {\mathcal {V}}} Q_{z,h}\in H(\mathop {\textrm{div}}\limits , \Omega )\) of the patch-wise contributions satisfies

This establishes the flux reconstruction \(Q_p\). The efficiency of the flux reconstruction will be based on the following general equivalence.

Lemma 9

(control of \(H(\mathop {\textrm{div}}\limits )\) minimization by residual [10, 31]) Given any vertex \(z\in {\mathcal {V}}\), a piecewise Raviart-Thomas function \(\sigma _z\in RT_{q}^{\textrm{pw}}(\mathcal {T}(z))\) and a piecewise polynomial \(r_z \in {P}_{q}\left( {\mathcal {T}}(z)\right) \) of degree \(q \in {\mathbb {N}}_0\), define the residual

for all \(v\in H^1(\Omega )\). If \(z\in {\mathcal {V}}(\Omega )\) is an interior vertex, then suppose additionally that \(\textrm{Res}_z(1) = 0\). Then

holds for a constant \(C_{\textrm{s}}\) that exclusively depends on the shape-regularity (and is in particular independent of the polynomial degree q).

Proof

The assertion follows from [10, Theorem 7] in \(n=2\) dimensions and [31, Corollaries 3.3, 3.6, and 3.8] in \(n=3\) dimensions.\(\square \)

Remark 1

The patch-wise construction of \(Q_p\) in Sect. 4.2 typically generates local data oscillation \(\textrm{osc}_{k+p}(\varphi _z f,{\mathcal {T}}(z))\) in the error analysis as in the proof of Theorem 3 in Sect. 4.6 below or, e.g., [30, Theorem 3.17]. A straightforward computation \(\textrm{osc}_{k+p}(\varphi _z f,{\mathcal {T}}(z)) \le \textrm{osc}_{k+p-1}(f,{\mathcal {T}}(z))\) apparently leads to a loss of one degree in the data oscillation but Lemma 8 verifies

for any \(p, q\in {\mathbb N}_0\). This allows for efficiency of the data oscillations on the right-hand side of the efficiency estimate [30, Formula (3.42)] and leads to a corresponding refinement in [30, Theorem 3.17].

4.3 Local equivalence of stabilizations

The first improvement to the current HHO literature is the local equivalence of the two stabilizations \({\tilde{s}}_h\) from (37) and \(s_h\) from (2). The authors of this paper could not find any motivation for the alternative stabilization \({\tilde{s}}_h\) in the error analysis of [25, Section 4] and suggest to apply Theorem 4 below to [25, Theorem 4.7] to recover the results therein for the original HHO stabilization \(s_h\). Recall the local stabilization \({\tilde{s}}_T\) in \({\tilde{s}}_h(v_h,v_h) :=\sum _{T \in \mathcal {T}} {\tilde{s}}_T(v_h,v_h)\) from (37) and \(S_{TF}v_h = \Pi _{F,k} \left( v_{\mathcal {T}} + (1-\Pi _{T,k})R v_{h} \right) |_{T}-v_{\mathcal {F}}|_{F}\) in the definition of \(s_h\) from (2).

Theorem 4

(local equivalence of stabilizations) Any \(v_h = (v_{\mathcal T},v_{\mathcal {F}}) \in V_h\) and \(T \in \mathcal T\) satisfy

The constants \(C_10\) and \(C_11\) exclusively depend on the polynomial degree k and the shape regularity of \({\mathcal {T}}\).

Proof

The second inequality in (47) follows directly from a triangle inequality and an inverse estimate. Therefore, the proof focuses on the first inequality in (47). Given \(v_h = (v_\mathcal T,v_{\mathcal {F}}) \in V_h\) and \(T \in \mathcal T\), let \(\varphi _k :=(\Pi _k R v_h - v_\mathcal {T})|_{T} \in P_k(T)\). Since \(S_{TF}v_h = \Pi _{F,k}(R v_{h|T} - v_\mathcal {F}|_F - \varphi _k)\), the triangle inequality \(\Vert \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}})\Vert _{L^2(F)} \le \Vert S_{TF}v_h\Vert _{L^2(F)} + \Vert \varphi _k\Vert _{L^2(F)}\), the discrete trace inequality \(\Vert \varphi _k\Vert _{L^2(F)} \lesssim h_F^{-1/2}\Vert \varphi _k\Vert _{L^2(T)}\), and the shape-regularity \(h_F \approx h_T\) for all \(F \in {\mathcal {F}}(T)\) reveal

Since \(\Pi _0 \varphi _k = 0\) (from the design of \(Rv_h\)), a Poincaré inequality shows

On the one hand, an integration by parts provides

On the other hand, an integration by parts and the definition of R imply

Since \(\Delta \varphi _k \in P_k(T)\), the \(L^2\) projection \(\Pi _k\) on the right-hand side of (50) is redundant. Hence, the combination of (50)–(51) with \(\nabla \varphi _k \cdot \nu _{T|F} \in P_k(F)\) for all \(F \in {\mathcal {F}}(T)\) results in

A Cauchy inequality on the right-hand side of (52), a discrete trace inequality, and \(S_{TF}v_h = \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}|F} - \varphi _k)\) for all \(F \in {\mathcal {F}}(T)\) lead to

Since \({\tilde{s}}_T(v_h,v_h) = \sum _{F \in {\mathcal {F}}(T)} h_F^{-1}\Vert \Pi _{F,k}(R v_{h|T} - v_{\mathcal {F}})\Vert _{L^2(F)}^2 + h_T^{-2}\Vert \varphi _k\Vert _{L^2(T)}^2\), the combination of (48)–(49) with (53) concludes the proof of (47).\(\square \)

4.4 Efficiency of the stabilization

The second improvement to the HHO literature is a quasi-best approximation estimate along the lines of the seminal paper [32, Theorem 4.10]. In combination with Theorem 4, this, in particular, provides the efficiency (38) of the stabilization up to data oscillation.

Theorem 5

(quasi-best approximation up to data oscillation) For any \(p\in {\mathbb {N}}_0\), the solution u to (3) and the discrete solution \(u_h\) to (18) satisfy

The constant \(C_12\) exclusively depends on k, p, and the shape regularity of \({\mathcal {T}}\).

Proof

Given \(k,p \in {\mathbb {N}}_0\), let \({\widetilde{u}} \in V\) solve the Poisson model problem \(-\Delta {\widetilde{u}} = \Pi _{k+p} f\) with the right-hand side \(\Pi _{k+p} f\). Section 4.3 in [32] constructs a stable enriching operator \({\mathcal {J}}: V_h \rightarrow V\) with local bubble-functions from [38]. A modification, where the polynomial degree \(k-1\) in [32, Eq. (4.16)] is replaced by \(k+p\), leads to a right-inverse \({\mathcal {J}}: V_h \rightarrow V\) of the interpolation \(\textrm{I}: V_h \rightarrow V\) with the stability \(|\!|\!|{\mathcal {J}} v_h|\!|\!|^2 \lesssim a_h(v_h,v_h)\) and the additional \(L^2\) orthogonality \({\mathcal {J}} v_h - v_{\mathcal {T}} \perp P_{k+p}({\mathcal {T}})\) for all \(v_h = (v_{\mathcal {T}},v_{\mathcal {F}}) \in V_h\). The extra orthogonality allows for

Consequently, the smoother \({\mathcal {J}}\) leads to a discrete solution \(u_h=(u_\mathcal {T},u_\mathcal {F})\in V_h\) in the modified HHO discretization of [32] as a quasi-best approximation of the above solution \({{\tilde{u}}}\). The point is that \(u_h\in V_h\) coincides with the original HHO solution \(u_h\) from (18). The arguments from the proof of [32, Theorem 4.10] reveal the quasi-best approximation

also for the above modified smoother \({\mathcal {J}}\). This, the triangle inequalities \(|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\le |\!|\!|u - {\tilde{u}}|\!|\!| + |\!|\!|{\tilde{u}} - R u_h|\!|\!|_\textrm{pw}\) and \(|\!|\!|{\tilde{u}} - v_{k+1}|\!|\!|_\textrm{pw}\le |\!|\!|u - {\tilde{u}}|\!|\!| + |\!|\!|u - v_{k+1}|\!|\!|_\textrm{pw}\) for any \(v_{k+1} \in P_{k+1}({\mathcal {T}})\), and the standard perturbation bound \(|\!|\!|u - {\tilde{u}}|\!|\!| \le C_P\textrm{osc}_{k+p}(f,\mathcal {T})\) conclude the proof. \(\square \)

4.5 Stabilization-free efficiency of averaging

The main result of this subsection establishes the stabilization-free efficiency of the nodal averaging technique.

Theorem 6

(averaging is efficient) Let \(u \in V\) resp. \(u_h \in V_h\) solve (3) resp. (18). Then \(R u_h\) and \({\mathcal {A}} R u_h\) satisfy, for any \(p \in {\mathbb {N}}_0\), that

The constant \(C_12\) exclusively depends on k, p, and the shape regularity of the triangulation \({\mathcal {T}}\).

The proof can follow the proof of [25, Theorem 4.7] but additionally utilizes the two significant improvements from Sects. 4.3–4.4 that allow a stabilization-free efficiency in (54).

Proof

Theorem 4.7 in [25] shows \(|\!|\!|(1 - {\mathcal {A}})R u_h|\!|\!|_\textrm{pw}^2 \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}^2 + {\tilde{s}}_h(u_h,u_h)\). This, the equivalence \({\tilde{s}}_h(u_h,u_h) \approx s_h(u_h,u_h)\) of stabilizations (from Theorem 4), and the efficiency \(s_h(u_h,u_h) \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}^2 + \textrm{osc}_{k+p}^2(f,{\mathcal {T}})\) (from Theorem 5) imply the assertion. \(\square \)

Remark 2

(p-robustness) The \(H^1(\Omega )\)-conforming post-processing of \(R u_h\) from [31] provides an efficiency constant independent of the polynomial degree k. The right-hand side of [31, Corollary 4.2] involves the stabilization-related term \(\sum _{F \in \mathcal {F}}h_F^{-1}\Vert \Pi _{F,0}[R u_h]_F\Vert _{L^2(F)}^2\). It remains an open question whether this term can be controlled p-robustly by \(|\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{k}(f,\mathcal {T})\) (with a multiplicative constant independent of the polynomial degree k).

4.6 Proof of Theorem 3

Let \(p \in {\mathbb {N}}_0\) and \(r = 0\) if \(k = 0\) and \(r = k+p\) if \(k \ge 1\) as in (35) be given. Recall the abbreviation \(G_h= \nabla _\textrm{pw}R u_h \in P_k({\mathcal {T}})\) with the discrete solution \(u_h \in V_h\) to (18) and let \(Q_p = \sum _{z \in {\mathcal {V}}} Q_{z,h}\in H(\mathop {\textrm{div}}\limits , \Omega )\) denote the flux reconstruction from Sect. 4.2. The proof establishes (36) in five steps.

Step 1 provides the GUB \(|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{eq},p}({\mathcal {T}})\). This can follow from the paradigm of [1, 2, 30] as outlined below. The choice \({G}:=G_h\) and \(w :={\mathcal {A}} R u_h\) in (8) and a triangle inequality lead to

Since \(\mathop {\textrm{div}}\limits Q_p + \Pi _{r} f = 0\) vanishes in \(\Omega \) by (42), a piecewise Poincaré inequality shows \(|\!|\!|f + \mathop {\textrm{div}}\limits Q_p|\!|\!|_* \le C_P\textrm{osc}_{r}(f,\mathcal {T})\). This, the bound \(|\!|\!|\textrm{div}(Q_p - G_h)|\!|\!|_* \le \Vert Q_p - G_h\Vert \) from the definition of \(|\!|\!|\bullet |\!|\!|_*\), and the previously displayed formula result in the reliability \(|\!|\!|u - R u_h|\!|\!|_\textrm{pw}\le \eta _{\textrm{eq,p}}({\mathcal {T}})\).

Step 2 establishes \(\textrm{osc}_{k-1}(f,{\mathcal {T}}) \lesssim \eta _{\textrm{eq},p}({\mathcal {T}})\). Lemma 8 provides

This, Step 1 and \(\textrm{osc}_{r}(f, \mathcal {T})\lesssim \eta _{\textrm{eq}, p}(\mathcal {T})\) from (34) conclude the proof of Step 2.

Step 3 reveals, for any \(q\in {\mathbb N}_0\), the efficiency of the flux reconstruction

for any polynomial degree \(k \ge 1\). The case \(k=0\) follows in Step 4 below. Recall \(f_z = \Pi _{k+p}(\varphi _z f - G_h\cdot \nabla \varphi _z)\) from (39) for any vertex \(z\in {\mathcal {V}}\) in the construction of \(Q_p\) from Sect. 4.2 and set \(\sigma _z :={\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h)\in RT_{k+p}^{\textrm{pw}}(\mathcal {T}(z))\). The piecewise Raviart-Thomas interpolation \({\mathcal {I}}_{\textrm{RT}}:H^1(\mathcal {T};\mathbb R^n)\rightarrow RT_{k+p}^\textrm{pw}(\mathcal {T})\) satisfies the well-known commuting diagram properties [7, Section 2.5.1]

for any facet \(F\in \mathcal {F}(T)\) of a simplex \(T\in \mathcal {T}\) and the normal trace \(\gamma _T:H^1(T;\mathbb R^n)\rightarrow L^2(\partial T)\) with \(\gamma _T \sigma :=\sigma \cdot \nu _T\) for \(\sigma \in H^1(T;\mathbb R^n)\). This and elementary algebra with the product rule \({\mathop {\textrm{div}}\limits }_{\textrm{pw}}(\varphi _z G_h)=\varphi _z{\mathop {\textrm{div}}\limits }_{\textrm{pw}}G_h+ \nabla \varphi _z\cdot G_h\in P_k({\mathcal {T}}(z))\) imply

Recall the residual \(\textrm{Res}_z(v)\) with \(v\in H^1(\Omega )\) for the vertex \(z\in {\mathcal {V}}\) from (44). The commuting diagram property (57)–(58) establish the identity

This, a piecewise integration by parts, and the property (20) verify \(\textrm{Res}(1) = (G_h,\nabla \varphi _z)_{L^2(\omega (z))} - (f,\varphi _z)_{L^2(\omega (z))}=0\) for any interior vertex \(z \in {\mathcal {V}}(\Omega )\). Hence, Lemma 9 applies for any vertex \(z\in {\mathcal {V}}\) and provides

for the local contributions \(Q_{z,h}\) of \(Q_p = \sum _{z \in {\mathcal {V}}} Q_{z,h}\) from (41). The identity \({\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) = \varphi _z G_h\) for \(p\ge 1\) from \(\varphi _zG_h\in P_{k+1}(\mathcal {T}(z))\subset RT_{k+p}^{\textrm{pw}}(\mathcal {T}(z))\) allows for a k- and p-robust efficiency control of the flux reconstruction error

with a constant \(C_14\) that solely depends on the shape regularity of \({\mathcal {T}}\). This is deemed noteworthy and motivates two different approaches for the bound of the residual \(\textrm{Res}_z(v)\) on the right-hand side of (59) for \(p=0\) and \(p \ge 1\) below.

Step 3.1 provides (56) for \(p=0\). Given any normalized \(v \in H^1_*(\omega (z))\) with \(\Vert \nabla v\Vert _{L^2(\omega (z))}=1\), the product \({\mathcal {I}}_{RT}(\varphi _z G_h) \cdot \nu _F\,v\) vanishes on any boundary facet \(F \in {\mathcal {F}}(\partial \omega (z))\) of the patch \(\omega (z)\). This, the commuting diagram property (57), and (58) with \(\varphi _z{\mathop {\textrm{div}}\limits }_{\textrm{pw}}G_h\in P_k({\mathcal {T}}(z))\) in the definition of the residual \(\textrm{Res}_z(v)\) from (44) verify

The shape regularity of \({\mathcal {T}}\), a Poincaré inequality for interior vertices \(z \in {\mathcal {V}}(\Omega )\), and a Friechrichs inequality for boundary vertices \(z \in {\mathcal {V}}(\partial \Omega )\) provide

Given any facet \(F \in {\mathcal {F}}(z) \cap {\mathcal {F}}(\Omega )\) in the facet spider \(F\in \mathcal {F}(z)\) with facet patch \(\omega (F)\), a trace inequality thus shows

Cauchy inequalities on the right-hand side of (61), the stability of the \(L^2\) projection, \(\Vert \varphi _z\Vert _{L^\infty (\omega _z)} = 1\), and (62)–(63) prove that \({\textrm{Res}}_z(v)\) is controlled by

The efficiency of those residual terms follows from the proof of Theorem 2 in Sect. 3. In combination with (43) and (59), this results in the global efficiency (56) for \(p = 0\) and concludes the proof of Step 3.1.

Step 3.2 provides (56) for \(p\ge 1\). Given any normalized \(v \in H^1_*(\omega (z))\) with \(\Vert \nabla v\Vert _{L^2(\omega (z))}=1\), (58) and the identity \({\mathcal {I}}_{\textrm{RT}}(\varphi _z G_h) = \varphi _z G_h\) in the definition of the residual \(\textrm{Res}_z(v)\) from (44) verify that \(\textrm{Res}_z(v)\) is equal to

This, a piecewise integration by parts, and the product rule \(\nabla (\varphi _z v)=\varphi _z\nabla v+ v\nabla \varphi _z\) reveal

The weak formulation (3) with the test function \(\varphi _z v \in H^1_0(\omega (z))\subset V\) shows \((\nabla u, \nabla (\varphi _z v))_{L^2(\omega (z))} = (f, \varphi _z v)_{L^2(\omega (z))}\). Consequently, (64) implies

This, a Cauchy-Schwarz, and a piecewise Poincaré inequality with the normalization \(\Vert \nabla v\Vert _{L^2(\omega (z))} = 1\) provide

Since \(\Vert \nabla \varphi _z\Vert _{L^\infty (\omega (z))}\approx h^{-1}\), the Leibniz rule, a triangle inequality, and (62) show that

can be bounded by a constant independent of h. Thus, the combination of (43) with (59) and (65) results in the efficiency of \(\Vert Q_p - G_h\Vert \) in (60) with a k- and p-robust constant \(C_{14}\). Since \(\varphi _z \Pi _{k-1} f \in P_k(T) \subset P_{k+p}(T)\), the Pythagoras theorem and \(\Vert \varphi _z\Vert _{L^\infty (\omega (z))} = 1\) show \(\Vert (1 - \Pi _{k+p})(\varphi _z f)\Vert _{L^2(T)} \le \Vert \varphi _z(1 - \Pi _{k-1}) f\Vert _{L^2(T)} \le \Vert (1-\Pi _{k-1})f\Vert _{L^2(T)}\) for all \(T \in {\mathcal {T}}(z)\), whence

This and Lemma 8 implies the efficiency

of the local oscillations from (65) for any \(q\in {\mathbb {N}}_0\) as in Remark 1. This and (60) prove (56) for \(p \ge 1\) and conclude the proof of Step 3.2.

Step 4 affirms the efficiency of the flux reconstruction

for the polynomial degree \(k = 0\) and any \(q\in {\mathbb {N}}_0\). Let \({\widetilde{u}} \in V\) solve the Poisson model problem \(-\Delta {\widetilde{u}} = \Pi _0 f\) with piecewise constant right-hand side \(\Pi _0 f \in P_0({\mathcal {T}})\). A careful inspection reveals that all arguments from Step 3 apply to the case \(k = 0\) for f replaced by \(\Pi _0 f\). This leads to \(\Vert Q_p - G_h\Vert \le C_{15}|\!|\!|{\widetilde{u}} - R u_h|\!|\!|_\textrm{pw}\) with a constant \(C_{15}\) that solely depends on the shape of \({\mathcal {T}}\) (because the data oscillations on the right-hand side of (60) vanish). Therefore, the triangle inequality \(|\!|\!|{\widetilde{u}} - R u_h|\!|\!|_\textrm{pw}\le |\!|\!|u - {\widetilde{u}}|\!|\!| + |\!|\!|u - R u_h|\!|\!|_\textrm{pw}\), the standard bound \(|\!|\!|u - {\widetilde{u}}|\!|\!| \le C_P\textrm{osc}_{0}(f,{\mathcal {T}})\), and \(C_P = 1/\pi < 1\) on simplicial domains lead to (66). This and the efficiency of the data oscillations from Lemma 8 conclude the proof of Step 4. Notice that the constant \({C_{15}}\) is independent of the parameter p.

Step 5 finishes the proof. On the one hand, the reliability \(|\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{k-1}(f,{\mathcal {T}}) \lesssim \eta _{\textrm{eq},p}\) is established in Step 1–2. On the other hand, the efficiency of the averaging operator from Theorem 6, Lemma 8, and the efficiency of the flux reconstruction in (56) for \(k \ge 1\) and in (66) for \(k = 0\) imply the efficiency \(\eta _{\textrm{eq,p}} \lesssim |\!|\!|u - R u_h|\!|\!|_\textrm{pw}+ \textrm{osc}_{q}(f,{\mathcal {T}})\) for any \(q \in {\mathbb {N}}_0\). This concludes the proof of Theorem 3.

This section provides numerical evidence for optimal convergence and a comparison of the stabilization-free GUBs \(\eta _{\textrm{res}}\) and \(\eta _{\textrm{eq,p}}\) from the Sects. 3 and 4 with the original error estimator \(\eta _{\textrm{HHO}}\) from [25, Theorem 4.3] for the HHO method in three 2D benchmarks.

5.1 A posteriori error estimation with explicit constants

All triangulations in this section consist of right-isosceles triangles with the Poincaré constant \(C_P = (\sqrt{2}\pi )^{-1}\) [34]. With the estimates \(C_\textrm{1}\le C_\mathcal {T}\), \(C_{\textrm{H}}\le 1\), and \(C_\textrm{2}\le C_\mathcal {F}\) from Example 1, the residual-based error estimator from (22) reads

The equilibrated GUB from Theorem 3,

depends on the post-processed quantity \(Q_p^\Delta :=Q_p-\nabla _{\textrm{pw}}Ru_h\) with the Raviart-Thomas function \(Q_p\in RT_{k+p}(\mathcal {T})\) of degree \(k+p\) for \(p\in {\mathbb {N}}_0\) from Sect. 4.2. An algorithmic description of the computation of \(Q_p^\Delta \) for arbitrary p follows in Appendix A. Theorem 3 shows that \(\eta _{\textrm{eq}, p}\) is efficient and reliable for all \(p\in {\mathbb {N}}_0\). The residual-based error estimator \(\eta _{\textrm{HHO}}\) from [25, Theorem 4.3] reads

with the operator \(R_{T, F}^k\) from [25, Eq. (2.59)] for the original HHO stabilization [25, Eq. (2.22)] that induces the global stabilization \(s_h\). The constant \(C_{\partial T}=12C_P(C_P+C_{\textrm{tr}})\le 2.0315\) for right-isosceles triangles bounds the trace of \(f-\Pi _0 f\) for \(f\in H^1(T)\) and improves on the estimate \(C_{\partial T}\le (C_P(h_T|\partial T|/|T|)(1+C_P))^2=2.524\) from [25, Subs. 4.1.1].

Lemma 10

(Poincaré-type inequality on trace) Given a simplex \(T \subset {\mathbb {R}}^n\), any \(f\in H^1(T)\) and \(C_{\partial T}:=(n+1)h_T^2|T|^{-1} C_P(C_P + 2C_{\textrm{tr}}/n)\) satisfy

Proof

Abbreviate \({\tilde{\ell }}(F):=(n+1)h_Fh_T^2|T|^{-1}\) for the facet \(F\in \mathcal {F}(T)\) of T and apply Lemma 3 to the singleton triangulation \(\{T\}\) and \(f-\Pi _0f\in H^1(T)\). This and the Poincaré inequality \(\Vert f-\Pi _{0, T}f\Vert _{L^2(T)}\le C_P h_T\Vert \nabla f\Vert _{L^2(T)}\) show

The assertion (67) follows from the observation that \({\tilde{\ell }}(F)\) is maximized on the facet \(F\in \mathcal {F}(T)\) with \(h_F=h_T\).\(\square \)

5.2 Implementation and adaptive algorithm

Our implementation of the HHO method in MATLAB uses nodal bases for the spaces \(P_k(\mathcal {F}), P_k(\mathcal {T})\), and \(P_{k+1}(\mathcal {T})\) and the direct solver mldivide (behind the \-operator) for the discrete system of equations representing (18). For implementation details on the HHO method itself we refer to [25, Appendix B]. The integration of polynomial expressions is carried out exactly. The errors in approximating non-polynomial expressions, such as exact solutions u and source terms f, by polynomials of sufficiently high degree are expected to be very small and are neglected for simplicity.

Algorithm 1 displays the standard adaptive algorithm (AFEM) [16, 22] driven by the refinement indicators, for any triangle \(T\in \mathcal {T}\),

with the modified jump \([\bullet ]_F:=\bullet \times n_F\) along a boundary side \(F\in \mathcal {F}(\partial \Omega )\) and the newest-vertex-bisection (NVB). The sum of this over all triangles is, up to some multiplicative constants, equivalent to \(\eta _\textrm{res}^2({\mathcal {T}})\).

AFEM algorithm

Full size image

5.3 High oscillations on the unit square

This benchmark on the unit square \(\Omega :=(0,1)^2\) considers the Laplace equation \(-\Delta u=f\) with source term f matching the smooth exact solution

Figure 1 displays the energy norm \(|\!|\!|e|\!|\!|_\textrm{pw}\) of the error \(e:=u-Ru_h\) for uniform and adaptive refinement by Algorithm 1 on the left. The smooth solution allows for optimal convergence rates \((k+1)/2\) in the number ndof of degrees of freedom, while the adaptive mesh sequence leads to a lower energy error with respect to ndof. The GUBs \(\eta _{\textrm{res}}, \eta _{\textrm{eq}, p}\), and \(\eta _{\textrm{HHO}}\) are efficient and therefore equivalent to the energy error \(|\!|\!|e|\!|\!|_\textrm{pw}\), see Fig. 1 on the right for \(k=0, 2\).

Figure 2 shows the efficiency indices \(EF(\eta ):=\eta /|\!|\!|e|\!|\!|_\textrm{pw}\) for the residual-based GUBs \(\eta =\eta _{\textrm{res}}, \eta _{\textrm{HHO}}\) and the equilibration-based GUB \(\eta =\eta _{\textrm{eq}, p}\) for \(p=0,1\). Higher values of p for a more expensive postprocessing in \(\eta _{\textrm{eq}, p}\) do not significantly improve on \(\eta _{\textrm{eq}, 1}\).

Convergence history of the energy error \(|\!|\!|e|\!|\!|_{\textrm{pw}}\) (left) and the GUBs \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}, \eta _{\textrm{eq}, 0}\) (right) on the square domain

Full size image

History of the overestimation factor \(EF(\eta )=\eta /|\!|\!|e|\!|\!|_\textrm{pw}\) for the residual-based error estimators \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}\) (left) and \(\eta _{\textrm{eq}, 0}, \eta _{\textrm{eq}, 1}\) (right) on the unit square

Full size image

5.4 Analytical solution for the slit domain

The source term \(f\in L^2(\Omega )\) in the second benchmark on the slit domain \(\Omega :=(0,1)^2\!\setminus \!([0,1)\times \{0\})\) matches the singular solution (in polar coordinates)

The singularity of u at the origin (0, 0) leads to reduced convergence rates 1/4 under uniform refinement, regardless of the polynomial degree k. Figure 3 shows that the adaptive algorithm recovers optimal rates and verifies the equivalence of the GUBs \(\eta _{\textrm{res}}, \eta _{\textrm{eq}, p}\), and \(\eta _{\textrm{HHO}}\) to the energy error \(|\!|\!|e|\!|\!|_\textrm{pw}\).

The efficiency indices in Fig. 4 show a strong overestimaton by \(\eta _{\textrm{res}}\) in the preasymptotic regime (undisplayed) with values \(EF(\eta _{\textrm{res}})>60\). However, asymptotically the quotients \(EF(\eta )=\eta /|\!|\!|e|\!|\!|_\textrm{pw}\) for the two residual-based GUBs \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}\) differ only by a factor 10, while the equilibrated GUBs \(\eta _{\textrm{eq}, p}\) provide the closest values to 1.

Convergence history of the energy error \(|\!|\!|e|\!|\!|_{\textrm{pw}}\) (left) and the GUBs \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}, \eta _{\textrm{eq}, 0}\) (right) on the slit domain

Full size image

History of the overestimation factor \(EF(\eta )=\eta /|\!|\!|e|\!|\!|_\textrm{pw}\) for the residual-based error estimators \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}\) (left) and \(\eta _{\textrm{eq}, 0}, \eta _{\textrm{eq}, 1}\) (right) on the slit domain

Full size image

5.5 Corner singularity in the L-shaped domain

The third benchmark problem is set in the L-shaped domain \(\Omega = (-1, 1)^2\setminus [0, 1)^2\) with constant right-hand side \(f\equiv 1\) with an exact solution \(u \in H^{1+s}(\Omega )\) for all \(0 \le s < 2/3\) [24, Theorem 14.6]. Figure 5 displays the convergence history of the error \(e:=u-Ru_h\) and compares the adaptive scheme, Algorithm 1, driven by the refinement indicators \(\eta _{\textrm{res}}(T)\) from (68) and

for \(T\in \mathcal {T}\) that are induced from the GUB \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}\), and \(\eta _{\textrm{eq}, 0}\). Here, the norm \(|\!|\!|e|\!|\!|_\textrm{pw}\) of the distance e from the discrete solution \(Ru_h\in P_{k+1}(\mathcal {T})\) over \(\mathcal {T}\) to the unknown solution \(u\in H^1(\Omega )\) is approximated by \(|\!|\!|\widehat{u}-Ru_h|\!|\!|_\textrm{pw}\), where \({{\widehat{u}}}\in P_{k+1}({\widehat{\mathcal {T}}})\) is the HHO approximation of u on an adaptive refinement \({\widehat{\mathcal {T}}}\) of \(\mathcal {T}\) with at least \(2|\mathcal {T}|\le |{\widehat{\mathcal {T}}}|\) elements. The three adaptive schemes (Algorithm 1, driven by \(\eta _{\textrm{res}}(T), \eta _{\textrm{HHO}}(T)\), or \(\eta _{\textrm{eq}, 0}(T)\)) recover optimal rates of convergence and lead to similar local refinement of the adaptive mesh sequences as in Fig. 6.

Convergence history plot of the energy error \(|\!|\!|e|\!|\!|_\textrm{pw}\) on the L-shaped domain with uniform and adaptive refinement with AFEM, driven by \(\eta _{\textrm{res}}(T)\), (left) and for AFEM, driven by \(\eta _{\textrm{res}}(T), \eta _{\textrm{HHO}}(T)\), and \(\eta _{\textrm{eq}, 0}(T)\), (right)

Full size image

Adaptive triangulations on the L-shaped domain for \(k=3\) from AFEM, driven by \(\eta _{\textrm{res}}(T)\) (left, \(|\mathcal {T}|=882\)), driven by \(\eta _{\textrm{HHO}}(T)\) (middle, \(|\mathcal {T}|=907\)), and driven by \(\eta _{\textrm{eq}, 0}(T)\) (right, \(|\mathcal {T}|=919\))

Full size image

5.6 Conclusion

The adaptive mesh-refining algorithm recovers optimal convergence rates in all three benchmarks. This holds for AFEM driven by any of the three refinement indicators derived from the GUB \(\eta _{\textrm{res}}, \eta _{\textrm{HHO}}\), and \(\eta _{\textrm{eq}, p}\). The generated mesh sequences from the adaptive schemes, driven by the different estimators, display a very similar concentration of the local mesh-refinement as in Fig. 6. All three benchmarks verify that the considered error estimators are GUB with reliability constant 1, while the post-processing in the equilibrated GUB \(\eta _{\textrm{eq}, p}\) produces minimal overestimation with significant additional computational costs.