# A posteriori error estimates for discontinuous Galerkin methods on polygonal and polyhedral meshes

Andrea Cangiani, Zhaonan Dong, Emmanuil H Georgoulis

► To cite this version:

Andrea Cangiani, Zhaonan Dong, Emmanuil H Georgoulis. A posteriori error estimates for discontinuous Galerkin methods on polygonal and polyhedral meshes. SIAM Journal on Numerical Analysis, 2023, 61 (5), pp.2352-2380. 10.1137/22M1516701. hal-03753221v3

HAL Id: hal-03753221

https://inria.hal.science/hal-03753221v3

Submitted on 12 Jul 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A POSTERIORI ERROR ESTIMATES

## FOR

## DISCONTINUOUS GALERKIN METHODS
ON POLYGONAL AND POLYHEDRAL MESHES

ANDREA CANGIANI*, ZHAONAN DONG†, AND EMMANUIL H. GEORGOULIS‡

**Abstract.** We present a new residual-type energy-norm *a posteriori* error analysis for interior penalty discontinuous Galerkin (dG) methods for linear elliptic problems. The new error bounds are also applicable to dG methods on meshes consisting of elements with very general polygonal/polyhedral shapes. The case of simplicial and/or box-type elements is included in the analysis as a special case. In particular, for the upper bounds, arbitrary number of very small faces are allowed on each polygonal/polyhedral element, as long as certain mild shape regularity assumptions are satisfied. As a corollary, the present analysis generalizes known *a posteriori* error bounds for dG methods, allowing in particular for meshes with arbitrary number of irregular hanging nodes per element. The proof hinges on a new conforming recovery strategy in conjunction with a Helmholtz decomposition. The resulting *a posteriori* error bound involves jumps on the tangential derivatives along elemental faces. Local lower bounds are also proven for a number of practical cases. Numerical experiments are also presented, highlighting the practical value of the derived *a posteriori* error bounds as error estimators.

**Key words.** discontinuous Galerkin; *a posteriori* error bound; polygonal/polyhedral meshes; polytopic elements; irregular hanging nodes.

**AMS subject classifications.** 65N30, 65M60, 65J10

**1. Introduction.** Recent years have witnessed an extensive activity in the development of various Galerkin methods posed on meshes consisting of general polygonal/polyhedral (henceforth, collectively referred to as *polytopic*) elements. A central question arising is the derivation of computable error bounds for such discretisations, so that the extreme geometric flexibility of such meshes can be harnessed.

Residual-type *a posteriori* error bounds for interior penalty dG methods on composite/polytopic meshes appeared in [30, 21]. Also, in the context of virtual element methods, corresponding bounds are proven in [14, 19], while for the weak Galerkin approach, an *a posteriori* error analysis can be found in [39]. In addition, corresponding results for the hybrid high-order method can be found in [24]. All aforementioned results are proven for shape-regular polytopic meshes, under the additional assumption that the diameters of elemental faces are of comparable size to the element diameter. Although the latter may be a reasonable assumption in the context of standard, simplicial meshes, it can be rather restrictive for general polytopic elements. This is because general polygons/polyhedra with more than *d* faces can be simultaneously shape-regular yet containing *small* faces, i.e., faces whose diameter is arbitrarily small compared to the element diameter.

This work aims exactly at rectifying this restrictive state of affairs. We prove new energy-norm *a posteriori* upper error bounds for interior penalty discontinuous

<footnote>* SISSA, International School for Advanced Studies, via Bonomea 265, I-34136 Trieste, Italy (andrea.cangiani@sissa.it).</footnote>

<footnote>†1) Inria, 2 rue Simone If, 75589 Paris, France, 2) CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée, France (zhaonan.dong@inria.fr)</footnote>

<footnote>‡1) The Maxwell Institute for Mathematical Sciences and Department of Mathematics, School of Mathematical and Computer Sciences, Heriot-Watt University, Edinburgh EH14 4AS, United Kingdom (E.Georgoulis@hw.ac.uk), 2) Department of Mathematics, School of Applied Mathematical and Physical Sciences, National Technical University of Athens, Zografou 15780, Greece & 3) IACMFORTH, Greece.</footnote>
Galerkin (dG) methods posed on meshes containing polytopic elements, including with
arbitrary number of small faces, as long as certain mild shape-regularity assumptions
are satisfied. The case of simplicial and/or box-type elements is included in the
analysis as a special case. For accessibility, we restrict the discussion to a model
elliptic problem, noting, nevertheless, that various generalizations are possible with
minor modifications.

As a general principle, residual-based *a posteriori* error analysis of non-conforming
and, in particular, dG methods requires a recovery of the numerical solution into a
related conforming function. The pioneering work of Karakashian & Pascal [35] (see
also [34]) proposed the recovery of the dG solution by a nodal averaging operator for
which a crucial stability result was proven [35, Theorem 2.2]; cf., also [34, Theorem
2.1] for an extension. This construction allowed for the first rigorous *a posteriori* error
analysis of a dG method for elliptic problems. A number of related results followed,
improving various aspects of the theory; for instance, see [3, 33, 18, 46, 1, 8, 32, 23, 37].
A key reason for the aforementioned restrictive assumption that all elemental faces
are of comparable size to the element diameter in existing *a posteriori* error analysis
for polytopic dG methods [30, 21] is exactly the lack of availability of a stability result
corresponding to [35, Theorem 2.2] for polytopic element meshes containing elements
with small faces.

In this work, we crucially avoid the use of averaging operators. Instead, the
proof of the upper error bound hinges on a new recovery into $H^1$-conforming func-
tions, in conjunction with a Helmholtz decomposition. To complete the analysis, we
also require the existence of appropriate auxiliary simplicial meshes on which quasi-
interpolants are defined. This can be verified in practice using simple and efficient
algorithms. We provide two such algorithms, one based on a sub-mesh and one em-
ploying tools from computational geometry and, in particular, constrained Delaunay
triangulations [20, 43]. The resulting *a posteriori* error bound involves also jumps on
the tangential derivatives along elemental faces.

Local lower bounds are also proven for a number of practical cases, indicating
the optimality of the new estimators. The key challenge in the proof of the latter
is, again, the treatment of small/degenerating element faces and the construction of
respective bubble functions. In particular, local lower bounds for the element residuals
are proven allowing for arbitrarily small faces. Lower bounds for the flux residuals
are proven under more restrictive assumptions (see Assumption 4.6 below), allowing,
nevertheless, for arbitrarily small faces.

We note that the case of meshes consisting of simplicial and/or box-type elements
is included in the analysis as a special case. For such ‘classical’ mesh concepts, the
developments presented below provide a new class of *a posteriori* error bounds appli-
cable even to *hp*-version dG methods, allowing in particular for meshes with arbitrary
number of irregular hanging nodes per element.

The remainder of this work is structured as follows. In Section 2, we define the
elliptic model problem, the admissible meshes, and finite element spaces and the inte-
rior penalty dG method on polytopic meshes. We also prove some important technical
results regarding the construction of auxiliary meshes which will be instrumental in
the proof of a posteriori error bounds. In Section 3, we prove the *a posteriori* upper
error bound, using the aforementioned technical developments, while in Section 4 we
provide respective lower bounds for the energy-norm error for a number of practical
cases. Finally, in Section 5, we present some numerical experiments confirming the ro-
bustness and efficiency of the error bound and highlighting
its practical value as an error estimator.
**2. Model problem and numerical method.** For a Lipschitz domain $\omega \subset \mathbb{R}^d$, $d = 1, 2, 3$, we denote by $H^s(\omega)$ the Hilbertian Sobolev space of index $s \ge 0$ of real-valued functions defined on $\omega$, endowed with the seminorm $|\cdot|_{H^s(\omega)}$ and norm $\|\cdot\|_{H^s(\omega)}$. Furthermore, we let $L_p(\omega)$, $p \in [1, \infty]$, be the standard Lebesgue space on $\omega$, equipped with the norm $\|\cdot\|_{L_p(\omega)}$. In the case $p = 2$, we shall simply write $\|\cdot\|_{\omega}$ to denote the $L_2$-norm over $\omega$ and simplify this further to $\|\cdot\|$ when $\omega = \Omega$, the physical domain. Finally, $|\omega|$ denotes the $d$-dimensional Hausdorff measure of $\omega$.

**2.1. Model problem.** Let $\Omega$ be a bounded, simply connected, and open polygonal/polyhedral domain in $\mathbb{R}^d$, $d = 2, 3$. The boundary $\partial\Omega$ of $\Omega$ is split into two disjoint parts $\Gamma_D$ and $\Gamma_N$ with $|\Gamma_D| \ne 0$. For technical reasons, when $d = 3$ and $|\Gamma_N| \ne 0$, the interface between $\Gamma_D$ and $\Gamma_N$ is assumed that is made up of straight planar segments. We consider the linear elliptic problem: find $u \in H^1(\Omega)$, such that

$$ (2.1) \qquad \begin{aligned} -\nabla \cdot (a\nabla u) &= f && \text{in } \Omega, \\ u &= g_D && \text{on } \Gamma_D, \\ a\nabla u \cdot \mathbf{n} &= g_N && \text{on } \Gamma_N, \end{aligned} $$

with known $f \in L_2(\Omega)$, $g_D \in H^{1/2}(\Gamma_D)$ and $g_N \in L_2(\Gamma_N)$ and symmetric diffusion tensor $a \in [L_\infty(\Omega)]^{d \times d}$ such that

$$ (2.2) \qquad \alpha^*|\xi|^2 \ge \xi^\top a(x)\xi \ge \alpha_*|\xi|^2 > 0 \quad \forall \xi \in \mathbb{R}^d, \quad \text{a.e. } x \in \Omega, $$

for some constants $\alpha^*, \alpha_* > 0$. For simplicity of the presentation, we assume that $a$ is piecewise constant, although this is not an essential restriction for the validity of the developments below.

Setting $H_D^1 := \{v \in H^1(\Omega) : v = 0 \text{ on } \Gamma_D\}$, the weak formulation of (2.1) is: find $u \in H^1(\Omega)$, $u = g_D$ on $\Gamma_D$ such that

$$ (2.3) \qquad \int_{\Omega} a \nabla u \cdot \nabla v \, dx = \int_{\Omega} fv \, dx + \int_{\Gamma_N} g_N v \, ds, $$

for all $v \in H_D^1(\Omega)$. The well-posedness is guaranteed by the Lax-Milgram Lemma.

**2.2. Finite element spaces and trace operators.** We consider meshes $\mathcal{T}$ consisting of general polygonal (for $d = 2$) or polyhedral (for $d = 3$) mutually disjoint open elements $K \in \mathcal{T}$, henceforth termed collectively as *polytopic*, with $\cup_{K \in \mathcal{T}} \bar{K} = \Omega$. Given $h_K := \text{diam}(K)$, the diameter of $K \in \mathcal{T}$, we define the mesh-function $\mathbf{h} : \cup_{K \in \mathcal{T}} K \to \mathbb{R}_+$ by $\mathbf{h}|_K = h_K$, $K \in \mathcal{T}$. Further, we let $\Gamma := \cup_{K \in \mathcal{T}} \partial K$ denote the mesh skeleton and set $\Gamma_{\text{int}} := \Gamma \setminus \partial \Omega$. The mesh skeleton $\Gamma$ is decomposed into $(d-1)$-dimensional simplices $F$ denoting the mesh *faces*, shared by at most two elements. These are distinct from elemental *interfaces*, which are defined as the simply-connected components of the intersection between the boundary of an element and either a neighbouring element or $\partial \Omega$. As such, an interface between two elements may consist of more than one face, separated by hanging nodes/edges shared by those two elements only. This includes both ‘classical’ hanging nodes, typically created by local mesh refinement, and non-standard ones separating non-co-planar faces. The latter may be created, for instance, by a mesh agglomeration procedure; we refer to Fig. 1 for an illustration for $d = 2$.

The finite element space $S_{\mathcal{T}}$ with respect to $\mathcal{T}$ is defined by

$$ S_{\mathcal{T}} \equiv S_{\mathcal{T}}^p := \{u \in L_2(\Omega) : u|_K \in \mathcal{P}_p(K), K \in \mathcal{T}\}, $$
<figure><img src="image_2.png" /><figcaption>FIG. 1. Element $K \in \mathcal{T}$ with its face-wise neighbours; hanging nodes are highlighted by a bullet.</figcaption></figure>

for some $p \in \mathbb{N}$ with $\mathcal{P}_p(K)$ denoting the space of $d$-variate polynomials of total degree up to $p$ on $K$. We stress that the local elemental polynomial spaces employed within $S_{\mathcal{T}}$ are defined in the physical coordinate system, i.e., without mapping from a given reference or canonical frame. This approach allows to retain the full local approximation properties of the underlying finite element space. We refer to [13, 12] for a detailed discussion on the benefits and implementation issues resulting from this choice.

