remarkRemark hypothesisHypothesis claimClaim _A posteriori_ error estimates
for dGA. CANGIANI, Z. DONG, E.H. GEORGOULIS

#  _A posteriori_ error estimates
for
discontinuous Galerkin methods
on polygonal and polyhedral meshes

ANDREA CANGIANI SISSA, International School for Advanced Studies, via Bonomea
265, I-34136 Trieste, Italy<EMAIL_ADDRESS>ZHAONAN DONG 1)
Inria, 2 rue Simone Iff, 75589 Paris, France, 2) CERMICS, Ecole des Ponts,
77455 Marne-la-Vallée, France<EMAIL_ADDRESS>Emmanuil H. Georgoulis
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 (), 2) Department of
Mathematics, School of Applied Mathematical and Physical Sciences, National
Technical University of Athens, Zografou 15780, Greece & 3) IACM-FORTH,
Greece<EMAIL_ADDRESS>
###### 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.

###### keywords:

discontinuous Galerkin; _a posteriori_ error bound; polygonal/polyhedral
meshes; polytopic elements; irregular hanging nodes.

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 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 fo 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 functions, 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 employing 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.8 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 applicable 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 interior 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 robustness and efficiency of the of the derived _a posteriori_ 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\geq 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_{\text{\rm D}}$ and
$\Gamma_{\text{\rm N}}$ with $|\Gamma_{\text{\rm D}}|\neq 0$. For technical
reasons, when $d=3$ and $|\Gamma_{\text{\rm N}}|\neq 0$, the interface between
$\Gamma_{\text{\rm D}}$ and $\Gamma_{\text{\rm 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

(1) $\displaystyle-\nabla\cdot(a\nabla u)=$ $\displaystyle\ f$
$\displaystyle\text{in }\Omega,$ $\displaystyle\quad u=$ $\displaystyle\
g_{\text{\rm D}}$ $\displaystyle\text{on }\Gamma_{\text{\rm D}},$
$\displaystyle\quad a\nabla u\cdot\mbox{\boldmath$\rm{n}$}=$ $\displaystyle\
g_{\text{\rm N}}$ $\displaystyle\text{on }\Gamma_{\text{\rm N}},$

with known $f\in L_{2}(\Omega)$, $g_{\text{\rm D}}\in
H^{1/2}(\Gamma_{\text{\rm D}})$ and $g_{\text{\rm N}}\in
L_{2}(\Gamma_{\text{\rm N}})$ and symmetric diffusion tensor
$a\in[L_{\infty}(\Omega)]^{d\times d}$ such that

(2)
$\alpha^{*}|\xi|^{2}\geq\xi^{\top}a(x)\xi\geq\alpha_{*}|\xi|^{2}>0\quad\forall\xi\in\mathbb{R}^{d},\quad\text{a.e.}\quad
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^{1}_{\text{\rm D}}:=\\{v\in H^{1}(\Omega):v=0\text{ on
}\Gamma_{\text{\rm D}}\\}$, the weak formulation of (1) is: find $u\in
H^{1}(\Omega)$, $u=g_{\text{\rm D}}$ on $\Gamma_{\text{\rm D}}$ such that

(3) $\int_{\Omega}a{\nabla u}\cdot{\nabla
v}\,\mathrm{d}\mathbf{x}=\int_{\Omega}fv\,\mathrm{d}\mathbf{x}+\int_{\Gamma_{\text{\rm
N}}}g_{\text{\rm N}}v\,\mathrm{d}s,$

for all $v\in H^{1}_{\text{\rm D}}(\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}=\bar{\Omega}$. Given
$h_{K}:=\operatorname{diam}(K)$, the diameter of $K\in\mathcal{T}$, we define
the mesh-function
$\mbox{\boldmath$\rm{h}$}:\cup_{K\in\mathcal{T}}K\to\mathbb{R}_{+}$ by
$\mbox{\boldmath$\rm{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{\rm int}}:=\Gamma\backslash\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$.

Figure 1: Element $K\in\mathcal{T}$ with its face-wise neighbours; hanging
nodes are highlighted by a bullet.

$K$

The _finite element space_ $S_{\mathcal{T}}$ with respect to ${\mathcal{T}}$
is defined by

$S_{\mathcal{T}}\equiv S^{p}_{\mathcal{T}}:=\\{u\in
L_{2}(\Omega):u|_{K}\in\mathcal{P}_{p}(K),K\in\mathcal{T}\\},$

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{\rm int}}$.
For $v$ and $\rm{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}})$,
$\\{\\!\\!\\{\mbox{\boldmath$\rm{q}$}\\}\\!\\!\\}|_{F}:=\frac{1}{2}(\mbox{\boldmath$\rm{q}$}|_{F\cap
K_{i}}+\mbox{\boldmath$\rm{q}$}|_{F\cap K_{j}})$, respectively, and the _jump_
across $F$ by $[\\![v]\\!]:=v|_{F\cap K_{i}}-v|_{F\cap K_{j}}$,
$[\\![\mbox{\boldmath$\rm{q}$}]\\!]:=\mbox{\boldmath$\rm{q}$}|_{F\cap
K_{i}}-\mbox{\boldmath$\rm{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_{\text{\rm D}}$, with $F\subset\partial K_{i}$,
$K_{i}\in\mathcal{T}$, we set $\\{\\!\\!\\{v\\}\\!\\!\\}:=v_{i},$
$\\{\\!\\!\\{\mbox{\boldmath$\rm{q}$}\\}\\!\\!\\}:=\mbox{\boldmath$\rm{q}$}_{i},$
$[\\![v]\\!]:=v_{i},$ and
$[\\![\mbox{\boldmath$\rm{q}$}]\\!]:=\mbox{\boldmath$\rm{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_{{\rm T}}v$ the face-wise
tangential gradient operator acting on the traces of $v$ on $\Gamma$, noting
that $\nabla_{{\rm T}}v$ is double-valued on $\Gamma_{\text{\rm 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_{\text{\rm 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_{\text{\rm D}}$ and
$\Gamma_{\text{\rm 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  (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}\geq\tau^{-1}h_{K}$, centred at some point
$\mathbf{x}^{0}\equiv\mathbf{x}^{0}_{K}\in K$; see Figure 2 (right panel) for
an illustration when $d=2$.

  * (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}\geq\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.1.

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.12 below.

Assumption 2.3 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|<<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.3(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.9 below.

To the best of our knowledge, Assumption 2.3 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.

---

$K$${\bf x}_{1}^{0}$$K_{1}$${\bf x}_{2}^{0}$$K_{2}$$K$${\bf x}^{0}$

Figure 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
${\bf x}_{1}^{0}$ and ${\bf 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 ${\bf x}^{0}$.

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

######  Assumption 2.2 (local quasi-uniformity).

Each mesh ${\mathcal{T}}$ is locally quasiuniform, i.e., there exists a
constant $\rho>0$ such that $\rho^{-1}\leq h_{K_{i}}/h_{K_{j}}\leq\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.3.

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

######  Proof 2.4.

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^{\prime}\in{\omega}_{K}$, then
$K^{\prime}\subset B_{h_{K}+h_{K^{\prime}}}(\mathbf{x}^{0}_{K})\subseteq
B_{(1+\rho)h_{K}}(\mathbf{x}^{0}_{K})$ thanks to Assumption 2.2. On the other
hand, Assumption 2.3 implies that $K^{\prime}$ contains the ball
$B_{r_{K^{\prime}}}(\mathbf{x}^{0}_{K^{\prime}})$ with
$r_{K^{\prime}}\geq\tau^{-1}h_{K^{\prime}}$. Letting $C_{2}:=\pi$ and
$C_{3}:=4\pi/3$, we thus have $|K^{\prime}|\geq C_{d}r_{K^{\prime}}^{d}\geq
C_{d}\tau^{-d}h_{K^{\prime}}^{d}\geq C_{d}\tau^{-d}\rho^{-d}h_{K}^{d}$.
Therefore, $n(C_{d}\tau^{-d}\rho^{-d}h_{K}^{d})\leq
C_{d}(1+\rho)^{d}h_{K}^{d}$, or $n\leq\tau^{d}\rho^{d}(1+\rho)^{d}$, thereby
showing that the number of neighbours of $K$ is uniformly bounded as required.

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.3. Then, for all $v\in
H^{1}(\Omega)$, we have the trace estimate

(4) $\|v\|_{\partial K}^{2}\leq C_{{\rm
tr}}\left(\frac{\zeta}{h_{K}}\|v\|_{K}^{2}+\frac{h_{K}}{\zeta}\|\nabla
v\|_{K}^{2}\right),$

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

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

(5) $\|{v}\|_{\partial K}^{2}\leq\frac{C_{{\rm inv}}}{h_{K}}\|{v}\|_{K}^{2},$

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

######  Proof 2.6.

These estimates are special cases of the corresponding ones presented in [10].
The trace estimate (4) follows from [10, Lemma 4.7] in conjunction with
Assumption 2.3, while (5) follows from [10, Lemma 4.4] along with Assumption
2.3.

######  Lemma 2.7.

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

(6) $\|{v-\Pi_{0}v}\|_{{K}}\leq C_{PF}h_{K}\|{\nabla v}\|_{{K}}$

and

(7) $\|{v-\Pi_{0}v}\|_{{\partial K}}\leq\tilde{C}_{PF}\sqrt{h_{K}}\|{\nabla
v}\|_{{K}},$

with $\Pi_{0}:L_{2}(\Omega)\to S_{\mathcal{T}}^{0}$, the orthogonal
$L_{2}$-projection onto $S_{\mathcal{T}}^{0}$, the space of element-wise
constants; here ${C}_{PF},\tilde{C}_{PF}>0$ depend on $d$ and on $\tau$ only.

######  Proof 2.8.

The key technical difficulty is to show that ${C}_{PF}$ and $\tilde{C}_{PF}$
are independent of the shape of $K$. Under Assumption 2.3, 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 (6). Then,
(7) follows using the trace inequality (4).

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

######  Definition 2.9 (Auxiliary mesh).

Given the sequence of meshes ${\mathcal{T}}$, we name _auxiliary mesh_ a
corresponding sequence of conforming simplicial meshes
$\widehat{{\mathcal{T}}}$ satisfying: for each $\widehat{{\mathcal{T}}}$ and
$T\in\widehat{{\mathcal{T}}}$

  * (a)

(shape-regularity) the radius $r_{T}$ of the largest circle inscribed in $T$
is such that $r_{T}\geq\widehat{\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 $\widehat{\rho}^{-1}\leq
h_{T}/h_{K}\leq\widehat{\rho}$,

with $\widehat{\tau},\widehat{\rho}>1$ constants independent of the
discretization parameters.

An immediate consequence of the above definition is that, if
$\widehat{{\mathcal{T}}}$ is an auxiliary mesh sequence, then the number of
intersections of each $T\in\widehat{{\mathcal{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.3.

We note that the evaluation of the error estimator presented below does _not_
require 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.3 and 2.2.
Rather, the issue is to show that auxiliary meshes tightly close to the
polytopic meshes can be constructed in principle. 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 corresponding
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.3(b) for boundary
faces. Then, an auxiliary mesh can be simply constructed by joining $F$ to
$\mathbf{x}^{0}_{K}$, 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.10 (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
$\mathcal{C}$ such that:

  * (1)

the vertices of $S$ are on the boundary of $\mathcal{C}$;

  * (2)

if a vertex $P$ of the CDT is in the interior of $\mathcal{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
$\mathcal{C}$ of Definition 2.10 does not _enclose_ any other point in
$\mathcal{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}^{0}_{K}\\}_{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.11.

We expect the CDTs associated to ${\mathcal{T}}$ to always satisfy the
auxiliary mesh Definition 2.9 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.9 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].

---
Figure 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.9 are marked with a grey bullet.

######  Remark 2.12.

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.13 (Quasi-interpolant).

Let $\widehat{{\mathcal{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)\rightarrow H^{1}(\Omega)\cap
S^{1}_{\widehat{\mathcal{T}}}$, such that

(8) $\|{v-I_{h}v}\|_{{T}}+{h_{T}}\|{\nabla(v-I_{h}v)}\|_{{T}}\leq C_{\rm
sz}{h_{T}}\|{\nabla v}\|_{{{{\omega}}_{T}}},$

where ${\omega}_{T}$ denotes the patch of elements in
$\widehat{{\mathcal{T}}}$ with non-empty intersection with $\overline{T}$. The
constant $C_{\rm sz}>0$ depends only on the shape-regularity constant
$\widehat{\tau}$ of the auxiliary mesh $\widehat{{\mathcal{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\widehat{{\mathcal{T}}}$.

Moreover, we have

(9) $\displaystyle\|{v-I_{h}v}\|_{{\partial K}}^{2}\leq C_{\rm
in}{h_{K}}\|{\nabla v}\|_{{\widehat{\omega}_{K}}}^{2},$

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

######  Proof 2.14.

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

$\displaystyle\|{v-I_{h}v}\|_{{\partial K}}^{2}$ $\displaystyle\leq$
$\displaystyle C_{{\rm
tr}}\left(h_{K}^{-1}\|{v-I_{h}v}\|_{{K}}^{2}+{h_{K}}\|{\nabla(v-I_{h}v)}\|_{{K}}^{2}\right)$
$\displaystyle\leq$ $\displaystyle C_{{\rm
tr}}\sum_{T\in\widehat{{\mathcal{T}}}:{T}\cap{K}\neq\emptyset}\left(h_{K}^{-1}\|{v-I_{h}v}\|_{{T}}^{2}+{h_{K}}\|{\nabla(v-I_{h}v)}\|_{{T}}^{2}\right),$

and (9) follows from (8), depending on the shape-regularity of ${\mathcal{T}}$
through $C_{tr}$ and on the shape-regularity of $\widehat{{\mathcal{T}}}$
through (8).

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

(10) $\Gamma_{T}:=\Gamma\cap T=\Gamma_{T}^{{\rm c}}\cup\Gamma_{T}^{{\rm o}},$

with $\Gamma_{T}^{{\rm c}}$ the set of interfaces
$\gamma_{T}\in\Gamma_{T}\cap\Omega$ such that $\gamma_{T}\subset\partial K$
with $K$ an element of ${\mathcal{T}}$ whose centre of star-shapedness is a
vertex of $T$, and $\Gamma_{T}^{{\rm o}}:=\Gamma_{T}\setminus\Gamma_{T}^{{\rm
c}}$. Note that the number of interfaces in $\Gamma_{T}$ is bounded since the
number of intersections of $T$ with the elements $K\in\mathcal{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 (10) reflects increasing levels of difficulty, with
$\Gamma_{T}^{{\rm 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}^{{\rm 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.

######  Remark 2.15.

In the case of sub-mesh auxiliary meshes, we always have $\Gamma_{T}^{{\rm
o}}=\emptyset$. Instead, for constrained Delauney auxiliary meshes,
$\Gamma_{T}^{{\rm o}}\neq\emptyset$ in general; we refer to Figure 4 for some
examples.

$T$$T$$T$$\gamma_{T}$

Figure 4: Constrained Delauney auxiliary mesh examples. Different
configurations for the set $\Gamma_{T}$ of (10) for the auxiliary element $T$
marked in grey. Left: $\Gamma_{T}=\Gamma_{T}^{{\rm c}}$. Centre:
$\Gamma_{T}=\Gamma_{T}^{{\rm o}}$. Right: $\Gamma_{T}^{{\rm o}}$ is made of a
single interface $\gamma_{T}$ while the rest of $\Gamma_{T}$ belongs to
$\Gamma_{T}^{{\rm c}}$. The simplices required in the proof of Lemma 2.16 to
estimate the contributions from each edge in $\gamma_{T}$ partially overlap
(not shown).

######  Lemma 2.16.

Let Assumptions 2.3 and 2.2 hold and let $\widehat{{\mathcal{T}}}$ be an
auxiliary mesh, related to ${\mathcal{T}}$. Given $q\in\mathbb{N}$, for all
$T\in\widehat{{\mathcal{T}}}$ and $v\in\mathcal{P}_{q}(T)$, we have

(11) $\|{v}\|_{{\Gamma_{T}}}^{2}\leq\frac{\widehat{C}_{{\rm
inv}}}{h_{T}}\|{v}\|_{T}^{2}.$

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

######  Proof 2.17.

We consider $\Gamma_{T}^{{\rm c}}$ and $\Gamma_{T}^{{\rm o}}$ separately,
starting with $\Gamma_{T}^{{\rm c}}$. Exploiting the star-shapedness property
of $\Gamma_{T}^{{\rm c}}$ with respect to the vertices of $T$, which is
inherited from Assumption 2.3, the inverse inequality

$\|{v}\|_{{\Gamma_{T}^{{\rm
c}}}}^{2}\leq\frac{\tau{{(q+1)(q+d)}}}{h_{T}}\|{v}\|_{T}^{2},$

follows in the same way as (5).

Considering now the set $\Gamma_{T}^{{\rm o}}$, we observe that each cut
interface belonging to $\Gamma_{T}^{{\rm o}}$ may be partitioned into a set of
$(d-1)$-dimentional 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}^{{\rm 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 $\widehat{F}$ be one such ($d-1$)-dimentional simplex within
$\Gamma_{T}^{{\rm 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)\geq r$. Otherwise, $S$ must be contained in the region
$\\{{\bf x}\in{\mathbb{R}}^{d}:d({\bf 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_{\widehat{F}}$ by joining
$\widehat{F}$ with a vertex $V$ of $T$ such that $d(V,\widehat{F})\geq r$. We
thus have, cf. (5), the inverse estimate:

(12)
$\|{v}\|_{{\widehat{F}}}^{2}\leq\widehat{\tau}\frac{(q+1)(q+d)}{h_{T}}\|{v}\|_{T_{\widehat{F}}}^{2}\qquad\forall
v\in\mathcal{P}_{q}(T).$

Then, summing up over all $\widehat{F}\in\Gamma_{T}^{{\rm o}}$ we conclude

(13) $\|{v}\|_{{\Gamma_{T}^{{\rm
o}}}}^{2}\leq\frac{\widehat{\tau}{{(q+1)(q+d)}}}{h_{T}}\sum_{\widehat{F}\in\Gamma_{T}^{{\rm
o}}}\|{v}\|_{T_{\widehat{F}}}^{2}\leq n_{{\rm
o}}\frac{\widehat{\tau}{{(q+1)(q+d)}}}{h_{T}}\|{v}\|_{T}^{2},$

with $n_{{\rm o}}$ the number of overlaps of the simplices $T_{\widehat{F}}$,
$\widehat{F}\in\Gamma_{T}^{{\rm o}}$. In particular, if $\Gamma_{T}^{{\rm o}}$
is only made of boundary interfaces, then $n_{{\rm o}}\leq 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}^{{\rm c}}$ and $\Gamma_{T}^{{\rm o}}$.

######  Remark 2.18.

The constant $\widehat{C}_{{\rm inv}}$ appearing in (11) accounts for the
complexity of the mesh in terms of topology and shape, quantified by the
number $n_{{\rm 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}^{{\rm 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}^{{\rm o}}$ is either
empty or it contains a single boundary edge, i.e., no overlaps are required,
resulting in the ‘ideal’ constant $\widehat{C}_{{\rm inv}}=\tau{{(q+1)(q+d)}}$
for each $T\in\widehat{{\mathcal{T}}}$.

###  2.4 Discontinuous Galerkin method

Let $\mathcal{V}:=S_{\mathcal{T}}+H^{1}(\Omega)$. The symmetric interior
penalty discontinuous Galerkin method reads: find $u_{h}\in S_{\mathcal{T}}$
such that

(14) $B(u_{h},v_{h})=\ell(v_{h})\qquad\text{for all }\quad v_{h}\in
S_{\mathcal{T}},$

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

(15) $\displaystyle B(w,v):=$
$\displaystyle\int_{\Omega}a\nabla_{h}w\cdot\nabla_{h}v\,\mathrm{d}\mathbf{x}+\int_{\Gamma\backslash\Gamma_{\text{\rm
N}}}\sigma[\\![w]\\!][\\![v]\\!]\,\mathrm{d}{s},$
$\displaystyle-\int_{\Gamma\backslash\Gamma_{\text{\rm
N}}}(\\{\\!\\!\\{a(\mbox{\boldmath$\rm{\Pi}$}\nabla
w)\cdot\mbox{\boldmath$\rm{n}$}\\}\\!\\!\\}[\\![v]\\!]+\\{\\!\\!\\{a(\mbox{\boldmath$\rm{\Pi}$}\nabla
v)\cdot\mbox{\boldmath$\rm{n}$}\\}\\!\\!\\}[\\![w]\\!])\,\mathrm{d}{s},$

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

$\displaystyle\ell(v):=\int_{\Omega}fv\,\mathrm{d}\mathbf{x}-\int_{\Gamma_{\text{\rm
D}}}g_{\rm D}\big{(}(a(\mbox{\boldmath$\rm{\Pi}$}\nabla
v))\cdot\mbox{\boldmath$\rm{n}$}-\sigma
v\big{)}\,\mathrm{d}s+\int_{\Gamma_{\text{\rm N}}}g_{\text{\rm
N}}v\,\mathrm{d}s,$

with $\mbox{\boldmath$\rm{\Pi}$}:[L_{2}(\Omega)]^{d}\to[S_{\mathcal{T}}]^{d}$
denoting the orthogonal $L_{2}$-projection operator onto the (vectorial)
finite element space, and $\sigma\in
L_{\infty}(\Gamma\backslash\Gamma_{\text{\rm N}})$ being the, so-called,
_discontinuity-penalization function_ given by

(16)
$\sigma({\mathbf{x}}):=\left\\{\begin{array}[]{ll}C_{\sigma}{\displaystyle\max_{K\in\\{K_{i},K_{j}\\}}\Big{\\{}\frac{\bar{a}_{K}C_{{\rm
inv}}}{h_{K}}\Big{\\}}},&~{}~{}\mathbf{x}\in F\in\Gamma_{\text{\rm
int}},~{}F\subseteq\partial K_{i}\cap\partial K_{j},\\\ \\\
C_{\sigma}\displaystyle\frac{\bar{a}_{K}C_{{\rm
inv}}}{h_{K}},&~{}~{}\mathbf{x}\in F\in\Gamma_{\text{\rm
D}},~{}F\subset\partial K,\end{array}\right.$

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_{{\rm inv}}$ for brevity; see
[12, 10] for details. Note that, using (3) and that $[\\![v]\\!]=0$ on
$\Gamma\backslash\Gamma_{\text{\rm N}}$ for all $v\in H^{1}_{\text{\rm D}}$,
we have $B(u,v)=\ell(v)$ for all $v\in H^{1}_{\text{\rm D}}$, with $u\in
H^{1}(\Omega)$ the solution to (3).

######  Remark 2.19.

To avoid further notational overhead, we opted in exposing the main results
for element-wise constant diffusion tensors, i.e., $a\in(S^{{\bf
0}}_{\mathcal{T}})^{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
$|\kern-0.7113pt\|{v}|\kern-0.7113pt\|:=\big{(}\|{\sqrt{a}\nabla_{h}v}\|^{2}+\|{\sqrt{\sigma}[\\![v]\\!]}\|_{{\Gamma\backslash\Gamma_{\text{\rm
N}}}}^{2}\big{)}^{1/2}$, we have the following result.

######  Lemma 2.20.

Under Assumption 2.3, there exists $C_{\sigma}>0$, such that

(17) $B(v,v)\geq C_{\rm
coer}|\kern-0.7113pt\|{v}|\kern-0.7113pt\|^{2}\quad\text{and}\quad B(w,v)\leq
C_{\rm
cont}|\kern-0.7113pt\|{w}|\kern-0.7113pt\|\,|\kern-0.7113pt\|{v}|\kern-0.7113pt\|\quad\text{for
all}\quad v\in\mathcal{V},$

respectively, with $C_{\rm coer},C_{\rm 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.3 and 2.2 and that an auxiliary
mesh according to Definition 2.10 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_{\mathcal{T}}$, defined by

(18) $B(u_{c},v)=B(u_{h},v)\qquad\forall v\in H^{1}_{\text{\rm D}}(\Omega),$

and $u_{c}=g_{\text{\rm D}}$ on $\Gamma_{\text{\rm 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
$\mbox{\boldmath$\rm{w}$}\in(L_{2}(\Omega))^{d}$, there exists $\xi\in
H^{1}_{\text{\rm D}}(\Omega)$ and $\phi\in[H^{1}(\Omega)]^{2d-3}$, $d=2,3$,
such that

(19) $a{\bf w}=a\nabla\xi+\emph{\text{curl}}\,\phi\quad\text{in }\,\Omega,$

and $\phi$ can be chosen so that

(20) $\emph{\text{curl}}\,\phi\cdot\mbox{\boldmath$\rm{n}$}=0\quad\text{on
}\,\Gamma_{\text{\rm N}}.$

Moreover, the following relations hold

(21) $\|{\sqrt{a}{\bf
w}}\|^{2}=\|{\sqrt{a}\nabla\xi}\|^{2}+\|{a^{-1/2}\emph{\text{curl}}\phi}\|^{2},$

and

(22) $\|{\nabla\phi}\|\leq C_{\Omega}\|{\emph{\text{curl}}\phi}\|,$

with a constant $C_{\Omega}>0$ only depending on $\Omega$.

######  Proof 3.3.

The proof of (19) 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
${\text{curl}}\,\phi$, from the symmetry of the diffusion tensor $a$, the
orthogonality (21) follows immediately. Finally, the proof of (22) can be
found in [16]; for $d=2$, we have $C_{\Omega}=1$.

######  Remark 3.4.

The Helmholtz decomposition can be generalised to multiply connected domains
[31]. However, concerning the validity in this setting of the relation (22),
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 (20) imposes a constraint on $\Gamma_{\text{\rm N}}$ for
$\phi\in[H^{1}(\Omega)]^{2d-3}$. Namely,
$\text{curl}\,\phi\cdot\mbox{\boldmath$\rm{n}$}_{F}=0$ on
$F\subset\Gamma_{\text{\rm N}}$, implying that $\phi\in[H^{1}(\Omega)]^{3}$
has constant components on each $(d-1)$-dimensional planar subset of
$\Gamma_{\text{\rm N}}$.

We apply the Helmholtz decomposition with
$\mbox{\boldmath$\rm{w}$}=\nabla_{h}e_{d}$. Hence $\xi\in H^{1}_{\text{\rm
D}}(\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

(23) $\displaystyle\|{\sqrt{a}\nabla_{h}e_{d}}\|^{2}=$
$\displaystyle\int_{\Omega}a\nabla_{h}e_{d}\cdot\nabla\xi\,\mathrm{d}\mathbf{x}+\int_{\Omega}\nabla_{h}e_{d}\cdot\text{curl}\,\phi\,\mathrm{d}\mathbf{x}.$

Since $u_{c}\in H^{1}(\Omega)$ with $u_{c}|_{\Gamma_{\text{\rm
D}}}=g_{\text{\rm D}}$ and $\xi\in H^{1}_{\text{\rm D}}(\Omega)$, (18) implies

$\int_{\Omega}a\nabla_{h}e_{d}\cdot\nabla\xi\,\mathrm{d}\mathbf{x}=-\int_{\Gamma_{\text{\rm
int}}}\\{\\!\\!\\{a(\mbox{\boldmath$\rm{\Pi}$}\nabla\xi)\cdot\mbox{\boldmath$\rm{n}$}\\}\\!\\!\\}[\\![u_{h}]\\!]\,\mathrm{d}{s}-\int_{\Gamma_{\text{\rm
D}}}(a(\mbox{\boldmath$\rm{\Pi}$}\nabla\xi)\cdot\mbox{\boldmath$\rm{n}$})(u_{h}-g_{\text{\rm
D}})\,\mathrm{d}{s}.$

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

$\displaystyle\int_{\Omega}a\nabla_{h}e_{d}\cdot\nabla\xi\,\mathrm{d}\mathbf{x}$
$\displaystyle\leq\|{\sigma^{-1/2}\\{\\!\\!\\{a(\mbox{\boldmath$\rm{\Pi}$}\nabla\xi)\\}\\!\\!\\}}\|_{{\Gamma\backslash\Gamma_{\text{\rm
N}}}}\Big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\Big{)}^{1/2}$ (24)
$\displaystyle\leq\Big{(}{C_{\sigma}}^{-1}\sum_{K\in\mathcal{T}}\|{{\mbox{\boldmath$\rm{\Pi}$}(\sqrt{a}\nabla\xi)}}\|_{{K}}^{2}\Big{)}^{1/2}\Big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\Big{)}^{1/2}$
$\displaystyle\leq(C_{\sigma})^{-1/2}\|{{\sqrt{a}\nabla\xi}}\|\Big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\Big{)}^{1/2}$
$\displaystyle\leq\|{{\sqrt{a}\nabla_{h}e_{d}}}\|\Big{(}\|{\sqrt{\sigma
C_{\sigma}^{-1}}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma C_{\sigma}^{-1}}(u_{h}-g_{\text{\rm
D}})}\|_{{\Gamma_{\text{\rm D}}}}^{2}\Big{)}^{1/2}.$

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

(25)
$\displaystyle\int_{\Omega}\nabla_{h}e_{d}\cdot\text{curl}\,\phi\,\mathrm{d}\mathbf{x}=\int_{\Omega}\nabla_{h}e_{d}\cdot\text{curl}(\phi-
I_{h}\phi)\,\mathrm{d}\mathbf{x}+\int_{\Omega}\nabla_{h}e_{d}\cdot\text{curl}\,I_{h}\phi\,\mathrm{d}\mathbf{x},$

with $I_{h}$ the (component-wise if $d=3$) quasi-interpolation operator of
Lemma 2.13.

Starting with the first term, observing that $\text{curl}\nabla
u_{c}=\mbox{\boldmath$\rm{0}$}$ and using the fact that
$\phi\in[H^{1}(\Omega)]^{2d-3}$ satisfying (20), implying that $\phi$ is a
constant function on each planar section of $\Gamma_{\text{\rm N}}$, and
choosing $\phi=I_{h}\phi$ on $\Gamma_{\text{\rm N}}$. Then we have

$\int_{\Omega}\nabla u_{c}\cdot\text{curl}(\phi-
I_{h}\phi)\,\mathrm{d}\mathbf{x}=\int_{\Gamma_{\text{\rm D}}}(\phi-
I_{h}\phi)\nabla_{{\rm T}}g_{\text{\rm D}}\,\mathrm{d}{s}.$

Applying integration by parts, observing that $\text{curl}\nabla
e_{d}|_{K}=\mbox{\boldmath$\rm{0}$}$, and using $(\phi-I_{h}\phi)$ is single
valued on each face, and $(\phi-I_{h}\phi)_{\Gamma_{\text{\rm N}}}=0$, yields

(26) $\displaystyle\sum_{K\in{\mathcal{T}}}\int_{K}\nabla
e_{d}\cdot\text{curl}(\phi-I_{h}\phi)\,\mathrm{d}\mathbf{x}$ $\displaystyle=$
$\displaystyle\ -\int_{\Gamma}(\phi-I_{h}\phi)[\\![\nabla_{{\rm
T}}u_{h}]\\!]\,\mathrm{d}{s}+\int_{\Gamma_{\text{\rm D}}}(\phi-
I_{h}\phi)\nabla_{{\rm T}}g_{\text{\rm D}}\,\mathrm{d}{s}$ $\displaystyle=$
$\displaystyle\ -\int_{\Gamma_{\text{\rm int}}}(\phi-
I_{h}\phi)[\\![\nabla_{{\rm
T}}u_{h}]\\!]\,\mathrm{d}{s}-\int_{\Gamma_{\text{\rm D}}}(\phi-
I_{h}\phi)\nabla_{{\rm T}}(u_{h}-g_{\text{\rm D}})\,\mathrm{d}{s}.$

Further, using (9) and, finally, (22) and (21), the right-hand side of (26)
can be further estimated from above by

(27) $\displaystyle\sum_{K\in\mathcal{T}}\|{h_{K}^{-1/2}(\phi-
I_{h}\phi)}\|_{{\partial K\backslash\Gamma_{\text{\rm
N}}}}\left(\|{\sqrt{h}_{K}[\\![\nabla_{{\rm T}}u_{h}]\\!]}\|_{{\partial
K\cap\Gamma_{\text{\rm int}}}}+\|{\sqrt{h}_{K}\nabla_{{\rm
T}}(u_{h}-g_{\text{\rm D}})}\|_{{\partial K\cap\Gamma_{\text{\rm D}}}}\right)$
$\displaystyle\leq\big{(}\sum_{K\in\mathcal{T}}C_{{\rm
tr}}C_{I}\|{\nabla\phi}\|_{{\widehat{\omega}_{K}}}^{2}\big{)}^{1/2}\big{(}\|{\sqrt{h}[\\![\nabla_{{\rm
T}}u_{h}]\\!]}\|_{{\Gamma_{\text{\rm int}}}}^{2}+\|{\sqrt{h}\nabla_{{\rm
T}}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm D}}}}^{2}\big{)}^{1/2}$
$\displaystyle\leq
C_{1}\|{\nabla\phi}\|_{{\Omega}}\big{(}\|{\sqrt{h}[\\![\nabla_{{\rm
T}}u_{h}]\\!]}\|_{{\Gamma_{\text{\rm int}}}}^{2}+\|{\sqrt{h}\nabla_{{\rm
T}}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm D}}}}^{2}\big{)}^{1/2}$
$\displaystyle\leq
C_{1}C_{\Omega}\sqrt{\alpha^{*}}\|{\sqrt{a}\nabla_{h}e_{d}}\|\big{(}\|{\sqrt{h}[\\![\nabla_{{\rm
T}}u_{h}]\\!]}\|_{{\Gamma_{\text{\rm int}}}}^{2}+\|{\sqrt{h}\nabla_{{\rm
T}}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm D}}}}^{2}\big{)}^{1/2},$

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

We now consider the second term in (25). Since
$I_{h}\phi\in[H^{1}(\Omega)]^{2d-3}$, we have
$\nabla\cdot\text{curl}\,I_{h}\phi=0$ on $\Omega$ and, hence,
$[\\![\text{curl}\,I_{h}\phi\cdot\mbox{\boldmath$\rm{n}$}]\\!]=0$ on
$\Gamma_{\text{\rm int}}$. Moreover, given that $I_{h}\phi$ is constant on
each component of $\Gamma_{\text{\rm N}}$, we also have
${\text{curl}\,I_{h}\phi\cdot\mbox{\boldmath$\rm{n}$}}=0$ on
$\Gamma_{\text{\rm N}}$. Then, integration by parts and working as above gives

$\displaystyle\int_{\Omega}\nabla_{h}e_{d}\cdot\text{curl}\,I_{h}\phi\,\mathrm{d}\mathbf{x}=$
$\displaystyle\sum_{K\in{\mathcal{T}}}\int_{\partial
K}e_{d}\mbox{\boldmath$\rm{n}$}_{K}\cdot\text{curl}\,I_{h}\phi\,\mathrm{d}{s}$
$\displaystyle=$ $\displaystyle-\int_{\Gamma_{\text{\rm
int}}}[\\![u_{h}]\\!](\mbox{\boldmath$\rm{n}$}\cdot\text{curl}\,I_{h}\phi)\,\mathrm{d}{s}-\int_{\Gamma_{\text{\rm
D}}}(u_{h}-g_{\text{\rm
D}})(\mbox{\boldmath$\rm{n}$}\cdot\text{curl}\,I_{h}\phi)\,\mathrm{d}{s}.$

Next, applying to $I_{h}\phi$ the trace inverse inequality with respect to the
auxiliary mesh $\widehat{{\mathcal{T}}}$ given in (11), we obtain

$\displaystyle\int_{\Omega}\nabla_{h}e_{d}\cdot\text{curl}\,I_{h}\phi\,\mathrm{d}\mathbf{x}$
$\displaystyle\leq\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}\Big{(}\sum_{T\in\widehat{{\mathcal{T}}}}\|{{\sigma^{-1/2}}\text{curl}\,I_{h}\phi}\|_{{T\cap(\Gamma\backslash\Gamma_{\text{\rm
N}})}}^{2}\Big{)}^{1/2}$
$\displaystyle\leq\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}\Big{(}\sum_{T\in\widehat{{\mathcal{T}}}}{\widehat{C}_{{\rm
inv}}h_{T}^{-1}}\|{{\sigma^{-1/2}}\text{curl}\,I_{h}\phi}\|_{{T}}^{2}\Big{)}^{1/2}$
$\displaystyle\leq C_{2}\big{(}\|{\sqrt{\sigma
C_{\sigma}^{-1}}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma C_{\sigma}^{-1}}(u_{h}-g_{\text{\rm
D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}\|{\text{curl}\,I_{h}\phi}\|,$

for $C_{2}>0$ constant depending on $C_{{\rm inv}}$, $\widehat{C}_{{\rm
inv}}$, $\widehat{\rho}$, and on $\alpha^{*}$. Next, we use the stability of
$I_{h}$ from Lemma 2.13, together with (22) to deduce

(28) $\|{\text{curl}\,I_{h}\phi}\|_{{\Omega}}\leq
C_{c}\|{\nabla\phi}\|_{{\Omega}}\leq
C_{c}C_{\Omega}\sqrt{\alpha_{*}}\|{a^{-1/2}\nabla\phi}\|\leq
C_{c}C_{\Omega}\sqrt{\alpha_{*}}\|{\sqrt{a}\nabla_{h}e_{d}}\|.$

Hence, combing (3.1), (26), (LABEL:eq:nonc_second_term_II), and (28), we
arrive at the bound

$\displaystyle\|{\sqrt{a}\nabla_{h}e_{d}}\|^{2}$ $\displaystyle\leq C_{\rm
nc}\Big{(}\|{\sqrt{\sigma
C_{\sigma}^{-1}}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma C_{\sigma}^{-1}}(u_{h}-g_{\text{\rm
D}})}\|_{{\Gamma_{\text{\rm D}}}}^{2}$ (29) $\displaystyle\hskip
28.45274pt+\|{\sqrt{h}[\\![\nabla_{{\rm T}}u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{h}\nabla_{{\rm T}}(u_{h}-g_{\text{\rm
D}})}\|_{{\Gamma_{\text{\rm D}}}}^{2}\Big{)};$

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

######  Remark 3.5.

We stress that the mesh-size $h$ in (3.1) 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
$\operatorname{diam}(K)/\operatorname{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^{1}_{\text{\rm D}}$, we have

(30) $\displaystyle B(e,v)$
$\displaystyle=\ell(v)-B(u_{h},v)=\ell(\eta)-B(u_{h},\eta),$

with $\eta:=v-v_{h}$, for any $v_{h}\in S_{\mathcal{T}}$. Recalling that
$e=e_{c}+e_{d}$, since $e_{c}\in H^{1}_{\text{\rm D}}$ we can fix $v=e_{c}$ in
(30) to further deduce

(31) $\|{\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 (18). The right-hand side of (31) can now be bounded via standard
arguments [35]: integration by parts, application of [4, Eq. (3.3)], the
observation that $\mbox{\boldmath$\rm{\Pi}$}\nabla u_{h}=\nabla u_{h}$, and
elementary manipulations yield

$\displaystyle\ell(\eta)-B(u_{h},\eta)$
$\displaystyle=\int_{\Omega}(f+\nabla_{h}\cdot(a\nabla_{h}u_{h}))\eta\,\mathrm{d}\mathbf{x}-\int_{\Gamma_{\text{\rm
int}}}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]\\{\\!\\!\\{\eta\\}\\!\\!\\}\,\mathrm{d}{s}$
$\displaystyle-\int_{\Gamma_{\text{\rm N}}}(a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}-g_{\text{\rm
N}})\eta\,\mathrm{d}s-\int_{\Gamma_{\text{\rm
int}}}\sigma[\\![u_{h}]\\!][\\![\eta]\\!]\,\mathrm{d}{s}-\int_{\Gamma_{\text{\rm
D}}}\sigma(u_{h}-g_{\rm D})\eta\,\mathrm{d}s$ (32)
$\displaystyle+\int_{\Gamma_{\text{\rm
int}}}\\{\\!\\!\\{a(\mbox{\boldmath$\rm{\Pi}$}\nabla\eta)\cdot\mbox{\boldmath$\rm{n}$}\\}\\!\\!\\}[\\![u_{h}]\\!]\,\mathrm{d}{s}+\int_{\Gamma_{\text{\rm
D}}}\big{(}a(\mbox{\boldmath$\rm{\Pi}$}\nabla\eta)\cdot\mbox{\boldmath$\rm{n}$}\big{)}(u_{h}-g_{\rm
D})\,\mathrm{d}s=:\sum_{i=1}^{7}T_{i}.$

Setting $\eta=e_{c}-\Pi_{0}e_{c}$ and using (6), we have

(33) $T_{1}\leq
C_{PF}\alpha_{*}^{-1/2}\Big{(}\sum_{K\in{\mathcal{T}}}\|{h_{K}(f+\nabla\cdot(a\nabla
u_{h})}\|_{{K}}^{2}\Big{)}^{1/2}\|{\sqrt{a}\nabla e_{c}}\|_{{.}}$

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

$\displaystyle T_{2}+T_{3}$
$\displaystyle\leq\sum_{K\in\mathcal{T}}\sum_{F\subset\partial
K\cap\Gamma_{\text{\rm
int}}}\|{h_{K}^{-1/2}\eta}\|_{{F}}\|{\sqrt{h}_{K}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{F}}$ (34)
$\displaystyle+\sum_{K\in\mathcal{T}}\sum_{F\subset\partial
K\cap\Gamma_{\text{\rm N}}}\|{h_{K}^{-1/2}\eta}\|_{{F}}\|{\sqrt{h}_{K}(a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}-g_{\text{\rm N}})}\|_{{F}}$
$\displaystyle\leq\tilde{C}_{PF}\alpha_{*}^{-1/2}\|{\sqrt{a}\nabla
e_{c}}\|\Big{(}\|{\sqrt{h}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{h}(a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}-g_{\text{\rm N}})}\|_{{\Gamma_{\text{\rm
N}}}}^{2}\Big{)}^{1/2}.$

Similarly, using the definition of $\sigma$ from (16), we have

$\displaystyle T_{4}+T_{5}$
$\displaystyle\leq\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}\|{{\sigma^{1/2}}[\\![\eta]\\!]}\|_{{\Gamma\backslash\Gamma_{\text{\rm
N}}}}$ (35) $\displaystyle\leq\Big{(}\sum_{K\in\mathcal{T}}\max_{F\in\partial
K\backslash\Gamma_{\text{\rm N}}}\sigma\|{\eta}\|_{{\partial
K\backslash\Gamma_{\text{\rm
N}}}}^{2}\Big{)}^{1/2}\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{{\sqrt{\sigma}}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}$
$\displaystyle\leq\sqrt{C_{\sigma}\alpha^{*}/\alpha_{*}}C_{{\rm
inv}}\tilde{C}_{PF}\rho\|{\sqrt{a}\nabla
e_{c}}\|\big{(}\|{{\sigma}^{1/2}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}.$

Now, using the trace inverse estimate (5), the stability of the
$L_{2}$-projection operator and that $\nabla\eta|_{K}=\nabla e_{c}|_{K}$, we
deduce

$\displaystyle T_{6}+T_{7}$
$\displaystyle\leq\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}\|{{\sigma^{-1/2}}\\{\\!\\!\\{a(\mbox{\boldmath$\rm{\Pi}$}\nabla
e_{c})\cdot\mbox{\boldmath$\rm{n}$}\\}\\!\\!\\}}\|_{{\Gamma\backslash\Gamma_{\text{\rm
N}}}}$ $\displaystyle\leq\Big{(}\sum_{K\in\mathcal{T}}\max_{F\in\partial
K\backslash\Gamma_{\text{\rm
N}}}\sigma^{-1}\|{\sqrt{a}(\mbox{\boldmath$\rm{\Pi}$}(\sqrt{a}\nabla
e_{c}))\cdot\mbox{\boldmath$\rm{n}$}}\|_{{F}}^{2}\Big{)}^{1/2}$ (36)
$\displaystyle\hskip
28.45274pt\times\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}$ $\displaystyle\leq
C_{\sigma}^{-1/2}\rho\|{\sqrt{a}\nabla
e_{c}}\|\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}\big{)}^{1/2}.$

Hence, by collecting above bounds (33), (3.2), (3.2), (3.2) and (3.2), we
arrive at the following bound on the conforming error:

$\displaystyle\|{\sqrt{a}\nabla e_{c}}\|$ $\displaystyle\leq C_{\rm
co}\Big{(}\|{h(f+\nabla_{h}\cdot(a\nabla_{h}u_{h})}\|^{2}+\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-g_{\text{\rm D}})}\|_{{\Gamma_{\text{\rm
D}}}}^{2}$ (37) $\displaystyle\hskip 28.45274pt+\|{\sqrt{h}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{h}(a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}-g_{\text{\rm N}})}\|_{{\Gamma_{\text{\rm
N}}}}^{2}\Big{)}^{1/2},$

with $C_{\rm 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.6 (upper bound).

Let $u$ be the solution of (1) and let $u_{h}\in S_{\mathcal{T}}$ be its dG
approximation on a polytopic mesh satisfying Assumptions 2.3 and 2.2. Also let
an auxiliary mesh according to Definition 2.10 is given. Then, we have the
following a posteriori error bound

(38) $\displaystyle|\kern-0.7113pt\|{u-u_{h}}|\kern-0.7113pt\|^{2}$
$\displaystyle\leq C_{\rm 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

$\displaystyle R_{K,E}:=\|{h(\Pi f+\nabla\cdot(a\nabla u_{h}))}\|_{{K}},$
$\displaystyle R_{K,N}:=\big{(}\|{\sqrt{h}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{\partial K\cap\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{h}(a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}-\bar{g}_{\text{\rm N}})}\|_{{\partial
K\cap\Gamma_{\text{\rm N}}}}^{2}\big{)}^{1/2},$ $\displaystyle
R_{K,J}:=\big{(}\|{\sqrt{\sigma}[\\![u_{h}]\\!]}\|_{{\partial
K\cap\Gamma_{\text{\rm int}}}}^{2}+\|{\sqrt{\sigma}(u_{h}-\bar{g}_{\text{\rm
D}})}\|_{{\partial K\cap\Gamma_{\text{\rm D}}}}^{2}\big{)}^{1/2},$
$\displaystyle R_{K,T}:=\big{(}\|{\sqrt{h}[\\![\nabla_{{\rm
T}}u_{h}]\\!]}\|_{{\partial K\cap\Gamma_{\text{\rm
int}}}}^{2}+\|{\sqrt{h}\nabla_{{\rm T}}(u_{h}-\bar{g}_{\text{\rm
D}})}\|_{{\partial K\cap\Gamma_{\text{\rm D}}}}^{2}\big{)}^{1/2},$
$\displaystyle O_{K,E}:=\|{h(f-\Pi f)}\|_{{K}},\quad
O_{K,N}:=\|{\sqrt{h}(g_{\text{\rm N}}-\bar{g}_{\text{\rm N}})}\|_{{\partial
K\cap\Gamma_{\text{\rm N}}}},$ $\displaystyle
O_{K,J}:=\|{\sqrt{\sigma}(g_{\text{\rm D}}-\bar{g}_{\text{\rm
D}})}\|_{{\partial K\cap\Gamma_{\text{\rm D}}}},\quad
O_{K,T}:=\|{\sqrt{h}\nabla_{{\rm T}}(g_{\text{\rm D}}-\bar{g}_{\text{\rm
D}})}\|_{{\Gamma_{\text{\rm D}}}},$

with $C_{\rm up}$ depending on $C_{\rm co}$ and $C_{\rm 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
$\bar{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 3.7.

The proof follows immediately from the bounds (3.1) and (3.2), together with
the triangle inequality
$|\kern-0.7113pt\|{u-u_{h}}|\kern-0.7113pt\|\leq|\kern-0.7113pt\|{e_{c}}|\kern-0.7113pt\|+|\kern-0.7113pt\|{e_{d}}|\kern-0.7113pt\|$.

######  Remark 3.8.

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.9.

Theorem 3.6 has been proven under Assumption 2.3(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{{\mathcal{T}}}$ would correspond ‘closely’ to the primal mesh
${\mathcal{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.6. Of particular interest is the extend 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.3. Then, for each
$K\in\mathcal{T}$, $p\in\mathbb{N}$ and $v\in\mathcal{P}_{p}(K)$, the
following inverse inequality holds

(39) $\|{\nabla v}\|_{K}^{2}\leq C_{{\rm inv,}K}h_{K}^{-2}\|{v}\|_{K}^{2}.$

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

Next, for a generic $d-$dimensional simplex $T$, we denote its barycentric co-
ordinates 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.

---

Figure 5: Left: a triangle with one ‘small’ face. Right: a polytopic element
with many ‘small’ faces.

$F_{i}$$T$$K$$O$$F_{i}$$\tau_{i}$

Let $K\in\mathcal{T}$ and let $m_{K}$ be the number of its faces. Given that
$K$ is star-shaped by Assumption 2.3, 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
$h_{K}\geq\operatorname{diam}(\tau_{j})\geq r_{K}\geq\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

(40)
$h_{K}^{-1}\leq\|{\nabla\lambda_{j}}\|_{L_{\infty}({\tau_{j}})}=d|\partial\tau_{j}\cap\partial
K|/|\tau_{j}|\leq\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

(41) $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 (41).

######  Lemma 4.3.

For each $K\in{\mathcal{T}}$ satisfying the Assumption 2.3 and for each
$v\in\mathcal{P}_{p}(K)$, we have

(42) $\|{\nabla(b_{K}v)}\|_{{K}}^{2}\leq 2(\tau^{2}+C_{{\rm
inv,}K})h_{K}^{-2}\|{v}\|_{{K}}^{2},$

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

(43) $\|{v}\|_{{K}}^{2}\leq C_{b,K}\|{b_{K}^{1/2}v}\|_{{K}}^{2},$

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

######  Proof 4.4.

Using the triangle inequality, the bound (40), and the inverse inequality
(4.1), we have, respectively,

$\displaystyle\|{\nabla(b_{K}v)}\|_{{K}}^{2}$ $\displaystyle\leq 2\|{(\nabla
b_{K})v}\|_{{K}}^{2}+2\|{b_{K}(\nabla v)}\|_{{K}}^{2}$ $\displaystyle\leq
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}$ (44) $\displaystyle\leq
2\tau^{2}h_{K}^{-2}\sum_{i=1}^{m_{K}}\|{v}\|_{{\tau_{i}}}^{2}+2C_{{\rm
inv,}K}h_{K}^{-2}\|{v}\|_{{K}}^{2}=2(\tau^{2}+C_{{\rm
inv,}K})h_{K}^{-2}\|{v}\|_{{K}}^{2},$

which is the bound required in (42). We now prove the norm equivalence
relation (43). Recalling the norm equivalence relation for each
$v\in\mathcal{P}_{p}(T)$ on a simplex $T$ from [49, Section 3.6], we have

(45)
$\|{v}\|_{{T}}^{2}\leq[2(p+2)]^{2d}\Big{(}\frac{d}{d+1}\Big{)}^{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

(46) $\|{v}\|_{{T}}^{2}\leq C_{b,K}\|{(\lambda_{T}^{i})^{1/2}v}\|_{{T}}^{2},$

for each $i\in\\{0,\dots,d\\}$. Hence, the bound (43) is proven using the
definition of $b_{K}(x)$ in (41):

(47)
$\displaystyle\|{v}\|_{{K}}^{2}=\sum_{j=1}^{m_{K}}\|{v}\|_{{\tau_{j}}}^{2}$
$\displaystyle\leq\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}.$

######  Remark 4.5.

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 (42) and (43) are independent
of the number and measure of the faces of the element $K$.

######  Theorem 4.6 (Elemental residual lower bound).

Let $u$ be the solution of (1) and let $u_{h}\in S_{\mathcal{T}}$ be its dG
approximation under Assumption 2.3 and 2.2. Then, for each
$K\in{\mathcal{T}}$, we have

(48) $\|{h(\Pi f+\nabla\cdot(a\nabla u_{h}))}\|_{{K}}^{2}\leq
2C_{b,K}\big{(}2\big{(}\tau^{2}+C_{{\rm
inv,}K})(\alpha^{*})^{2}\|{\sqrt{a}\nabla e}\|_{{K}}^{2}+O_{K,E}^{2}\big{)}.$

######  Proof 4.7.

We fix $v\in H^{1}_{\rm D}(\Omega)$ as $v|_{K}=b_{K}{(\Pi
f+\nabla\cdot(a\nabla u_{h}))}$, where $b_{K}$ is the element bubble function
in (41), and extended to zero outside $K$. Using relations (30), (3.2), (42),
and $b_{K}\leq 1$, we obtain

$\displaystyle\int_{K}{b_{K}{(\Pi f+\nabla\cdot(a\nabla
u_{h}))}}^{2}\,\mathrm{d}\mathbf{x}=\int_{K}(\Pi
f-f)v\,\mathrm{d}\mathbf{x}+\int_{K}a\nabla e\cdot\nabla
v\,\mathrm{d}\mathbf{x}$ $\displaystyle\leq\|{\Pi
f-f}\|_{{K}}\|{v}\|_{{K}}+\|{\sqrt{a}\nabla e}\|_{{K}}\|{\sqrt{a}\nabla
v}\|_{{K}}$ (49) $\displaystyle\leq(\|{\Pi
f-f}\|_{{K}}+\sqrt{2\big{(}\tau^{2}+C_{{\rm
inv,}K})}\alpha^{*}h_{K}^{-1}\|{\sqrt{a}\nabla e}\|_{{K}}\big{)}\|{(\Pi
f+\nabla\cdot(a\nabla u_{h}))}\|_{{K}},$

Recalling (43), we have $\|{\Pi f+\nabla\cdot(a\nabla u_{h})}\|_{{K}}^{2}\leq
C_{b,K}\|{b_{K}^{1/2}(\Pi f+\nabla\cdot(a\nabla u_{h}))}\|_{{K}}^{2}$ from
which the result (48) already follows.

###  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.8.

The number of faces of every element $K\in\mathcal{T}$ is uniformly bounded.
For $d=3$ only, for every $F\in\Gamma_{\text{\rm int}}$, the radius $r_{F}$ of
the largest $(d-1)$-dimensional ball inscribed in $F$ satisfies
$r_{F}\geq\tau^{-1}h_{F}$, with $\tau$ as in Assumption 2.3.

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.9 (Face bubble).

Let $K\in\mathcal{T}$ and $F\subset\partial K$ a mesh face satisfying
Assumption 4.8. 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.10.

Let $K\in\mathcal{T}$ and $F\subset\partial K$ a mesh face under Assumption
2.3. 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

(50) $\|{v}\|_{{F}}^{2}\leq C_{b,F}\|{b_{F}^{1/2}v}\|_{{F}}^{2},$

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

(51) $\|{\nabla(b_{F}v)}\|_{{T_{F}^{K}}}^{2}\leq C_{{\rm
inv,}{T}}h_{K}h_{F}^{-2}\|{v}\|_{{F}}^{2},$

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

(52) $\|{{\rm curl}\,(b_{T}v)}\|_{{T_{F}^{K}}}^{2}\leq(d-1)C_{{\rm
inv,}{T}}h_{K}h_{F}^{-2}\|{v}\|_{{F}}^{2}.$

######  Proof 4.11.

The bound (50) is given in [35] and (51) follows immediately from (39), the
fact that $b_{F}\leq 1$ and the fact that $\|{v}\|_{{T_{F}^{K}}}^{2}\leq
h_{K}\|{v}\|_{{F}}^{2}$. Finally, the bound (52) is a trivial consequence of
(51) observing once again that $\|{{\rm curl}\,v}\|_{K}^{2}\leq(d-1)\|{\nabla
v}\|_{K}^{2}$.

######  Theorem 4.12 (flux residuals lower bound).

Let $u$ be the solution of (1) and let $u_{h}\in S_{\mathcal{T}}$ be its dG
approximation under Assumption 2.3, 2.2 and 4.8. Then, for each
$K\in{\mathcal{T}}$, we have

(53) $\displaystyle\|{{h}^{1/2}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{\partial K\cap\Gamma_{\text{\rm
int}}}}^{2}$ $\displaystyle\qquad\leq
6C_{b,F}\Big{(}2C_{b,K}\big{(}\tau^{2}+C_{{\rm
inv,}K})(\alpha^{*})^{2}\|{\sqrt{a}\nabla
e}\|_{{\omega_{K}}}^{2}+\sum_{K^{\prime}\in\omega_{K}}(1+2C_{b,K})O_{K^{\prime},E}^{2}$
$\displaystyle\qquad\qquad\qquad+\alpha^{*}C_{{\rm inv,}{T}}\sum_{F\in\partial
K\cap\Gamma_{\text{\rm int}}}h_{F}^{-{1}}h_{F^{\perp}}\|{\sqrt{a}\nabla
e}\|_{{\omega_{F}}}^{2}\Big{)},$

and

(54) $\|{h^{1/2}[\\![\nabla_{{\rm T}}u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}\leq 2C_{b,F}(d-1)\alpha_{*}C_{{\rm inv,}{T}}\sum_{F\in\partial
K\cap\Gamma_{\text{\rm 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^{\prime}}$ with $K^{\prime}$ the element
neighbouring $K$ across $F$ and
$h_{F^{\perp}}:=\max\\{h_{K},h_{K^{\prime}}\\}$.

######  Proof 4.13.

In view of proving (53), we first consider any $F\in\partial K$ with
$F\in\Gamma_{\text{\rm int}}$. Further, we fix $v\in H^{1}(\omega_{F})$ as the
constant extension of $[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]|_{F}$ in the direction normal to $F$,
so that $b_{F}v\in H^{1}_{0}(\omega_{F})$, cf. Lemma 4.10. Then, testing the
error equation (30) with $b_{F}v$ extended to zero on the whole of $\Omega$,
we get

$\displaystyle\int_{\omega_{F}}a\nabla_{h}e\cdot\nabla(b_{F}v)\,\mathrm{d}\mathbf{x}=$
$\displaystyle\int_{\omega_{F}}(f-\Pi
f)(b_{F}v)\,\mathrm{d}\mathbf{x}+\int_{\omega_{F}}{(\Pi f+\nabla\cdot(a\nabla
u_{h}))}(b_{F}v)\,\mathrm{d}\mathbf{x}$ $\displaystyle-\int_{F}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]^{2}b_{F}\,\mathrm{d}{s}.$

From this, using (51) and the fact that,
$\|{b_{F}v}\|_{{K}}^{2}=\|{b_{F}v}\|_{{T_{F}^{K}}}^{2}\leq
h_{K}\|{[\\![a\nabla u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{F}}^{2}$,
the same bound being true on $K^{\prime}$, we obtain

$\displaystyle\|{[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]b_{F}^{1/2}}\|_{{F}}^{2}\leq$
$\displaystyle\sum_{\mathcal{K}\in\\{K,K^{\prime}\\}}\left[(\|{\Pi
f-f}\|_{{T_{F}^{\mathcal{K}}}}+\|{(\Pi f+\nabla\cdot(a\nabla
u_{h}))}\|_{{T_{F}^{\mathcal{K}}}})\|{b_{F}v}\|_{{T_{F}^{\mathcal{K}}}}\right.$
$\displaystyle\left.\quad\quad\qquad+\|{\sqrt{a}\nabla
e}\|_{{T_{F}^{\mathcal{K}}}}\|{\sqrt{a}\nabla(b_{F}v)}\|_{{T_{F}^{\mathcal{K}}}}\right]$
$\displaystyle\leq$
$\displaystyle\sum_{\mathcal{K}\in\\{K,K^{\prime}\\}}\left[h_{\mathcal{K}}^{1/2}(\|{\Pi
f-f}\|_{{T_{F}^{\mathcal{K}}}}+\|{(\Pi f+\nabla\cdot(a\nabla
u_{h}))}\|_{{T_{F}^{\mathcal{K}}}})\right.$
$\displaystyle\quad\quad\left.+(\alpha^{*}C_{{\rm
inv,}{T}}h_{F}^{-1}h_{\mathcal{K}})^{1/2}h_{F}^{-1/2}\|{\sqrt{a}\nabla
e}\|_{{T_{F}^{\mathcal{K}}}}\right]\|{[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{F}}.$

This, together with (50), gives

$\displaystyle\|{h^{1/2}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{F}}^{2}\leq$ $\displaystyle
2C_{b,F}\\!\\!\sum_{\mathcal{K}\in\\{K,K^{\prime}\\}}\left[(\|{h(\Pi
f-f)}\|_{{T_{F}^{\mathcal{K}}}}+\|{h(\Pi f+\nabla\cdot(a\nabla
u_{h}))}\|_{{T_{F}^{\mathcal{K}}}})\right.$
$\displaystyle\qquad\quad\quad\left.+(\alpha^{*}C_{{\rm
inv,}{T}})^{1/2}h_{F}^{-1}h_{\mathcal{K}}\|{\sqrt{a}\nabla
e}\|_{{T_{F}^{\mathcal{K}}}}\right]^{2}.$

Summing over all internal faces of $K$, noting carefully that the involved
domains do not overlap, we finally obtain

$\displaystyle\|{{h}^{1/2}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{\partial K\cap\Gamma_{\text{\rm
int}}}}^{2}\leq$ $\displaystyle 6C_{b,F}\Big{(}\|{h(\Pi
f-f)}\|_{{\omega_{K}}}^{2}+\|{h(\Pi f+\nabla\cdot(a\nabla
u_{h}))}\|_{{\omega_{K}}}^{2}$ (55) $\displaystyle+\alpha^{*}C_{{\rm
inv,}{T}}\sum_{F\in\partial K\cap\Gamma_{\text{\rm
int}}}(h_{F}^{-{1}}h_{F^{\perp}})^{2}\|{\sqrt{a}\nabla
e}\|_{{\omega_{F}}}^{2}\Big{)},$

as $h_{F^{\perp}}:=\max\\{h_{K},h_{K^{\prime}}\\}$, from which (53) now
follows by employing (48).

The proof concerning the tangential jump residual is similar. Given
$F\in\Gamma_{\text{\rm int}}$, we fix $v\in H^{1}(\omega_{F})$ as the constant
extension of $[\\![\nabla_{{\rm T}}u_{h}]\\!]|_{F}$ in the direction normal to
$F$, so that $b_{F}v\in H^{1}_{0}(\omega_{F})$. Using the fact that
$\text{curl}\,\nabla u=\mbox{\boldmath$\rm{0}$}$, we have the key observation

(56) $\int_{\omega_{F}}\text{curl}\,(b_{F}v)\cdot\nabla
u\,\mathrm{d}\mathbf{x}=0.$

Integration by parts and (52), give

$\displaystyle\|{[\\![\nabla_{{\rm T}}u_{h}]\\!]b_{F}^{1/2}}\|_{{F}}^{2}=$
$\displaystyle\int_{\omega_{F}}\text{curl}\,(b_{F}v)\cdot\nabla
u_{h}\,\mathrm{d}\mathbf{x}=\int_{\omega_{F}}\text{curl}\,(b_{F}v)\cdot\nabla(u_{h}-u)\,\mathrm{d}\mathbf{x}$
$\displaystyle\leq\alpha_{*}\sum_{\mathcal{K}=K,K^{\prime}}\|{\sqrt{a}\nabla
e}\|_{{T_{F}^{\mathcal{K}}}}\|{\text{curl}\,(b_{T}v)}\|_{{T_{F}^{\mathcal{K}}}}$
$\displaystyle\leq\alpha_{*}\big{(}(d-1)C_{{\rm
inv,}{T}}h_{\mathcal{K}}\big{)}^{1/2}h_{F}^{-1}\|{[\\![\nabla_{{\rm
T}}u_{h}]\\!]}\|_{{F}}\sum_{\mathcal{K}\in\\{K,K^{\prime}\\}}\|{\sqrt{a}\nabla
e}\|_{{T_{F}^{\mathcal{K}}}}.$

Hence, by using (50), we obtain

$\|{h^{1/2}[\\![\nabla_{{\rm T}}u_{h}]\\!]}\|_{{F}}^{2}\leq
2\alpha_{*}(d-1)C_{{\rm
inv,}{T}}C_{b,F}(h_{F}^{-1}h_{F^{\perp}})^{2}\|{\sqrt{a}\nabla
e}\|_{{\omega_{F}}}^{2}.$

The required lower bound on the jump of the tangential gradient now follows by
summing over all internal faces belonging to $K$.

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.12
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.9 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^{\prime}}}$. 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.14.

Under the assumption of Theorem 4.12, if, moreover, for every
$K\in\mathcal{T}$ and $F\in\partial K$ it holds $h_{F}\geq\tau^{-1}h_{K}$,
then

(57) $\displaystyle\|{{h}^{1/2}[\\![a\nabla
u_{h}\cdot\mbox{\boldmath$\rm{n}$}]\\!]}\|_{{\partial K\cap\Gamma_{\text{\rm
int}}}}^{2}$ $\displaystyle\qquad\leq
6C_{b,F}\Big{(}(2C_{b,K}\big{(}\tau^{2}+C_{{\rm
inv,}K})(\alpha^{*})^{2}+\tau\alpha^{*}C_{{\rm inv,}K})\|{\sqrt{a}\nabla
e}\|_{{\omega_{K}}}^{2}$
$\displaystyle\qquad\qquad\qquad+\sum_{K^{\prime}\in\omega_{K}}(1+2C_{b,K})O_{K^{\prime},E}^{2}\Big{)},$

and

(58) $\|{h^{1/2}[\\![\nabla_{{\rm T}}u_{h}]\\!]}\|_{{\Gamma_{\text{\rm
int}}}}^{2}\leq 2C_{b,F}(d-1)\tau\alpha_{*}C_{{\rm inv,}K}\|{\sqrt{a}\nabla
e}\|_{{\omega_{K}}}^{2}.$

######  Remark 4.15.

Corollary 4.14 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.14 together with Theorem 4.6 and Theorem 3.6 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_{\mathcal{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}\leq C_{{\rm
inv,}{T}}h_{K}^{-1}\|{v}\|_{{F}}^{2},$

holds true, the lower bounds of the flux residuals (53) and (54) 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 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.

|
---|---
Figure 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.

We consider the problem (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 $\mathcal{O}(DoFs^{-\frac{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.8 may not be necessary.

|
---|---
Figure 7: Example 1. Convergence history of the error and error estimator of
Equation 38 for $p=1,2,3,4$ (left) and respective effectivity measured as the
ratio estimator over error (right).

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}$.

| $p=1$ | $p=2$
---|---|---
$\\#$ elem | $R_{E}$ | $R_{J}$ | $R_{N}$ | $R_{T}$ | $R_{E}$ | $R_{J}$ | $R_{N}$ | $R_{T}$
$114$ | $64\%$ | $18\%$ | $9\%$ | $9\%$ | $56\%$ | $26\%$ | $9\%$ | $9\%$
$498$ | $59\%$ | $23\%$ | $9\%$ | $9\%$ | $45\%$ | $36\%$ | $9\%$ | $10\%$
$2063$ | $36\%$ | $37\%$ | $13\%$ | $14\%$ | $45\%$ | $33\%$ | $10\%$ | $12\%$
$8912$ | $32\%$ | $37\%$ | $14\%$ | $17\%$ | $43\%$ | $31\%$ | $11\%$ | $15\%$
Table 1: 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}$.

Further, to highlight the importance of 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.5, 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}\leq 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
(38) 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 8: Example 1. Convergence history of the error, error estimator of
Equation 38 and the classical estimator for $p=1,2$.

###  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 (38)
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:

$\text{SOLVE}\longrightarrow\text{ESTIMATE}\longrightarrow\text{MARK}\longrightarrow\text{REFINE/AGGLOMERATE}.$

We consider the problem (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

$\displaystyle u=$ $\displaystyle
r^{2/3}\sin(2\psi/3)+\exp(-1000((x-0.5)^{2}+(y-0.25)^{2}))$
$\displaystyle+\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
$\mathcal{O}(DoFs^{-\frac{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.

|  |
---|---|---
|  |
Figure 9: Error, estimator (top) and effectivity index (bottom) of the new adaptive polytopic dG method with $p=1,2,3$. |  |
---|---|---
|  |
Figure 10: The sequence of meshes generated by the adaptive polytopic dG algorithm with $p=3$. |
---|---
Figure 11: Example 2. Error and estimator of the new adaptive polytopic dG
method and adaptive dG method without agglomeration for $p=2,3$.

## 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.
  * [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äßle, 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.
  * [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 ncvem 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. Veeser, 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. Veeser 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.