Let $K_i$ and $K_j$ be two adjacent elements of $\mathcal{T}$ sharing a face $F \subset \partial K_i \cap \partial K_j \subset \Gamma_{\text{int}}$. For $v$ and $q$ element-wise continuous scalar- and vector-valued functions, respectively, we define the *average* across $F$ by $\{[v]\}|_F := \frac{1}{2}(v|_{F \cap K_i} + v|_{F \cap K_j})$, $\{[q]\}|_F := \frac{1}{2}(q|_{F \cap K_i} + q|_{F \cap K_j})$, respectively, and the *jump* across $F$ by $[v] := v|_{F \cap K_i} - v|_{F \cap K_j}$, $\{[q]\} := q|_{F \cap K_i} - q|_{F \cap K_j}$, using the convention $i > j$ in the element numbering to determine the sign. On a boundary face $F \subset \Gamma_D$, with $F \subset \partial K_i$, $K_i \in \mathcal{T}$, we set $\{[v]\} := v_i$, $\{[q]\} := q_i$, $[v] := v_i$, and $\{[q]\} := q_i$, respectively.

For $v \in S_{\mathcal{T}}$ we denote by $\nabla_h v$ the element-wise gradient; namely, $(\nabla_h v)|_K := \nabla(v|_K)$ for all $K \in \mathcal{T}$. Also, we denote by $\nabla_{\mathcal{T}} v$ the face-wise tangential gradient operator acting on the traces of $v$ on $\Gamma$, noting that $\nabla_{\mathcal{T}} v$ is double-valued on $\Gamma_{\text{int}}$. With a slight abuse of notation, we use the same symbol to denote the tangential gradient of boundary functions such as the Dirichlet datum $g_D$.

**2.3. Mesh assumptions, inverse inequalities, and approximation results.** Each mesh $\mathcal{T}$ is required to conform to the problem data in the following basic way. First, $\mathcal{T}$ must represent exactly the domain, namely $\cup_{K \in \mathcal{T}} K = \Omega$, and be consistent with the subdivision of $\partial\Omega$ into $\Gamma_D$ and $\Gamma_N$. Moreover, we require resolution of multiscale features of the domain, such as complex boundaries and bottlenecks. Note that, in the context of polytopic meshes, such resolution is not intrinsic in that multiscale geometrical features can be represented by relatively ‘large’ elements with ‘small’ faces. Hence we assume that the local mesh size of each mesh $\mathcal{T}$ is comparable to the local finest scale of $\Omega$. It is clear that such saturation-type assumption can always be satisfied, possibly after a finite number of refinements of an original coarse mesh. Further, we require the following general polytopic mesh regularity assumption.

**ASSUMPTION 2.1 (Mesh regularity).** We assume that each mesh $\mathcal{T}$ satisfies the following mesh regularity conditions. For each $K \in \mathcal{T}$ it holds:

(a) $K$ is star-shaped with respect to an inscribed ball of radius $r_K \ge \tau^{-1}h_K$, centred at some point $\mathbf{x}^0 = \mathbf{x}_K^0 \in K$; see Figure 2 (right panel) for an illustration when $d=2$.
<figure><img src="image_3.png" /><figcaption>FIG. 2. Left panel: A polygonal element $K$ with 7 nodes is subdivided into two polygons $K_1$ and $K_2$ star-shaped with respect to the points $\mathbf{x}_1^0$ and $\mathbf{x}_2^0$, respectively. Right panel: a patch of polygonal elements. The element $K$ has 12 nodes and is star-shaped with respect to a ball centred in $\mathbf{x}^0$.</figcaption></figure>

(b) Each face $F \subset \partial K \cap \partial \Omega$ is star-shaped with respect to a $(d-1)$-dimensional ball of radius $r_F \ge \tau^{-1}h_K$.

Here, $\tau > 1$ is a constant independent of the discretization parameters. In what follows, we assume that the centres of the inscribed balls are selected to be chosen so that $\tau$ is as minimal.

**Remark 2.2.** All results below generalize immediately to meshes containing polytopic elements that are finite unions of star-shaped polytopes; see Figure 2 (left panel) for an example. The minimal modifications required in the proofs are detailed in Remark 2.10 below.

Assumption 2.1 allows for very general element shapes, including non-convex polytopes with arbitrary number of degenerating faces i.e. element faces $F \subset \partial K$ with $|F| \ll h_K^{d-1}$. Notable examples of acceptable subdivisions comprise elements with a bounded number of possibly degenerate hanging nodes and elements with ‘many’ faces obtained by agglomeration of very fine triangulations of the problem domain. Note, however, that Assumption 2.1(b) forbids degenerating/shape-irregular elemental faces on the boundary. This restriction can be relaxed to some extent at the expense of introducing unknown/hard-to-estimate constants, as discussed in Remark 3.7 below.

To the best of our knowledge, Assumption 2.1 allows for the most general polytopic meshes for which a posteriori error bounds are proven, for any Galerkin discretization and for any PDE problem. Nevertheless, it is, perhaps inevitably, more restrictive compared to the respective ones required for stability and a priori error analysis of dG methods; see [11], [12, Section 4] and [10], for details. The key advantage of the present setting is that it allows us to be as explicit as possible in the constants involved in the a posteriori error bounds.

We also require the following, standard, local quasi-uniformity assumption.

**ASSUMPTION 2.3 (local quasi-uniformity).** Each mesh $\mathcal{T}$ is locally quasiuniform, i.e., there exists a constant $\rho > 0$ such that $\rho^{-1} \le h_{K_i}/h_{K_j} \le \rho$, whenever $K_i, K_j \in \mathcal{T}$ share a common face.

We now state and prove a few results stemming from the above shape-regularity and local quasiuniformity assumptions. A first, geometrical, consequence is that the number of interface neighbours of each $K \in \mathcal{T}$ is, in fact, bounded.

**LEMMA 2.4.** *Under Assumptions 2.1 and 2.3, the number of interface neighbours*
of each element of $\mathcal{T}$ is uniformly bounded.

**Proof.** Denote by $\omega_K$ the set of elements in $\mathcal{T}$ which are neighbours of $K$. We derive a rough upper-bound on the cardinality $n$ of $\omega_K$ as follows. If $K' \in \omega_K$, then $K' \subset B_{h_K+h_{K'}}(\mathbf{x}_K^0) \subseteq B_{(1+\rho)h_K}(\mathbf{x}_K^0)$ thanks to Assumption 2.3. On the other hand, Assumption 2.1 implies that $K'$ contains the ball $B_{r_{K'}}(\mathbf{x}_{K'}^0)$ with $r_{K'} \ge \tau^{-1}h_{K'}$. Letting $C_2 := \pi$ and $C_3 := 4\pi/3$, we thus have $|K'| \ge C_d r_{K'}^d \ge C_d \tau^{-d}h_{K'}^d \ge C_d \tau^{-d}\rho^{-d}h_K^d$. Therefore, $n(C_d\tau^{-d}\rho^{-d}h_K^d) \le C_d(1+\rho)^d h_K^d$, or $n \le \tau^d\rho^d(1+\rho)^d$, thereby showing that the number of neighbours of $K$ is uniformly bounded as required. $\square$

In particular, the presence of ‘many, small’ faces is allowed if these are grouped into a few interfaces only; see Figure 2 (right panel) for an example. Complex interfaces may be produced by agglomeration procedures used to perform numerical upscaling of complex domains described through very fine triangulations [5, 12, 40] or within adaptive algorithms to align the mesh to solution features and coefficients anisotropies [27, 12, 15]

**LEMMA 2.5.** Let $K \in \mathcal{T}$ satisfying Assumption 2.1. Then, for all $v \in H^1(\Omega)$, we have the trace estimate

$$ (2.4) \qquad \|v\|_{\partial K}^2 \le C_{\text{tr}} \left( \frac{\zeta}{h_K} \|v\|_K^2 + \frac{h_K}{\zeta} \|\nabla v\|_K^2 \right), $$

for any $\zeta > 0$ and $C_{\text{tr}}$ a positive constant only depending on $\tau$ and on $d$.

Also, for each $v \in \mathcal{P}_p(K)$, the inverse estimate

$$ (2.5) \qquad \|v\|_{\partial K}^2 \le \frac{C_{\text{inv}}}{h_K} \|v\|_K^2, $$

holds with $C_{\text{inv}} := \tau(p+1)(p+d)$.

**Proof.** These estimates are special cases of the corresponding ones presented in [10]. The trace estimate (2.4) follows from [10, Lemma 4.7] in conjunction with Assumption 2.1, while (2.5) follows from [10, Lemma 4.4] along with Assumption 2.1. $\square$

**LEMMA 2.6.** Given $K \in \mathcal{T}$ satisfying Assumption 2.1, for each $v \in H^1(K)$, $K \in \mathcal{T}$, we have the bounds

$$ (2.6) \qquad \|v - \Pi_0 v\|_K \le C_{PF} h_K \|\nabla v\|_K $$

and

$$ (2.7) \qquad \|v - \Pi_0 v\|_{\partial K} \le \tilde{C}_{PF} \sqrt{h_K} \|\nabla v\|_K, $$

with $\Pi_0 : L_2(\Omega) \to S_T^0$, the orthogonal $L_2$-projection onto $S_T^0$, the space of element-wise constants; here $C_{PF}$, $\tilde{C}_{PF} > 0$ depend on $d$ and on $\tau$ only.

**Proof.** The key technical difficulty is to show that $C_{PF}$ and $\tilde{C}_{PF}$ are independent of the shape of $K$. Under Assumption 2.1, for $K \in \mathcal{T}$ we can apply the Poincare-Friedrichs inequalities proven in [50, Theorem 3.5] and [47, Proposition 2.10], with explicit dependence on the shape-regularity constant $\tau$ and dimension $d$, yielding (2.6). Then, (2.7) follows using the trace inequality (2.4). $\square$

A crucial technical aspect of the analysis below is the availability of a shape-regular, auxiliary triangulation defined as follows.

**DEFINITION 2.7 (Auxiliary mesh).** Given the sequence of meshes $\mathcal{T}$, we name auxiliary mesh a corresponding sequence of conforming simplicial meshes $\tilde{\mathcal{T}}$ satisfying: for each $\tilde{\mathcal{T}}$ and $T \in \tilde{\mathcal{T}}$
(a) (shape-regularity) the radius $r_T$ of the largest circle inscribed in $T$ is such that
$r_T \ge \hat{\tau}^{-1}h_T$, where $h_T$ denotes the diameter of $T$;

(b) (local mesh-size compatibility) if $K \in \mathcal{T}$ is such that $K \cap T \neq \emptyset$, it holds
$\hat{\rho}^{-1} \le h_T/h_K \le \hat{\rho}$,
with $\hat{\tau}, \hat{\rho} > 1$ constants independent of the discretization parameters.

An immediate consequence of the above definition is that, if $\hat{T}$ is an auxiliary
mesh sequence, then the number of intersections of each $T \in \hat{T}$ with the elements of
the corresponding polytopic mesh $\mathcal{T}$ is uniformly bounded in function of the shape-
regularity and local quasi-uniformity constants of both the polytopic and auxiliary
mesh. The proof of this fact follows along the lines of that of Lemma 2.4.

We note that the evaluation of the error estimator presented below does not re-
quire the construction of auxiliary meshes in practice. As long as their existence can
be assumed, the a posteriori error bound holds. Moreover, we do not expect such
assumption to limit in any possible way the configurations of very general polytopic
mesh sequences allowed by Assumptions 2.1 and 2.3. Rather, the issue is to show that
auxiliary meshes tightly close to the polytopic meshes can be constructed in princi-
ple. To this end, we present two possible algorithms for the construction of auxiliary
meshes which apply to progressively complicated primal mesh configurations. Both
algorithms are easily and cheaply implementable. As such, if desired, the correspond-
ing auxiliary mesh quality parameters may be computed in practice, thus permitting
the explicit evaluation of their impact on the a posteriori error bound.

**2.3.1. Auxiliary sub-mesh.** Assume that the mesh $\mathcal{T}$ is fully shape-regular in the sense that, for each $K \in \mathcal{T}$, each face $F \in \partial K$ satisfies the shape-regularity property which is stated in Assumption 2.1(b) for boundary faces. Then, an auxiliary mesh can be simply constructed by joining $F$ to $x_K^0$, the centre of star-shapedness of $K \in \mathcal{T}$, for each $K \in \mathcal{T}$ and $F \in \partial K$. This approach can be extended to the more general case in which every interface can be replaced by a shape-regular triangulated surface which does not compromise the shape-regularity of the neighbouring elements. For instance, any of the four circular interfaces appearing in Figure 2 (right panel) may be replaced by the segment joining its end-points. This will result in the auxiliary sub-mesh shown in Figure 3 (left panel).

**2.3.2. Constrained Delauney auxiliary mesh.** Auxiliary sub-meshes are in general not obvious to construct and, moreover, employing a sub-mesh is not possible when the element faces are shape-irregular and/or of arbitrarily small size with respect to the elemental size. In this case, it is necessary to consider auxiliary meshes which are not logically sub-meshes of the corresponding polytopic meshes. One possibility, designed to maximise shape-regularity while maintaining the auxiliary mesh as close as possible to the polytopic mesh in terms of local mesh-size, is to exploit the concept of constrained Delaunay triangulations, introduced in [20] for $d = 2$ and generalized to any $d$ in [43]; see also [26, 44]. We recall their definition, limiting ourselves to the case of interest, namely when the constraints are given by the mesh faces laying on the boundary of $\Omega$.

**Definition 2.8 (Constrained Delaunay triangulation).** Let $\mathcal{X} = \{\mathcal{P}, \mathcal{F}\}$ with $\mathcal{P}$ a set of points in $\Omega \subset \mathbb{R}^d$ and $\mathcal{F}$ the set of boundary faces of the mesh $\mathcal{T}$. A constrained Delaunay triangulation (CDT) associated to $\mathcal{X}$ is a triangulation of $\Omega$, which conforms to $\mathcal{F}$, has $\mathcal{P}$ as its set of internal vertices, and satisfies the following constrained Delaunay property: for every $k \in \{1, \dots, d\}$ and every $k$-dimensional simplex $S$ in the triangulation which is not on $\partial\Omega$, there exists a circle $C$ such that:
<figure><img src="image_4.png" /><figcaption>FIG. 3. A mesh with 5 elements (solid lines) and associated auxiliary meshes (dash lines). Left panel: auxiliary sub-mesh with 24 elements. Right panel: auxiliary constrained Delauney mesh with 16 elements. The dual nodes belonging to the set $\mathcal{P}$ of Definition 2.7 are marked with a grey bullet.</figcaption></figure>

(1) the vertices of $S$ are on the boundary of $C$;

(2) if a vertex $P$ of the CDT is in the interior of $C$, then the straight line connecting $P$ to at least one of the vertices of $S$ intersects $\partial\Omega$. (Then, we say that $P$ cannot be seen from one of the vertices of $S$.)

This definition generalizes the concept of Delaunay triangulations in that, if no constraints are given, it would coincide with the definition of Delaunay triangulations. Moreover, as Delaunay triangulations, CDTs maximizes the minimum angle among all triangulations generated by the cloud of points $\mathcal{P}$ and constrained by $\mathcal{F}$. The existence of CDTs is analyzed in [20, 43, 44]: constrained Delaunay triangulations always exists for $d=2$ while for $d=3$ they exist if any ridge formed by $\mathcal{F}$ is strongly Delaunay. A simplex is strongly Delaunay if the circle $C$ of Definition 2.8 does not enclose any other point in $X$. As shown in [43, 44], this condition can always be satisfied, possibly after the insertion of a finite number of regular nodes on non-strongly Delaunay ridges in the skeleton of $\mathcal{F}$. Moreover, once every boundary edge is strongly Delaunay, the restriction of the CDT on each interface is Delaunay.

Given the polytopic mesh $\mathcal{T}$, here we consider the constrained Delaunay triangulation of $\Omega$ with seeds $\mathcal{P} = \{\mathbf{x}_K^0\}_{K \in \mathcal{T}}$, possibly after the modification of $\mathcal{F}$ discussed above. We refer to Figure 3 (right panel) for an illustration.

**Remark 2.9.** We expect the CDTs associated to $\mathcal{T}$ to always satisfy the auxiliary mesh Definition 2.7 owing to their shape-regularity maximisation property and the fact that the seeds in $\mathcal{P}$ are well-distanced by assumption. However, due to the extreme generality of $\mathcal{T}$, proving this fact appears to be challenging and would result into overly-pessimistic estimation of the shape-regularity and quasi-uniformity constants. Specifically, the difficulty comes from the contrasting requirements of shape-regularity and local mesh-size compatibility, due to which an element of the CDT may overlap with elements of $\mathcal{T}$ which are not direct neighbours; see Figure 4 (right panel) for an example. Hence, in the *a posteriori* error analysis below, we have opted for keeping the requirements of Definition 2.7 as an assumption which can be verified economically and sharply in practice. Indeed, contrary to auxiliary sub-meshes, CDTs can always be constructed and their construction is, in fact, simpler in general. If desired, their qualitative parameters may be easily evaluated using well established and efficient algorithms [20, 43, 44].

**Remark 2.10.** In the case of meshes with elements made of finite unions of star-shaped sub-polytopes, auxiliary meshes should be constructed starting from the centres of the sub-polytopes instead. The results below still hold, as long as local quasi-
uniformity is assumed for the sub-mesh comprising the star-shaped sub-polytopes.

**2.3.3. Auxiliary mesh interpolation and inverse estimates.** We recall the following Scott-Zhang-type quasi-interpolation result [42, 41, 48, 17].

LEMMA 2.11 (Quasi-interpolant). Let $\hat{T}$ be a shape-regular simplicial subdivision of $\Omega$ not containing any hanging nodes. Then, there exists a quasi-interpolation operator $I_h: H^1(\Omega) \to H^1(\Omega) \cap S_{\hat{T}}^1$, such that

$$ (2.8) \qquad \|v - I_h v\|_T + h_T \|\nabla(v - I_h v)\|_T \le C_{sz} h_T \|\nabla v\|_{\omega_T}, $$

where $\omega_T$ denotes the patch of elements in $\hat{T}$ with non-empty intersection with $\bar{T}$. The constant $C_{sz} > 0$ depends only on the shape-regularity constant $\hat{\tau}$ of the auxiliary mesh $\hat{T}$. If the function $v$ has nonhomogeneous piecewise linear trace on $\partial\Omega$, we have $I_h v|_F = v|_F$, for all $F \subset \partial T \cap \partial\Omega$, $T \in \hat{T}$.

Moreover, we have

$$ (2.9) \qquad \|v - I_h v\|_{\partial K}^2 \le C_{in} h_K \|\nabla v\|_{\hat{\omega}_K}^2, $$

with $C_{in} > 0$ depending only on the shape-regularity of $T$ and of $\hat{T}$. Here, $\hat{\omega}_K = \cup \{\omega_T : T \cap K \neq \emptyset, T \in \hat{T}\}$.

*Proof.* The proof of (2.8) can be found in [42, 41, 48, 17] for various levels of generality. Noting that (2.9) refers to the trace on the skeleton of of the original polytopic mesh, we apply the trace inequality (2.4) with $\zeta = 1$,

$$ \begin{aligned} \|v - I_h v\|_{\partial K}^2 &\le C_{tr} (h_K^{-1} \|v - I_h v\|_K^2 + h_K \|\nabla(v - I_h v)\|_K^2) \\ &\le C_{tr} \sum_{T \in \hat{T}: T \cap K \neq \emptyset} (h_K^{-1} \|v - I_h v\|_T^2 + h_K \|\nabla(v - I_h v)\|_T^2), \end{aligned} $$

and (2.9) follows from (2.8), depending on the shape-regularity of $T$ through $C_{tr}$ and on the shape-regularity of $\hat{T}$ through (2.8). $\square$

The next polynomial inverse estimation result, relating $L_2$-norms on subsets of the mesh skeleton $\Gamma$ of $T$ with $L_2$-norms over elements of the auxiliary mesh $\hat{T}$ will be important for the analysis below. In this context, for each $T \in \hat{T}$, we consider the set of cut interfaces obtained by the intersection of $\Gamma$ with the simplex $T$, which we characterise as follows:

$$ (2.10) \qquad \Gamma_T := \Gamma \cap T = \Gamma_T^c \cup \Gamma_T^o, $$

with $\Gamma_T^c$ the set of interfaces $\gamma_T \in \Gamma_T \cap \Omega$ such that $\gamma_T \subset \partial K$ with $K$ an element of $T$ whose centre of star-shapedness is a vertex of $T$, and $\Gamma_T^o := \Gamma_T \setminus \Gamma_T^c$. Note that the number of interfaces in $\Gamma_T$ is bounded since the number of intersections of $T$ with the elements $K \in T$ is bounded; however, each such interface may be made of an arbitrary number of cut faces, due to the complexity of the intersecting polytopic elements.

The subdivision in (2.10) reflects increasing levels of difficulty, with $\Gamma_T^o$ collecting complex interfaces for which the proof of the inverse estimate is more challenging. As usual, the proof rests in employing simplices obtained by joining each face composing $\Gamma_T$ with the most appropriate vertex of $T$ and summing up all contributions. When $\Gamma_T^o \neq \emptyset$, such simplices may overlap. In this case, the constant of the resulting inverse estimate depends on the number of such overlaps, and thus reflect the complexity of $\Gamma_T$. The complex interface $\gamma_T$ highlighted in Figure 4 (right panel) provides an example in which this eventuality may occur.
<figure><img src="image_5.png" /><figcaption>FIG. 4. Constrained Delauney auxiliary mesh examples. Different configurations for the set $\Gamma_T$ of (2.10) for the auxiliary element T marked in grey. Left: $\Gamma_T = \Gamma_T^c$. Centre: $\Gamma_T = \Gamma_T^o$. Right: $\Gamma_T^o$ is made of a single interface $\gamma_T$ while the rest of $\Gamma_T$ belongs to $\Gamma_T^c$. The simplices required in the proof of Lemma 2.13 to estimate the contributions from each edge in $\gamma_T$ partially overlap (not shown).</figcaption></figure>

**Remark 2.12.** In the case of sub-mesh auxiliary meshes, we always have $\Gamma_T^o = \emptyset$. Instead, for constrained Delauney auxiliary meshes, $\Gamma_T^o \neq \emptyset$ in general; we refer to Figure 4 for some examples.

**LEMMA 2.13.** Let Assumptions 2.1 and 2.3 hold and let $\hat{T}$ be an auxiliary mesh, related to $T$. Given $q \in \mathbb{N}$, for all $T \in \hat{T}$ and $v \in \mathcal{P}_q(T)$, we have

$$ (2.11) \qquad \|v\|_{\Gamma_T^c}^2 \le \frac{\hat{C}_{\text{inv}}}{h_T} \|v\|_T^2. $$

The constant $\hat{C}_{\text{inv}}$ depends on $d$, $q$, and on shape-regularity and local quasi-uniformity of both $T$ and $\hat{T}$ only. If $\Gamma_T \equiv \Gamma_T^c$ then $\hat{C}_{\text{inv}} = \tau(q+1)(q+d)$. If $\Gamma_T^o \neq \emptyset$, $\hat{C}_{\text{inv}}$ depends also on the number of overlaps $n_o$ required to cover the elements of $\Gamma_T^o$, cf. (2.13). In particular, if $\Gamma_T^o \subset \partial\Omega$, then $n_o \le d$.

**Proof.** We consider $\Gamma_T^c$ and $\Gamma_T^o$ separately, starting with $\Gamma_T^o$. Exploiting the star-shapedness property of $\Gamma_T^o$ with respect to the vertices of $T$, which is inherited from Assumption 2.1, the inverse inequality

$$ \|v\|_{\Gamma_T^c}^2 \le \frac{\tau(q+1)(q+d)}{h_T} \|v\|_T^2, $$

follows in the same way as (2.5).

Considering now the set $\Gamma_T^o$, we observe that each cut interface belonging to $\Gamma_T^o$ may be partitioned into a set of $(d-1)$-dimensional simplices. Indeed, each interface inherits a set of, possibly cut, simplicial faces from $\Gamma$. If a face $F \in \Gamma$ is only partially contained in $\Gamma_T^o$, then $F \cap T$ is still an interval if $d=2$ while it can always be subdivided into four triangles if $d=3$. Let now $\hat{F}$ be one such $(d-1)$-dimensional simplex within $\Gamma_T^o$. We note that, if a simplex $T$ has inradius $r_T$, then, for any given intersecting hyperplane $Z$, there exists a vertex $V$ of $T$ such that $d(V, Z) \ge r$. Otherwise, $S$ must be contained in the region $\{\mathbf{x} \in \mathbb{R}^d : d(\mathbf{x}, Z) < r\}$, in contradiction with the fact that $T$ contains a (closed) ball of radius $r$. It follows that we can always construct a non degenerate simplex $T_{\hat{F}}$ by joining $\hat{F}$ with a vertex $V$ of $T$ such that $d(V, \hat{F}) \ge r$. We thus have, cf. (2.5), the inverse estimate:

$$ (2.12) \qquad \|v\|_{\hat{F}}^2 \le \tau \frac{(q+1)(q+d)}{h_T} \|v\|_{T_{\hat{F}}}^2 \quad \forall v \in \mathcal{P}_q(T). $$
Then, summing up over all $\hat{F} \in \Gamma_T^o$ we conclude

$$ (2.13) \qquad \|v\|_{\Gamma_T^o}^2 \le \frac{\hat{\tau}(q+1)(q+d)}{h_T} \sum_{\hat{F} \in \Gamma_T^o} \|v\|_{T_{\hat{F}}}^2 \le n_o \frac{\hat{\tau}(q+1)(q+d)}{h_T} \|v\|_{T}^2, $$

with $n_o$ the number of overlaps of the simplices $T_{\hat{F}}$, $\hat{F} \in \Gamma_T^o$. In particular, if $\Gamma_T^o$ is only made of boundary interfaces, then $n_o \le d$, as there may be at most $d$ such interfaces, each made of a single $(d-1)$-dimensional simplex. The required estimate now follows by summing up the contributions from $\Gamma_T^c$ and $\Gamma_T^o$. $\square$

**Remark 2.14.** The constant $\hat{C}_{\text{inv}}$ appearing in (2.11) accounts for the complexity of the mesh in terms of topology and shape, quantified by the number $n_o$ of overlap required to cover the mesh skeleton, see the mesh shown in Figure 4 (right) for an illustrative example. In typical practical cases, e.g., meshes stemming from standard algorithms such as Voronoi tessellations, as well as shape-regular adaptively generated meshes, we expect to have $\Gamma_T \equiv \Gamma_T^c$ for the vast majority of auxiliary elements. For instance, for the adaptively refined mesh with multiple hanging-nodes shown Figure 4 (left), $\Gamma_T^o$ is either empty or it contains a single boundary edge, i.e., no overlaps are required, resulting in the ‘ideal’ constant $\hat{C}_{\text{inv}} = \tau(q+1)(q+d)$ for each $T \in \hat{T}$.

**2.4. Discontinuous Galerkin method.** Let $\mathcal{V} := S_T + H^1(\Omega)$. The symmetric interior penalty discontinuous Galerkin method reads: find $u_h \in S_T$ such that

$$ (2.14) \qquad B(u_h, v_h) = \ell(v_h) \quad \text{for all} \quad v_h \in S_T, $$

whereby $B(\cdot, \cdot) : \mathcal{V} \times \mathcal{V} \to \mathbb{R}$ is defined by

$$ (2.15) \qquad \begin{aligned} B(w, v) := & \int_{\Omega} a \nabla_h w \cdot \nabla_h v \, dx + \int_{\Gamma \setminus \Gamma_N} \sigma \llbracket w \rrbracket \llbracket v \rrbracket \, ds, \\ & - \int_{\Gamma \setminus \Gamma_N} (\{\!\{a(\Pi \nabla w) \cdot \mathbf{n}\}\!\}\llbracket v \rrbracket + \{\!\{a(\Pi \nabla v) \cdot \mathbf{n}\}\!\}\llbracket w \rrbracket) \, ds, \end{aligned} $$

for $w, v \in \mathcal{V}$, and $\ell(\cdot) : \mathcal{V} \to \mathbb{R}$ by

$$ \ell(v) := \int_{\Omega} fv \, dx - \int_{\Gamma_D} g_D((a(\Pi \nabla v)) \cdot \mathbf{n} - \sigma v) \, ds + \int_{\Gamma_N} g_N v \, ds, $$

with $\Pi : [L_2(\Omega)]^d \to [S_T]^d$ denoting the orthogonal $L_2$-projection operator onto the (vectorial) finite element space, and $\sigma \in L_\infty(\Gamma \setminus \Gamma_N)$ being the, so-called, *discontinuity-penalization function* given by

$$ (2.16) \qquad \sigma(\mathbf{x}) := \begin{cases} C_\sigma \max_{K \in \{K_i, K_j\}} \left\{ \frac{\bar{a}_K C_{\text{inv}}}{h_K} \right\}, & \mathbf{x} \in F \in \Gamma_{\text{int}}, F \subseteq \partial K_i \cap \partial K_j, \\ C_\sigma \frac{\bar{a}_K C_{\text{inv}}}{h_K}, & \mathbf{x} \in F \in \Gamma_D, F \subset \partial K, \end{cases} $$

with $C_\sigma$ a positive constant and $\bar{a}_K := |\sqrt{a}|_2^2|_K$, $K \in \mathcal{T}$; here $|\cdot|_2$ denotes the natural matrix-$l_2$-norm. The known dependence of the penalty on the local polynomial degree is included in $C_{\text{inv}}$ for brevity; see [12, 10] for details. Note that, using (2.3) and that $\llbracket v \rrbracket = 0$ on $\Gamma \setminus \Gamma_N$ for all $v \in H_D^1$, we have $B(u, v) = \ell(v)$ for all $v \in H_D^1$, with $u \in H^1(\Omega)$ the solution to (2.3).
**Remark 2.15.** To avoid further notational overhead, we opted in exposing the main results for element-wise constant diffusion tensors, i.e., $a \in (S_T^0)^{d \times d}$, and for the classical interior penalty dG method. With minor modifications, the results below can also be extended to more general coefficients. Moreover, we expect that a corresponding analysis to what is presented below holds also for the interior penalty dG variants from [29, 25].

Upon defining the dG-norm by $|||v||| := (||\sqrt{a}\nabla_h v||^2 + ||\sqrt{\sigma}[[v]]||^2_{\Gamma \setminus \Gamma_N})^{1/2}$, we have the following result.

**LEMMA 2.16.** *Under Assumption 2.1, there exists $C_\sigma > 0$, such that*

(2.17) $B(v, v) \ge C_{\text{coer}} |||v|||^2$ and $B(w, v) \le C_{\text{cont}} |||w||| |||v|||$ for all $v \in \mathcal{V}$, respectively, with $C_{\text{coer}}, C_{\text{cont}} > 0$, independent of $h$ of $p$, and $K \in \mathcal{T}$.

We refer to [12, 10] for the proof and the explicit definition of $C_\sigma$. A priori error bounds are also available [12, 10].

**3. A posteriori error analysis.** The following analysis requires Assumptions 2.1 and 2.3 and that an auxiliary mesh according to Definition 2.8 is given.

We decompose the error into two components:

$$e = u - u_h = (u - u_c) + (u_c - u_h) =: e_c + e_d,$$

whereby $u_c \in H^1(\Omega)$ is the recovery of the discrete solution $u_h \in S_T$, defined by

$$(3.1) \qquad B(u_c, v) = B(u_h, v) \quad \forall v \in H_D^1(\Omega),$$

and $u_c = g_D$ on $\Gamma_D$. The existence and uniqueness of $u_c$ is guaranteed by the Lax-Milgram Lemma.

**Remark 3.1.** The construction of $u_c$ is known in the theory of finite element methods and has been used in various contexts, e.g., in [28] for the design of equilibrated flux a posteriori error estimators and in [45] for the analysis of domain decomposition preconditioners. A crucial reason of using this recovery instead of the averaging operator as in [35], is that it is essentially independent of the mesh geometry and topology; this is clearly helpful in the present context of very general polytopic meshes.

**3.1. Bounding the non-conforming error $e_d$.** Inspired by [22, 16], cf. also [6, 9], we decompose the nonconforming error $e_d = u_c - u_h$ further via a Helmholtz decomposition.

**LEMMA 3.2.** *Given that $\Omega$ is simply connected, for any $w \in (L_2(\Omega))^d$, there exists $\xi \in H_D^1(\Omega)$ and $\phi \in [H^1(\Omega)]^{2d-3}$, $d = 2, 3$, such that*

$$(3.2) \qquad a\mathbf{w} = a\nabla\xi + \operatorname{curl}\phi \quad \text{in } \Omega,$$

and $\phi$ can be chosen so that

$$(3.3) \qquad \operatorname{curl}\phi \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_N.$$

Moreover, the following relations hold

$$(3.4) \qquad \|\sqrt{a}\mathbf{w}\|^2 = \|\sqrt{a}\nabla\xi\|^2 + \|a^{-1/2}\operatorname{curl}\phi\|^2,$$

and

$$(3.5) \qquad \|\nabla\phi\| \le C_\Omega \|\operatorname{curl}\phi\|,$$

with a constant $C_\Omega > 0$ only depending on $\Omega$.
Proof. The proof of (3.2) is given in [22, Theorem 3.1] for the $d=2$ case and is extended to $d=3$ in [16]. Since $\nabla\xi$ is orthogonal to curl $\phi$, from the symmetry of the diffusion tensor $a$, the orthogonality (3.4) follows immediately. Finally, the proof of (3.5) can be found in [16]; for $d=2$, we have $C_\Omega = 1$. $\square$

Remark 3.3. The Helmholtz decomposition can be generalised to multiply connected domains [31]. However, concerning the validity in this setting of the relation (3.5), which is fundamental to our analysis, we are only aware of the recent preprint [7]. For this reason, we prefer to limit the current analysis to the simply connected setting leaving possible extensions to future work.

Condition (3.3) imposes a constraint on $\Gamma_N$ for $\phi \in [H^1(\Omega)]^{2d-3}$. Namely, curl $\phi \cdot n_F = 0$ on $F \subset \Gamma_N$, implying that $\phi \in [H^1(\Omega)]^3$ has constant components on each $(d-1)$-dimensional planar subset of $\Gamma_N$.

We apply the Helmholtz decomposition with $\mathbf{w} = \nabla_h e_d$. Hence $\xi \in H_D^1(\Omega)$ and $\phi \in [H^1(\Omega)]^{2d-3}$ are such that $a\nabla_h e_d = a\nabla\xi + \text{curl}\phi$, and we have

$$ (3.6) \qquad \|\sqrt{a}\nabla_h e_d\|^2 = \int_{\Omega} a\nabla_h e_d \cdot \nabla\xi \,dx + \int_{\Omega} \nabla_h e_d \cdot \text{curl}\phi \,dx. $$

Since $u_c \in H^1(\Omega)$ with $u_c|_{\Gamma_D} = g_D$ and $\xi \in H_D^1(\Omega)$, (3.1) implies

$$ \int_{\Omega} a \nabla_h e_d \cdot \nabla \xi \, dx = - \int_{\Gamma_{\text{int}}} \{\!\{ a(\mathbf{\Pi} \nabla \xi) \cdot \mathbf{n} \}\!\} \!\left[\!\left[ u_h \right]\!\right] \, ds - \int_{\Gamma_D} (a(\mathbf{\Pi} \nabla \xi) \cdot \mathbf{n})(u_h - g_D) \, ds. $$

Hence, using the Cauchy-Schwarz inequality, the trace inverse estimate (2.5), the definition of $\sigma$, and the orthogonality (3.4), we have, respectively,

$$ \begin{align} & \int_{\Omega} a \nabla_h e_d \cdot \nabla \xi \, dx \nonumber \\ & \le \|\sigma^{-1/2} \{\!\{ a(\mathbf{\Pi} \nabla \xi) \}\!\} \|_{\Gamma \setminus \Gamma_N} (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \tag{3.7} \\ & \le (C_{\sigma}^{-1} \sum_{K \in \mathcal{T}} \|\mathbf{\Pi}(\sqrt{a}\nabla\xi)\|_K^2)^{1/2} (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \nonumber \\ & \le (C_{\sigma})^{-1/2} \|\sqrt{a}\nabla\xi\| (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \nonumber \\ & \le \|\sqrt{a}\nabla_h e_d\| (\|\sqrt{\sigma}C_{\sigma}^{-1}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}C_{\sigma}^{-1}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2}. \nonumber \end{align} $$

To bound the second term on the right-hand side of (3.6), we first decompose

$$ (3.8) \qquad \int_{\Omega} \nabla_h e_d \cdot \text{curl}\phi \, dx = \int_{\Omega} \nabla_h e_d \cdot \text{curl}(\phi - I_h\phi) \, dx + \int_{\Omega} \nabla_h e_d \cdot \text{curl} I_h\phi \, dx, $$

with $I_h$ the (component-wise if $d=3$) quasi-interpolation operator of Lemma 2.11.

Starting with the first term, observing that $\text{curl}\nabla u_c = 0$ and using the fact that $\phi \in [H^1(\Omega)]^{2d-3}$ satisfying (3.3), implying that $\phi$ is a constant function on each planar section of $\Gamma_N$, and choosing $\phi = I_h\phi$ on $\Gamma_N$. Then we have

$$ \int_{\Omega} \nabla u_c \cdot \text{curl}(\phi - I_h \phi) \, dx = \int_{\Gamma_D} (\phi - I_h \phi) \nabla_T g_D \, ds. $$
Applying integration by parts, observing that $\mathrm{curl}\nabla e_d|_K = 0$, and using $(\phi - I_h\phi)$ is single valued on each face, and $(\phi - I_h\phi)_{\Gamma_N} = 0$, yields

$$ (3.9) \qquad \sum_{K \in \mathcal{T}} \int_K \nabla e_d \cdot \mathrm{curl}(\phi - I_h \phi) \, \mathrm{d}\mathbf{x} \\ = - \int_{\Gamma} (\phi - I_h \phi) \llbracket \nabla_T u_h \rrbracket \, \mathrm{d}s + \int_{\Gamma_D} (\phi - I_h \phi) \nabla_T g_D \, \mathrm{d}s \\ = - \int_{\Gamma_{\mathrm{int}}} (\phi - I_h \phi) \llbracket \nabla_T u_h \rrbracket \, \mathrm{d}s - \int_{\Gamma_D} (\phi - I_h \phi) \nabla_T (u_h - g_D) \, \mathrm{d}s. $$

Further, using (2.9) and, finally, (3.5) and (3.4), the right-hand side of (3.9) can be further estimated from above by

$$ \begin{aligned} & \sum_{K \in \mathcal{T}} \|h_K^{-1/2}(\phi - I_h \phi)\|_{\partial K \setminus \Gamma_N} (\|\sqrt{h}_K \llbracket \nabla_T u_h \rrbracket\|_{\partial K \cap \Gamma_{\mathrm{int}}} + \|\sqrt{h}_K \nabla_T (u_h - g_D)\|_{\partial K \cap \Gamma_D}) \\ & \le (\sum_{K \in \mathcal{T}} C_{\mathrm{tr}} C_I \|\nabla \phi\|_{\tilde{\omega}_K}^2)^{1/2} (\|\sqrt{h} \llbracket \nabla_T u_h \rrbracket\|_{\Gamma_{\mathrm{int}}}^2 + \|\sqrt{h} \nabla_T (u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \\ & \le C_1 \|\nabla \phi\|_{\Omega} (\|\sqrt{h} \llbracket \nabla_T u_h \rrbracket\|_{\Gamma_{\mathrm{int}}}^2 + \|\sqrt{h} \nabla_T (u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \\ & \le C_1 C_{\Omega} \sqrt{\alpha^*} \|\sqrt{a} \nabla_h e_d\| (\|\sqrt{h} \llbracket \nabla_T u_h \rrbracket\|_{\Gamma_{\mathrm{int}}}^2 + \|\sqrt{h} \nabla_T (u_h - g_D)\|_{\Gamma_D}^2)^{1/2}, \end{aligned} $$

for $C_1$ a constant depending only on the shape-regularity constant $\hat{\tau}$ of the auxiliary mesh $\hat{T}$, and on the local quasi-uniformity constants $\rho$ and $\hat{\rho}$.

We now consider the second term in (3.8). Since $I_h\phi \in [H^1(\Omega)]^{2d-3}$, we have $\nabla \cdot \mathrm{curl} I_h\phi = 0$ on $\Omega$ and, hence, $\llbracket \mathrm{curl} I_h\phi \cdot \mathbf{n} \rrbracket = 0$ on $\Gamma_{\mathrm{int}}$. Moreover, given that $I_h\phi$ is constant on each component of $\Gamma_N$, we also have $\mathrm{curl} I_h\phi \cdot \mathbf{n} = 0$ on $\Gamma_N$. Then, integration by parts and working as above gives

$$ \int_{\Omega} \nabla_h e_d \cdot \mathrm{curl} I_h \phi \, \mathrm{d}\mathbf{x} = \sum_{K \in \mathcal{T}} \int_{\partial K} e_d \mathbf{n}_K \cdot \mathrm{curl} I_h \phi \, \mathrm{d}s \\ = - \int_{\Gamma_{\mathrm{int}}} \llbracket u_h \rrbracket (\mathbf{n} \cdot \mathrm{curl} I_h \phi) \, \mathrm{d}s - \int_{\Gamma_D} (u_h - g_D)(\mathbf{n} \cdot \mathrm{curl} I_h \phi) \, \mathrm{d}s. $$

Next, applying to $I_h\phi$ the trace inverse inequality with respect to the auxiliary mesh $\hat{T}$ given in (2.11), we obtain

$$ \begin{aligned} & \int_{\Omega} \nabla_h e_d \cdot \mathrm{curl} I_h \phi \, \mathrm{d}\mathbf{x} \\ & \le (\|\sqrt{\sigma} \llbracket u_h \rrbracket\|_{\Gamma_{\mathrm{int}}}^2 + \|\sqrt{\sigma}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \left(\sum_{T \in \hat{T}} \|\sigma^{-1/2} \mathrm{curl} I_h \phi\|_{T \cap (\Gamma \setminus \Gamma_N)}^2\right)^{1/2} \\ & \le (\|\sqrt{\sigma} \llbracket u_h \rrbracket\|_{\Gamma_{\mathrm{int}}}^2 + \|\sqrt{\sigma}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \left(\sum_{T \in \hat{T}} \hat{C}_{\mathrm{inv}} h_T^{-1} \|\sigma^{-1/2} \mathrm{curl} I_h \phi\|_T^2\right)^{1/2} \\ & \le C_2 (\|\sqrt{\sigma C_{\sigma}^{-1}} \llbracket u_h \rrbracket\|_{\Gamma_{\mathrm{int}}}^2 + \|\sqrt{\sigma C_{\sigma}^{-1}}(u_h - g_D)\|_{\Gamma_D}^2)^{1/2} \|\mathrm{curl} I_h \phi\|, \end{aligned} $$

for $C_2 > 0$ constant depending on $C_{\mathrm{inv}}$, $\hat{C}_{\mathrm{inv}}$, $\hat{\rho}$, and on $\alpha^*$. Next, we use the stability of $I_h$ from Lemma 2.11, together with (3.5) to deduce

$$ (3.11) \qquad \|\mathrm{curl} I_h \phi\|_{\Omega} \le C_c \|\nabla \phi\|_{\Omega} \le C_c C_{\Omega} \sqrt{\alpha_*} \|a^{-1/2} \nabla \phi\| \le C_c C_{\Omega} \sqrt{\alpha_*} \|\sqrt{a} \nabla_h e_d\|. $$
Hence, combining (3.7), (3.9), (3.10), and (3.11), we arrive at the bound

$$ (3.12) \qquad \|\sqrt{a\nabla_h e_d}\|^2 \le C_{nc} \left( \|\sqrt{\sigma C_\sigma^{-1}}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma C_\sigma^{-1}}(u_h - g_D)\|_{\Gamma_D}^2 \right. \\ \qquad \qquad \qquad \qquad \left. + \|\sqrt{h}[\nabla_T u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{h}\nabla_T(u_h - g_D)\|_{\Gamma_D}^2 \right); $$

the constant $C_{nc} > 0$ depends on $C_{inv}$, $\hat{C}_{inv}$, $\hat{\rho}$, $C_\Omega$, and on $\alpha^*$, $\alpha_*$, but is independent from $h$ and the number and measure of the mesh faces.

Remark 3.4. We stress that the mesh-size $h$ in (3.12) is the local element diameter for $\mathcal{T}$, i.e.,independent of the measure of the faces and the number of faces per element. This new bound refines the, now classical, results in [35], by showing that the dG error has, in fact, two sources: the normal flux and the tangential gradient. By applying the inverse inequality on each face $F$, the $L_2$-norm of the tangential jump can be bounded from above by the $L_2$-norm of the jump term itself, thus recovering the bound in [35]. However, such bound would be proportional to $\text{diam}(K)/\text{diam}(F)$, for each face $F \in \Gamma$, which may be severely pessimistic for increasingly small faces $F$.

## 3.2. Bounding the conforming error $e_c$. For $v \in H_D^1$, we have

$$ (3.13) \qquad B(e, v) = \ell(v) - B(u_h, v) = \ell(\eta) - B(u_h, \eta), $$

with $\eta := v - v_h$, for any $v_h \in S_T$. Recalling that $e = e_c + e_d$, since $e_c \in H_D^1$ we can fix $v = e_c$ in (3.13) to further deduce

$$ (3.14) \qquad \|\sqrt{a\nabla e_c}\|^2 = B(e_c, e_c) = (\ell(\eta) - B(u_h, \eta)) - B(e_d, e_c) = \ell(\eta) - B(u_h, \eta), $$

from (3.1). The right-hand side of (3.14) can now be bounded via standard arguments [35]: integration by parts, application of [4, Eq. (3.3)], the observation that $\Pi\nabla u_h = \nabla u_h$, and elementary manipulations yield

$$ \begin{aligned} \ell(\eta) - B(u_h, \eta) = & \int_{\Omega} (f + \nabla_h \cdot (a \nabla_h u_h)) \eta \,dx - \int_{\Gamma_{\text{int}}} \llbracket a \nabla u_h \cdot \mathbf{n} \rrbracket \llbracket \eta \rrbracket \,ds \\ & - \int_{\Gamma_N} (a \nabla u_h \cdot \mathbf{n} - g_N) \eta \,ds - \int_{\Gamma_{\text{int}}} \sigma \llbracket u_h \rrbracket \llbracket \eta \rrbracket \,ds - \int_{\Gamma_D} \sigma(u_h - g_D) \eta \,ds \\ (3.15) \qquad & + \int_{\Gamma_{\text{int}}} \llbracket a(\Pi\nabla\eta) \cdot \mathbf{n} \rrbracket \llbracket u_h \rrbracket \,ds + \int_{\Gamma_D} (a(\Pi\nabla\eta) \cdot \mathbf{n})(u_h - g_D) \,ds =: \sum_{i=1}^7 T_i. \end{aligned} $$

Setting $\eta = e_c - \Pi_0 e_c$ and using (2.6), we have

$$ (3.16) \qquad T_1 \le C_{PF} \alpha_*^{-1/2} \left( \sum_{K \in \mathcal{T}} \|h_K(f + \nabla \cdot (a \nabla u_h))\|_K^2 \right)^{1/2} \|\sqrt{a} \nabla e_c\|. $$

Employing (2.7), along with standard manipulations, we also have

$$ \begin{aligned} T_2 + T_3 \le & \sum_{K \in \mathcal{T}} \sum_{F \subset \partial K \cap \Gamma_{\text{int}}} \|h_K^{-1/2} \eta\|_F \|\sqrt{h}_K \llbracket a \nabla u_h \cdot \mathbf{n} \rrbracket\|_F \\ (3.17) \qquad & + \sum_{K \in \mathcal{T}} \sum_{F \subset \partial K \cap \Gamma_N} \|h_K^{-1/2} \eta\|_F \|\sqrt{h}_K (a \nabla u_h \cdot \mathbf{n} - g_N)\|_F \\ & \le \tilde{C}_{PF} \alpha_*^{-1/2} \|\sqrt{a} \nabla e_c\| \left( \|\sqrt{h} \llbracket a \nabla u_h \cdot \mathbf{n} \rrbracket\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{h} (a \nabla u_h \cdot \mathbf{n} - g_N)\|_{\Gamma_N}^2 \right)^{1/2}. \end{aligned} $$
Similarly, using the definition of $\sigma$ from (2.16), we have

$$ T_4 + T_5 \le (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2)^{1/2} \|\sigma^{1/2}[\eta]\|_{\Gamma \setminus \Gamma_{\text{N}}} \\ (3.18) \qquad \le \left( \sum_{K \in \mathcal{T}} \max_{F \in \partial K \setminus \Gamma_{\text{N}}} \sigma \|\eta\|_{\partial K \setminus \Gamma_{\text{N}}}^2 \right)^{1/2} \left( \|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2 \right)^{1/2} \\ \le \sqrt{C_{\sigma} \alpha^*/\alpha_*} C_{\text{inv}} \tilde{C}_{PF} \rho \|\sqrt{a} \nabla e_c\| (\|\sigma^{1/2}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2)^{1/2}. $$

Now, using the trace inverse estimate (2.5), the stability of the $L_2$-projection operator and that $\nabla\eta|_K = \nabla e_c|_K$, we deduce

$$ T_6 + T_7 \le (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2)^{1/2} \|\sigma^{-1/2}\{a(\mathbf{\Pi}\nabla e_c) \cdot \mathbf{n}\}\|_{\Gamma \setminus \Gamma_{\text{N}}} \\ \le \left( \sum_{K \in \mathcal{T}} \max_{F \in \partial K \setminus \Gamma_{\text{N}}} \sigma^{-1} \|\sqrt{a}(\mathbf{\Pi}(\sqrt{a}\nabla e_c)) \cdot \mathbf{n}\|_F^2 \right)^{1/2} \\ (3.19) \qquad \times (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2)^{1/2} \\ \le C_{\sigma}^{-1/2} \rho \|\sqrt{a} \nabla e_c\| (\|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2)^{1/2}. $$

Hence, by collecting above bounds (3.16), (3.17), (3.18), (3.19) and (3.15), we arrive at the following bound on the conforming error:

$$ \|\sqrt{a} \nabla e_c\| \le C_{\text{co}} \left( \|h(f + \nabla_h \cdot (a \nabla_h u_h))\|^2 + \|\sqrt{\sigma}[u_h]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - g_{\text{D}})\|_{\Gamma_{\text{D}}}^2 \right. \\ \left. (3.20) \qquad + \|\sqrt{h}[a \nabla u_h \cdot \mathbf{n}]\|_{\Gamma_{\text{int}}}^2 + \|\sqrt{h}(a \nabla u_h \cdot \mathbf{n} - g_{\text{N}})\|_{\Gamma_{\text{N}}}^2 \right)^{1/2}, $$

with $C_{\text{co}}$ depending on $C_{\sigma}$, $\rho$, $\tau$, the polynomial degree $p$, $C_{PF}$, and $\tilde{C}_{PF}$, but independent of $h$ and the number and measure of the elemental faces.

We are now ready to present the a posteriori error upper bound.

**THEOREM 3.5** (upper bound). Let $u$ be the solution of (2.1) and let $u_h \in S_{\mathcal{T}}$ be its $dG$ approximation on a polytopic mesh satisfying Assumptions 2.1 and 2.3. Also let an auxiliary mesh according to Definition 2.8 is given. Then, we have the following a posteriori error bound

$$ (3.21) \qquad \|u - u_h\|^2 \le C_{\text{up}} \sum_{K \in \mathcal{T}} (R_K^2 + O_K^2), $$

with the local estimator $R_K^2 = R_{K,E}^2 + R_{K,N}^2 + R_{K,J}^2 + R_{K,T}^2$, and the data oscillation $O_K^2 = O_{K,E}^2 + O_{K,N}^2 + O_{K,J}^2 + O_{K,T}^2$, given by

$$ \begin{align*} R_{K,E} &:= \|h(\Pi f + \nabla \cdot (a \nabla u_h))\|_{K}, \\ R_{K,N} &:= (\|\sqrt{h}[a \nabla u_h \cdot \mathbf{n}]\|_{\partial K \cap \Gamma_{\text{int}}}^2 + \|\sqrt{h}(a \nabla u_h \cdot \mathbf{n} - \bar{g}_{\text{N}})\|_{\partial K \cap \Gamma_{\text{N}}}^2)^{1/2}, \\ R_{K,J} &:= (\|\sqrt{\sigma}[u_h]\|_{\partial K \cap \Gamma_{\text{int}}}^2 + \|\sqrt{\sigma}(u_h - \bar{g}_{\text{D}})\|_{\partial K \cap \Gamma_{\text{D}}}^2)^{1/2}, \\ R_{K,T} &:= (\|\sqrt{h}[\nabla_T u_h]\|_{\partial K \cap \Gamma_{\text{int}}}^2 + \|\sqrt{h}\nabla_T(u_h - \bar{g}_{\text{D}})\|_{\partial K \cap \Gamma_{\text{D}}}^2)^{1/2}, \\ O_{K,E} &:= \|h(f - \Pi f)\|_{K}, \quad O_{K,N} := \|\sqrt{h}(g_{\text{N}} - \bar{g}_{\text{N}})\|_{\partial K \cap \Gamma_{\text{N}}}, \\ O_{K,J} &:= \|\sqrt{\sigma}(g_{\text{D}} - \bar{g}_{\text{D}})\|_{\partial K \cap \Gamma_{\text{D}}}, \quad O_{K,T} := \|\sqrt{h}\nabla_T(g_{\text{D}} - \bar{g}_{\text{D}})\|_{\Gamma_{\text{D}}}, \end{align*} $$

with $C_{\text{up}}$ depending on $C_{\text{co}}$ and $C_{\text{nc}}$ only, but is independent of $h$ and of the number and measure of the elemental faces; here, for any for $K \in \mathcal{T}$, such that $\partial K \cap \Gamma_S \neq \emptyset$
with $S \in \{D, N\}$, we set $\tilde{g}_S|_{\partial K \cap \Gamma_S} \in \mathcal{P}_{p_k}(\partial K \cap \Gamma_S)$ with $g_S$ denoting an approximation of the Dirichlet and Neumann data, respectively.

*Proof.* The proof follows immediately from the bounds (3.12) and (3.20), together with the triangle inequality $\|u - u_h\| \le \|e_c\| + \|e_d\|$. □

*Remark 3.6.* In the above, we followed a known approach in splitting the estimator into a ‘residual part’ and a ‘data oscillation part’, assuming that $f \in L_2(\Omega)$ and for sufficiently smooth boundary data. In this setting the data oscillation error is typically dominated by the residual estimators. However, if the forcing data $f \in H^{-1}(\Omega)$, then data oscillation may dominate the error [38]. It would be an interesting future development to investigate the approach from [38] in the context of discontinuous Galerkin methods.

*Remark 3.7.* Theorem 3.5 has been proven under Assumption 2.1(b) which disallows boundary faces with arbitrarily small size relative to the local mesh size. This assumption is reasonable in as much resolution of the problem domain is required in order for the numerical solution to incorporate the boundary conditions. Nevertheless, this assumption can be relaxed in the case of Dirichlet boundary conditions as follows. Noting that the latter is only required to construct the interpolant of the divergence-free component $\phi$ of the non-conforming error, which is not constrained on the Dirichlet boundary. Thus, the interpolant may be constructed for an extension $\tilde{\phi}$, (e.g., as a Stein-type extension operator) defined on an extended domain $\tilde{\Omega} \supset \Omega$ whose respective mesh $\tilde{T}$ would correspond ‘closely’ to the primal mesh $T$ and is constructed so that it may contain no small boundary faces. The resulting bounds, however, would depend on the, typically unknown, boundedness constant of the extension operator.

**4. Lower bounds.** We now derive lower bounds for the a posteriori error estimator of Theorem 3.5. Of particular interest is the extension to which the efficiency of the estimator can be shown to be independent of the number and of the relative sizes of $(d-1)$-dimensional faces in the mesh. The situation differs for the elemental residual and face jump residuals; for clarity, we deal with them separately.

**4.1. Elemental residual.** Lower bounds for the elemental residual can be derived under no further assumptions on the mesh. The analysis is based on a new element bubble function and some auxiliary results.

**LEMMA 4.1** ([10, Corollary 4.24]). Let $\mathcal{T}$ satisfy the Assumption 2.1. Then, for each $K \in \mathcal{T}$, $p \in \mathbb{N}$ and $v \in \mathcal{P}_p(K)$, the following inverse inequality holds

$$ (4.1) \qquad \|\nabla v\|_K^2 \le C_{\text{inv},K} h_K^{-2} \|v\|_K^2. $$

with $C_{\text{inv},K}$ a positive constant depending only on $d$, $p$ and $\tau$. Note also the trivial inequality $\|\text{curl } v\|_K^2 \le (d-1)\|\nabla v\|_K^2$.

Next, for a generic $d$-dimensional simplex $T$, we denote its barycentric coordinates by $\lambda_T^i$, $d = 0, \dots, d$, and denote by $F_i$, $i = 0, \dots, d$ the corresponding $(d-1)$-dimensional simplicial face of $T$ such that $\lambda_T^i|_{F_i} = 0$. Note that $\|\nabla \lambda_T^i\|_{L_\infty(T)} = d|F_i|/|T|$, since $\nabla \lambda_T^i$ is constant. Importantly, the maximum norm is determined by the distance of the $i$-th vertex from the face $F_i$, but it is *independent of the measure* of face $F_i$, see Figure 5 (left) for an illustration.

Let $K \in \mathcal{T}$ and let $m_K$ be the number of its faces. Given that $K$ is star-shaped by Assumption 2.1, we can construct a non-overlapping subdivision of $K$ into $m_K$ simplicial sub-elements $\tau_j$ by joining the face $F_j$, $j = 1, \dots, m_K$, of $K$ with the centre of the largest ball inscribed in $K$; see Figure 5 (right) for an illustration. Note that
<figure><img src="image_6.png" /><figcaption>FIG. 5. Left: a triangle with one 'small' face. Right: a polytopic element with many 'small' faces.</figcaption></figure>

$h_K \ge \text{diam}(\tau_j) \ge r_K \ge \tau^{-1}h_K$. Moreover, letting $\lambda_j := \lambda_{\tau_j}^i$ with $i$ such that $\lambda_{\tau_j}^i$ is the barycentric coordinate of $\tau_j$ corresponding to the vertex of $\tau_j$ which is internal to $K$, it follows that

$$ (4.2) \qquad h_K^{-1} \le \|\nabla \lambda_j\|_{L_{\infty}(\tau_j)} = d|\partial \tau_j \cap \partial K|/|\tau_j| \le \tau h_K^{-1}. $$

**DEFINITION 4.2** (Element bubble). Let $K \in \mathcal{T}$ and let $m_K$ be the number of its faces. With the above notation, the element bubble function $b_K$ is defined as

$$ (4.3) \qquad b_K|_{\tau_j} = \lambda_j, $$

for $j = 1, \dots, m_K$.

By construction, $b_K$ is a continuous piecewise polynomial function with zero trace and with values in $[0, 1]$ on $K$. Next, we will derive some important properties of the new bubble function (4.3).

**LEMMA 4.3.** For each $K \in \mathcal{T}$ satisfying the Assumption 2.1 and for each $v \in P_p(K)$, we have

$$ (4.4) \qquad \|\nabla(b_K v)\|_K^2 \le 2(\tau^2 + C_{\text{inv},K})h_K^{-2}\|v\|_K^2, $$

with $\tau$ as in Assumption 2.1 and $C_{\text{inv},K}$ is in (4.1), and

$$ (4.5) \qquad \|v\|_K^2 \le C_{b,K} \|b_K^{1/2} v\|_K^2, $$

with $C_{b,K} := [2(p+2)]^{2d \frac{d+1}{(d-1)!}}$.

*Proof.* Using the triangle inequality, the bound (4.2), and the inverse inequality (4.1), we have, respectively,

$$ \begin{aligned} \|\nabla(b_K v)\|_K^2 &\le 2\|(\nabla b_K)v\|_K^2 + 2\|b_K(\nabla v)\|_K^2 \\ &\le 2 \sum_{j=1}^{m_K} \|\nabla \lambda_j\|_{L_{\infty}(\tau_j)}^2 \|v\|_{\tau_j}^2 + 2\|b_K\|_{L_{\infty}(K)}^2 \|\nabla v\|_K^2 \end{aligned} $$

$$ (4.6) \qquad \le 2\tau^2 h_K^{-2} \sum_{i=1}^{m_K} \|v\|_{\tau_i}^2 + 2C_{\text{inv},K}h_K^{-2} \|v\|_K^2 = 2(\tau^2 + C_{\text{inv},K})h_K^{-2} \|v\|_K^2, $$
which is the bound required in (4.4). We now prove the norm equivalence relation (4.5). Recalling the norm equivalence relation for each $v \in \mathcal{P}_p(T)$ on a simplex $T$ from [49, Section 3.6], we have

$$ (4.7) \qquad \|v\|_T^2 \le [2(p+2)]^{2d} \left(\frac{d}{d+1}\right)^{d+1} \frac{1}{(d-1)!} \|(\psi_T)^{1/2}v\|_T^2, $$

with $\psi_T := (d+1)^{d+1}(\prod_{i=0}^d \lambda_T^i)$. Then, by using $\|\lambda_T^i\|_{L_\infty(T)} = 1$, $i = 0, \dots, d$, we deduce

$$ (4.8) \qquad \|v\|_T^2 \le C_{b,K} \|(\lambda_T^i)^{1/2} v\|_T^2, $$

for each $i \in \{0, \dots, d\}$. Hence, the bound (4.5) is proven using the definition of $b_K(x)$ in (4.3):

$$ (4.9) \qquad \|v\|_K^2 = \sum_{j=1}^{m_K} \|v\|_{\tau_j}^2 \le \sum_{j=1}^{m_K} C_{b,K} \|(\lambda_j)^{1/2} v\|_{\tau_j}^2 = C_{b,K} \|b_K^{1/2} v\|_K^2. \quad \square $$

**Remark 4.4.** The new element bubble function on polytopic meshes in the above Lemma is different from the classical element bubble function on simplices. In particular, we note that the important relations (4.4) and (4.5) are independent of the number and measure of the faces of the element $K$.

**THEOREM 4.5** (Elemental residual lower bound). Let $u$ be the solution of (2.1) and let $u_h \in S_T$ be its $dG$ approximation under Assumption 2.1 and 2.3. Then, for each $K \in T$, we have

$$ (4.10) \qquad \|h(\Pi f + \nabla \cdot (a \nabla u_h))\|_K^2 \le 2C_{b,K}(2(\tau^2 + C_{\text{inv},K})(\alpha^*)^2)\|\sqrt{a\nabla e}\|_K^2 + O_{K,E}^2. $$

*Proof.* We fix $v \in H_D^1(\Omega)$ as $v|_K = b_K(\Pi f + \nabla \cdot (a \nabla u_h))$, where $b_K$ is the element bubble function in (4.3), and extended to zero outside $K$. Using relations (3.13), (3.15), (4.4), and $b_K \le 1$, we obtain

$$ \begin{align} \int_K b_K (\Pi f + \nabla \cdot (a \nabla u_h))^2 \, d\mathbf{x} &= \int_K (\Pi f - f) v \, d\mathbf{x} + \int_K a \nabla e \cdot \nabla v \, d\mathbf{x} \nonumber \\ &\le \|\Pi f - f\|_K \|v\|_K + \|\sqrt{a \nabla e}\|_K \|\sqrt{a \nabla v}\|_K \tag{4.11} \\ &\le (\|\Pi f - f\|_K + \sqrt{2(\tau^2 + C_{\text{inv},K}) \alpha^* h_K^{-1} \|\sqrt{a \nabla e}\|_K}) \|(\Pi f + \nabla \cdot (a \nabla u_h))\|_K, \nonumber \end{align} $$

Recalling (4.5), we have $\|\Pi f + \nabla \cdot (a \nabla u_h)\|_K^2 \le C_{b,K} \|b_K^{1/2}(\Pi f + \nabla \cdot (a \nabla u_h))\|_K^2$ from which the result (4.10) already follows. $\square$

**4.2. Flux residuals.** In view of proving the lower bound for the flux residuals, we require the number of faces of each element to be uniformly bounded. Furthermore, in the case $d=3$ we shall assume that each face $F$ is shape-regular. Note that such assumptions still allows for arbitrarily small faces.

**ASSUMPTION 4.6.** The number of faces of every element $K \in T$ is uniformly bounded. For $d=3$ only, for every $F \in \Gamma_{\text{int}}$, the radius $r_F$ of the largest $(d-1)$-dimensional ball inscribed in $F$ satisfies $r_F \ge \tau^{-1}h_F$, with $\tau$ as in Assumption 2.1.

Note that the above assumption does not forbid the size of a mesh face to be arbitrarily smaller than that of the elements it belongs to.

To construct the face bubble function, we consider the standard face bubble functions supported in a pair of simplices contained in the neighbouring elements.
**DEFINITION 4.7 (Face bubble).** Let $K \in \mathcal{T}$ and $F \subset \partial K$ a mesh face satisfying Assumption 4.6. Define $T_F^K \subset K$ to be the simplex having $F$ as a face and opposite vertex the point at distance $h_F$ from $F$ along the segment joining the barycentre of $F$ with the centre of star-shapedness of $K$. The face bubble function $b_F$ is defined on $K$ as the standard bubble function of $T_F^K$, cf. [2, 35], extended by zero to the rest of $K$.

**LEMMA 4.8.** Let $K \in \mathcal{T}$ and $F \subset \partial K$ a mesh face under Assumption 2.1. Let $v \in \mathcal{P}_p(F)$ and denote by $v$ also the constant extension in $T_F^K$ of $v$ in the direction normal to $F$. We have

$$ (4.12) \qquad \|v\|_F^2 \le C_{b,F} \|b_F^{1/2} v\|_F^2, $$

with $C_{b,F} := [2(p+2)]^{2(d-1)(d-2)}$, and

$$ (4.13) \qquad \|\nabla(b_F v)\|_{T_F^K}^2 \le C_{\text{inv},T} h_K h_F^{-2} \|v\|_F^2, $$

with $C_{\text{inv},T}$ the constant of the inverse inequality (4.1) in the case of $d-1$-dimensional simplices. Moreover, we have

$$ (4.14) \qquad \|\text{curl}(b_T v)\|_{T_F^K}^2 \le (d-1)C_{\text{inv},T} h_K h_F^{-2} \|v\|_F^2. $$

*Proof.* The bound (4.12) is given in [35] and (4.13) follows immediately from (4.1), the fact that $b_F \le 1$ and the fact that $\|v\|_{T_F^K}^2 \le h_K \|v\|_F^2$. Finally, the bound (4.14) is a trivial consequence of (4.13) observing once again that $\|\text{curl} v\|_K^K^2 \le (d-1)\|\nabla v\|_K^K$.□

**THEOREM 4.9 (flux residuals lower bound).** Let $u$ be the solution of (2.1) and let $u_h \in S_{\mathcal{T}}$ be its $dG$ approximation under Assumption 2.1, 2.3 and 4.6. Then, for each $K \in \mathcal{T}$, we have

$$ (4.15) \qquad \begin{aligned} & \|h^{1/2} \llbracket a \nabla u_h \cdot \mathbf{n} \rrbracket\|_{\partial K \cap \Gamma_{\text{int}}}^2 \\ & \le 6C_{b,F} \left( 2C_{b,K}(\tau^2 + C_{\text{inv},K})(\alpha^*)^2 \|\sqrt{a\nabla e}\|_{\omega_K}^2 + \sum_{K' \in \omega_K} (1+2C_{b,K})O_{K',E}^2 \right. \\ & \qquad \left. + \alpha^* C_{\text{inv},T} \sum_{F \in \partial K \cap \Gamma_{\text{int}}} h_F^{-1} h_{F\perp} \|\sqrt{a\nabla e}\|_{\omega_F}^2 \right), \end{aligned} $$

and

$$ (4.16) \qquad \|h^{1/2} \llbracket \nabla_T u_h \rrbracket\|_{\Gamma_{\text{int}}}^2 \le 2C_{b,F}(d-1)\alpha_* C_{\text{inv},T} \sum_{F \in \partial K \cap \Gamma_{\text{int}}} h_F^{-1} h_{F\perp} \|\sqrt{a\nabla e}\|_{\omega_F}^2; $$

here, $\omega_K$ is the patch of elements neighbouring $K$, $\omega_F = T_F^K \cup T_F^{K'}$ with $K'$ the element neighbouring $K$ across $F$ and $h_{F\perp} := \max\{h_K, h_{K'}\}$.

*Proof.* In view of proving (4.15), we first consider any $F \in \partial K$ with $F \in \Gamma_{\text{int}}$. Further, we fix $v \in H^1(\omega_F)$ as the constant extension of $\llbracket a\nabla u_h \cdot \mathbf{n} \rrbracket|_F$ in the direction normal to $F$, so that $b_F v \in H_0^1(\omega_F)$, cf. Lemma 4.8. Then, testing the error equation (3.13) with $b_F v$ extended to zero on the whole of $\Omega$, we get

$$ \int_{\omega_F} a \nabla_h e \cdot \nabla (b_F v) \, d\mathbf{x} = \int_{\omega_F} (f - \Pi f)(b_F v) \, d\mathbf{x} + \int_{\omega_F} (\Pi f + \nabla \cdot (a \nabla u_h))(b_F v) \, d\mathbf{x} \\ - \int_F \llbracket a \nabla u_h \cdot \mathbf{n} \rrbracket^2 b_F \, ds. $$
By construction, we have $|\omega_F| \sim h_F^d$ for the patches $\omega_F$. The terms involving norms over $\omega_F$ in Theorem 4.9 reflect and account for the presence of relatively small faces. Indeed, linking the size of $\omega_F$ to that of $F$, instead of that of the element $K$, allows very large ratios $h_K/h_F$.

If, on the other hand, the size of each of the element's face is comparable to that of the element itself, then we may modify the construction of the face bubble of Definition 4.7 by moving the opposite vertex all the way to the centre of star-shapedness of $K$. In such a case, $h_F^{-1}h_{F\perp} \|\sqrt{a}\nabla e\|_{\omega_F} \lesssim \|\sqrt{a}\nabla e\|_{K\cup K'}$. Thus, for meshes with potentially many but regular hanging nodes, the new flux-residuals' lower bounds revert to the classical ones, as encapsulated in the following corollary.

**COROLLARY 4.10.** Under the assumption of Theorem 4.9, if, moreover, for every $K \in \mathcal{T}$ and $F \in \partial K$ it holds $h_F \ge \tau^{-1}h_K$, then

$$ (4.19) \qquad \begin{aligned} & \|h^{1/2} \llbracket a\nabla u_h \cdot \mathbf{n} \rrbracket\|_{\partial K \cap \Gamma_{\text{int}}}^2 \\ & \le 6C_{b,F} \left( (2C_{b,K}(\tau^2 + C_{\text{inv},K})(\alpha^*)^2 + \tau\alpha^*C_{\text{inv},K}) \|\sqrt{a}\nabla e\|_{\omega_K}^2 \right. \\ & \qquad \left. + \sum_{K' \in \omega_K} (1 + 2C_{b,K})O_{K',E}^2 \right), \end{aligned} $$

and

$$ (4.20) \qquad \|h^{1/2} \llbracket \nabla_T u_h \rrbracket\|_{\Gamma_{\text{int}}}^2 \le 2C_{b,F}(d-1)\tau\alpha_* C_{\text{inv},K} \|\sqrt{a}\nabla e\|_{\omega_K}^2. $$

**Remark 4.11.** Corollary 4.10 holds in the setting of fully shape-regular meshes allowing for the sub-mesh auxiliary mesh construction, cf. Section 2.3.1. Hence, Corollary 4.10 together with Theorem 4.5 and Theorem 3.5 establishes the reliability and efficiency of the classical residual error estimator for general fully shape-regular polytopic meshes with multiple hanging nodes. The analysis, also in this case, differs from the classical one in that the finite element space used to control the non-conforming error, being based on the auxiliary mesh, is not a subspace of the discrete solution space $S_T$. Moreover, the resulting error bound is as explicit as the classical bound because, in this case, the auxiliary mesh quality is fully controlled by that of the polytopic mesh. The element bubble construction is also new and accounts for the polytopic nature of the mesh.

On the other hand, controlling the flux residuals in the extreme case of possibly unbounded non shape-regular interfaces requires further new ideas. Whenever it is possible to construct a face bubble function $b_F$ such that the following bound

$$ \|\nabla(b_F v)\|_{T_F^K}^2 \le C_{\text{inv},T} h_K^{-1} \|v\|_F^2, $$

holds true, the lower bounds of the flux residuals (4.15) and (4.16) will be independent of $h_F^{-1}h_{F\perp}$. An alternative approach could be to consider bubble functions constructed on a neighbouring set of structured elements. Then, the bubble functions $b_F$ will be independent of the individual face size and the number of elements. For instance, this is the approach used in [36] to derive a lower bound of the flux residual of the FEM employing structured anisotropic triangular meshes. However, the construction of such face bubble functions for the general-shaped polytopic meshes considered in this work is highly non-trivial.

**5. Numerical experiments.** We present two numerical examples testing the new a posteriori error estimator. With the first example we test the impact of polygonal elements with a large number of small faces on the effectivity index. With the
<figure><img src="image_7.png" /><figcaption>FIG. 6. Example 1. Sample meshes with 114 (left) and 498 (right) polygonal elements obtained by agglomeration of a fine triangular mesh made of one million elements.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7. Example 1. Convergence history of the error and error estimator of Equation 3.21 for $p = 1, 2, 3, 4$ (left) and respective effectivity measured as the ratio estimator over error (right).</figcaption></figure>

second, we test the performance of the estimator within a non-standard adaptive algorithm. In all cases, we set $C_{\sigma} = 10$.

**5.1. Example 1.** We construct a sequence of polygonal meshes containing 114, 498, 2063, 8912, and 32768 elements obtained by successive agglomeration of a very fine triangular background mesh made of $10^6$ elements. Each of the polygonal elements contains at least 50 edges, see Figure 6 for an illustration.

We consider the problem (2.1), with $a = I_{2 \times 2}$ and $\Omega := (-1, 1)^2$. The Dirichlet boundary conditions and the source term $f$ are determined by the exact solution $u = \sin(\pi x) \sin(\pi y)$. Numerical results for $p = 1, 2, 3, 4$ are displayed in Figure 7. The observed convergence rate of both the error and the estimator is $O(DoFs^{-p/2})$, i.e., optimal in terms of the total number of degrees of freedom $DoFs$. Moreover, the effectivity index is bounded between 1.2 and 2.6, hence showing that efficiency is not affected by the complexity of the element shapes. This numerical observation reflects that Assumption 4.6 may not be necessary.

Next, we compare the percentage contribution of the different components to the total estimator. Setting $R_X := (\sum_{K \in \mathcal{T}} R_{K,X}^2)^{1/2}$ for $X \in \{E, N, J, T\}$, in Table 1 we provide the percentage of total element residual $R_E$, total jump residual $R_J$, total
jump of the normal flux residual $R_N$, and total jump of the tangential flux residual $R_T$ for $p = 1, 2$. For the coarse meshes of Figure 7, the element residual dominates the total estimator, followed by $R_J$. For finer meshes, we observe significant contribution by $R_E$, $R_J$, and combined $R_T$ and $R_N$. Further, to highlight the importance of

<table><caption>TABLE 1</caption><thead><tr><th rowspan="2"># elem</th><th colspan="4">p = 1</th><th colspan="4">p = 2</th></tr><tr><th>R<sub>E</sub></th><th>R<sub>J</sub></th><th>R<sub>N</sub></th><th>R<sub>T</sub></th><th>R<sub>E</sub></th><th>R<sub>J</sub></th><th>R<sub>N</sub></th><th>R<sub>T</sub></th></tr></thead><tbody><tr><td>114</td><td>64%</td><td>18%</td><td>9%</td><td>9%</td><td>56%</td><td>26%</td><td>9%</td><td>9%</td></tr><tr><td>498</td><td>59%</td><td>23%</td><td>9%</td><td>9%</td><td>45%</td><td>36%</td><td>9%</td><td>10%</td></tr><tr><td>2063</td><td>36%</td><td>37%</td><td>13%</td><td>14%</td><td>45%</td><td>33%</td><td>10%</td><td>12%</td></tr><tr><td>8912</td><td>32%</td><td>37%</td><td>14%</td><td>17%</td><td>43%</td><td>31%</td><td>11%</td><td>15%</td></tr></tbody></table>

Example 1. Estimator's components percentage contributions to the dG error for $p = 1, 2$: element residual $R_E$, jump residual $R_J$, jump of the normal flux residual $R_N$, and jump of the tangential flux residual $R_T$.

the presence of $R_T$ in the new estimator presented in this work, we compare it to the a posteriori error estimator that can be derived using standard techniques from [35], which does not contain the tangential flux jump $R_T$. As already mentioned in Remark 3.4, it is immediate to bound $R_{K,T}$ from $R_{K,J}$ through an inverse estimate on each face $F$, giving $R_{K,T} \lesssim \frac{h_K}{\rho_F} R_{K,J}$, with $\rho_F$ denoting the inscribed radius of the face $F$. Clearly, for small faces, we have $\rho_F \le h_F \ll h_K$, showcasing that $R_{K,T}$ is theoretically sharper than $R_{J,T}$. In Figure 8, we present the error, the estimator from (3.21) and the classic estimator whereby $R_{K,T}$ is replaced by $\frac{h_K}{\rho_F} R_{K,J}$, for $p = 1, 2$. The superiority of the estimator presented in this work is evident for coarse meshes with large ratio $\frac{h_K}{h_F}$. We note that the jump terms account for more than 80% of the classical error estimator, thus indicating that the term $R_{J,T}$ is indeed responsible for the relative over-estimation of the error. This confirms the theoretical intuition and showcases the practicality of the estimator proven in this work.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Example 1. Convergence history of the error, error estimator of Equation 3.21 and the classical estimator for $p = 1, 2$.</figcaption></figure>

**5.2. Example 2.** We test a new adaptive algorithm driven by the error estimator from Section 3. Starting from a relatively coarse simplicial mesh, we use the estimator (3.21) to mark simplicial elements for refinement through a bulk-chasing criterion
(also known as Dörfler marking), and also mark pairs of elements for agglomeration based on the size of the jump residual terms on elemental interfaces. Refinement of simplicial elements is performed via a newest vertex bisection algorithm. In the agglomeration step, general, polygonal meshes will be generated. In successive iterations, polygonal meshes which are marked for refinement are subdivided into either a finer polygonal mesh or a simplicial mesh, depending on their level of agglomeration. For simplicity, we do not consider the data oscillation terms. The adaptive algorithm can thus be described as:

SOLVE → ESTIMATE → MARK → REFINE/AGGLOMERATE.

We consider the problem (2.1) with $a = I_{2 \times 2}$ on $\Omega := (-1, 1)^2 \setminus (0, 1) \times (-1, 0)$. The Dirichlet boundary conditions and the source term $f$ are determined by the exact solution

$$u = r^{2/3} \sin(2\psi/3) + \exp(-1000((x-0.5)^2 + (y-0.25)^2)) + \exp(-1000((x-0.5)^2 + (y-0.75)^2)),$$

which has a point singularity at the origin. We test the adaptive dG algorithm described above with $p = 1, 2, 3$, with Dörfler's marking strategy 25% for refinement and the maximum marking strategy 5% for agglomeration. We point out that the agglomeration step is driven by the jump terms $R_{F,N}$, $R_{F,J}$, and $R_{F,T}$ for all faces $F \subseteq \partial K_i \cap \partial K_j$ on the meshes interface between element $K_i$ and $K_j$.

The performance of the proposed adaptive algorithm is showcased in Figure 9. The convergence rates of both error and estimator are optimal in terms of the total number of degrees of freedom (DoFs), which is $O(DoFs^{-p/2})$. Further, on coarse mesh levels, the mesh agglomeration dominates the mesh refinement. Consequently, the number of DoFs is initially reduced by the adaptive algorithm while the error is also reduced. Another important observation is that the coarse mesh level's effectivity index still seems quite reasonable, namely 2 to 3 times greater than the asymptotic value. This is in spite of the presence of some very large polygonal elements next to small shape-regular triangles. These can be seen, for example, in the sequence of meshes produced by the adaptive algorithm with $p = 3$ shown in Figure 10.

Finally, we perform a comparison between the adaptive algorithm with and without agglomeration for $p = 2, 3$, using the same marking parameter and stopping criterion. The results are presented in Figure 11. Clearly, agglomeration helps reducing DoFs almost without influencing the accuracy on coarse meshes. As the meshes are refined, the advantage of polygonal elements is gradually reduced, as expected.

REFERENCES

[1] M. AINSWORTH, A posteriori error estimation for discontinuous Galerkin finite element approximation, SIAM J. Numer. Anal., 45 (2007), pp. 1777–1798.

[2] M. AINSWORTH AND J. T. ODEN, A posteriori error estimation in finite element analysis, Pure and Applied Mathematics (New York), Wiley-Interscience [John Wiley & Sons], New York, 2000.

[3] M. AINSWORTH AND R. RANKIN, Fully computable error bounds for discontinuous Galerkin finite element approximations on meshes with an arbitrary number of levels of hanging nodes, SIAM J. Numer. Anal., 47 (2010), pp. 4112–4141.

[4] D. N. ARNOLD, F. BREZZI, B. COCKBURN, AND L. D. MARINI, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal., pp. 1749–1779.

[5] F. BASSI, L. BOTTI, A. COLOMBO, AND S. REBAY, Agglomeration based discontinuous Galerkin discretization of the Euler and Navier-Stokes equations, Comput. & Fluids, 61 (2012), pp. 77–85.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Error, estimator (top) and effectivity index (bottom) of the new adaptive polytopic dG method with $p = 1, 2, 3$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. The sequence of meshes generated by the adaptive polytopic dG algorithm with $p = 3$.</figcaption></figure>

[6] R. BECKER, P. HANSBO, AND M. G. LARSON, *Energy norm a posteriori error estimation for discontinuous Galerkin methods*, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 723–733.

[7] F. BERTRAND, C. CARSTENSEN, B. GRÄSSLE, AND N. T. TRAN, *Stabilization-free hho a posteriori error control*, arXiv preprint arXiv:2207.01038, (2022).

[8] A. BONITO AND R. H. NOCHETTO, *Quasi-optimal convergence rate of an adaptive discontinuous Galerkin method*, SIAM J. Numer. Anal., 48 (2010), pp. 734–771.

[9] Z. CAI, X. YE, AND S. ZHANG, *Discontinuous Galerkin finite element methods for interface problems: a priori and a posteriori error estimations*, SIAM J. Numer. Anal., 49 (2011), pp. 1761–1787.
<figure><img src="image_12.png" /><figcaption>FIG. 11. *Example 2.* Error and estimator of the new adaptive polytopic dG method and adaptive dG method without agglomeration for $p = 2, 3$.</figcaption></figure>

[10] A. CANGIANI, Z. DONG, AND E. GEORGOULIS, hp-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements, Math. Comp., 91 (2022), pp. 1–35.

[11] A. CANGIANI, Z. DONG, AND E. H. GEORGOULIS, hp-version space-time discontinuous galerkin methods for parabolic problems on prismatic meshes, SIAM J. Sci. Comput., 39 (2017), pp. A1251–A1279.

[12] A. CANGIANI, Z. DONG, E. H. GEORGOULIS, AND P. HOUSTON, hp-Version Discontinuous Galerkin Methods on Polygonal and Polyhedral Meshes, Springer, 2017.

[13] A. CANGIANI, E. H. GEORGOULIS, AND P. HOUSTON, hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci., pp. 2009–2041.

[14] A. CANGIANI, E. H. GEORGOULIS, T. PRYER, AND O. J. SUTTON, A posteriori error estimates for the virtual element method, Numer. Math., 137 (2017), pp. 857–893.

[15] A. CANGIANI, E. H. GEORGOULIS, AND O. J. SUTTON, Adaptive non-hierarchical Galerkin methods for parabolic problems with application to moving mesh and virtual element methods, Math. Models Methods Appl. Sci., 31 (2021), pp. 711–751.

[16] C. CARSTENSEN, S. BARTELS, AND S. JANSCHE, A posteriori error estimates for nonconforming finite element methods, Numer. Math., 92 (2002), pp. 233–256.

[17] C. CARSTENSEN AND S. A. FUNKEN, Constants in Clément-interpolation error and residual based a posteriori error estimates in finite element methods, East-West J. Numer. Math., 8 (2000), pp. 153–175.

[18] C. CARSTENSEN, T. GUDI, AND M. JENSEN, A unifying theory of a posteriori error control for discontinuous Galerkin FEM, Numer. Math., 112 (2009), pp. 363–379.

[19] C. CARSTENSEN, R. KHOT, AND A. K. PANI, A priori and a posteriori error analysis of the lowest-order ncvm for second-order linear indefinite elliptic problems, Numer. Math., 151 (2022), pp. 551–600.

[20] L. P. CHEW, Constrained Delaunay triangulations, vol. 4, 1989, pp. 97–108. Computational geometry (Waterloo, ON, 1987).

[21] J. CUI, F. GAO, Z. SUN, AND P. ZHU, A posteriori error estimate for discontinuous Galerkin finite element method on polytopal mesh, Numer. Methods Partial Differential Equations, 36 (2020), pp. 601–616.

[22] E. DARI, R. DURAN, C. PADRA, AND V. VAMPA, A posteriori error estimators for nonconforming finite element methods, ESAIM Math. Model. Numer. Anal., 30 (1996), pp. 385–400.

[23] A. DEMLOW AND E. H. GEORGOULIS, Pointwise a posteriori error control for discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal., 50 (2012), pp. 2159–2181.

[24] D. A. DI PIETRO AND R. SPECOGNA, An a posteriori-driven adaptive mixed high-order method with application to electrostatics, J. Comput. Phys., 326 (2016), pp. 35–55.

[25] Z. DONG AND E. H. GEORGOULIS, Robust interior penalty discontinuous Galerkin methods, J. Sci. Comput., 92 (2022), pp. Paper No. 57, 23.

[26] H. EDELSBRUNNER, Geometry and topology for mesh generation, Cambridge Monographs on Applied and Computational Mathematics, Cambridge University Press, 1 ed., 2001.

[27] S. R. ELIAS, G. D. STUBLEY, AND G. D. RAITHBY, An adaptive agglomeration method for additive correction multigrid, Int. J. Numer. Methods Engrg., 40 (1997), pp. 887–903.

[28] A. ERN AND M. VOHRALÍK, Polynomial-degree-robust a posteriori estimates in a unified setting
for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations, SIAM J. Numer. Anal., 53 (2015), pp. 1058-1081.

[29] E. H. GEORGOULIS AND A. LASIS, *A note on the design of hp-version interior penalty discontinuous Galerkin finite element methods for degenerate problems*, IMA J. Numer. Anal., 26 (2006), pp. 381-390.

[30] S. GIANI AND P. HOUSTON, *hp-adaptive composite discontinuous Galerkin methods for elliptic problems on complicated domains*, Numer. Methods Partial Differential Equations, 30 (2014), pp. 1342-1367.

[31] V. GIRAULT AND P. A. RAVIART, *Finite element methods for Navier-Stokes equations*, vol. 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1986. Theory and algorithms.

[32] R. H. W. HOPPE, G. KANSCHAT, AND T. WARBURTON, *Convergence analysis of an adaptive interior penalty discontinuous Galerkin method*, SIAM J. Numer. Anal., 47 (2008/09), pp. 534-550.

[33] P. HOUSTON, D. SCHÖTZAU, AND T. P. WIHLER, *Energy norm a posteriori error estimation of hp-adaptive discontinuous Galerkin methods for elliptic problems*, Math. Models Methods Appl. Sci., 17 (2007), pp. 33-62.

[34] O. A. KARAKASHIAN AND F. PASCAL, *Convergence of adaptive discontinuous Galerkin approximations of second-order elliptic problems*, SIAM J. Numer. Anal., pp. 641-665 (electronic).

[35] ——, *A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems*, SIAM J. Numer. Anal., 41 (2003), pp. 2374-2399.

[36] N. KOPTEVA, *Lower a posteriori error estimates on anisotropic meshes*, Numer. Math., 146 (2020), pp. 159-179.

[37] C. KREUZER AND E. H. GEORGOULIS, *Convergence of adaptive discontinuous Galerkin methods*, Math. Comp., 87 (2018), pp. 2611-2640.

[38] C. KREUZER AND A. VEESEL, *Oscillation in a posteriori error estimation*, Numer. Math., 148 (2021), pp. 43-78.

[39] H. LI, L. MU, AND X. YE, *A posteriori error estimates for the weak Galerkin finite element methods on polytopal meshes*, Commun. Comput. Phys., 26 (2019), pp. 558-578.

[40] M. M. CORTI, P. ANTONIETTI, L. DEDÉ, AND A. QUARTERONI, *Numerical modelling of the brain poromechanics by high-order discontinuous galerkin methods*, arXiv:2210.02272, (2022).

[41] J. M. MELENK, *hp-interpolation of nonsmooth functions and an application to hp-a posteriori error estimation*, SIAM J. Numer. Anal., 43 (2005), pp. 127-155.

[42] L. R. SCOTT AND S. ZHANG, *Finite element interpolation of nonsmooth functions satisfying boundary conditions*, Math. Comp., 54 (1990), pp. 483-493.

[43] J. R. SHEWCHUK, *A condition guaranteeing the existence of higher-dimensional constrained delaunay triangulations.*, in Proceedings of the Fourteenth Annual Symposium on Computational Geometry., Association for Computing Machinery, New York, 1998, pp. 76-85.

[44] ——, *General-dimensional constrained Delaunay and constrained regular triangulations*. I. Combinatorial properties, Discrete Comput. Geom., 39 (2008), pp. 580-637.

[45] I. SMEARS, *Nonoverlapping domain decomposition preconditioners for discontinuous Galerkin approximations of Hamilton-Jacobi-Bellman equations*, J. Sci. Comput., 74 (2018), pp. 145-174.

[46] S. K. TOMAR AND S. I. REPIN, *Efficient computable error bounds for discontinuous Galerkin approximations of elliptic problems*, J. Comput. Appl. Math., 226 (2009), pp. 358-369.

[47] A. VEESEL AND R. VERFÜRTH, *Poincaré constants for finite element stars*, IMA J. Numer. Anal., 32 (2012), pp. 30-47.

[48] R. VERFÜRTH, *Error estimates for some quasi-interpolation operators*, M2AN Math. Model. Numer. Anal., 33 (1999), pp. 695-713.

[49] ——, *A posteriori error estimation techniques for finite element methods*, Numerical Mathematics and Scientific Computation, Oxford University Press, Oxford, 2013.

[50] W. ZHENG AND H. QI, *On Friedrichs-Poincaré-type inequalities*, J. Math. Anal. Appl., 304 (2005), pp. 542-551.