ADAPTIVE ROBUST INTERIOR PENALTY
DISCONTINUOUS GALERKIN METHODS

ZHAONAN DONG, EMMANUIL H. GEORGOULIS, AND NATALIA KOPTEVA

ABSTRACT. This work is concerned with the development of an adaptive algo-
rithm for the so-called robust interior penalty discontinuous Galerkin (RIPDG)
method for an elliptic model problem, based on new rigorous *a posteriori* error
bounds. The RIPDG methodology offers a combined robust behaviour with
respect to sharp local variation in mesh-size, polynomial degree and diffusion
contrast simultaneously. Thus, special attention is given on the dependence
of the resulting *a posteriori* error bounds in these quantities. The proof em-
ploys a special $H^1$- conforming reconstruction of the numerical solution, along
with a Helmholtz decomposition argument. A series of numerical experiments
showcasing the reliability of the proven *a posteriori* error bounds are included.

1. INTRODUCTION

Discontinuous Galerkin (DG) methods are now established as a general-purpose
methodology for the numerical approximation of partial differential equation (PDE)
problems of various types. A significantly attractive attribute of DG methods is
the extreme flexibility in the mesh design they admit (see, e.g., [14, 9, 10, 13, 11]
and the references therein). This is particularly pertinent in the context of adaptive
algorithms.

At the same time, the proof of *a posteriori* error bounds that drive mathematically justifiable adaptive methodologies is typically considerably more technical for discontinuous Galerkin methods than their conforming counterparts. The derivation of various *a posteriori* error bounds and the convergence of adaptive algorithms for discontinuous Galerkin methods for elliptic problems has received considerable interest during the last 20 years or so; see [29, 30, 27, 1, 2, 22, 6, 12, 33] and other works. A key technical tool for the proof of such bounds is the reconstruction/recovery of a piecewise discontinuous numerical solution into the natural (conforming) function space of the exact PDE problem solution. In the celebrated work of Karakashian & Pascal [29], the DG solution is recovered into a globally continuous, element-wise polynomial via nodal averaging. Later works involve to the use of Helmholtz decomposition [1, 2] or, more recently, via ‘orthogonal’ $H^1$-reconstructions [20, 12, 17].

With the exception of [22, 33], the above developments concern the classical
interior penalty DG (IPDG) method for elliptic problems of Arnold [4]. Ern &
Stephansen [22] and Kopteva & Rankin [33] consider IPDG methods involving
weighted averages in the respective numerical fluxes to address the issues of high
contrast diffusion and of high local variation in meshsize, respectively. The, so-
called, *robust interior penalty discontinuous Galerkin* (RIPDG) method for elliptic
problems, introduced in [18] is a variant of the classical IPDG admitting robust
behaviour with respect to sharp local variation in mesh-size, polynomial degree and
diffusion contrast simultaneously, by a careful selection of weighted averages in the
definition of the respective numerical fluxes. RIPDG unifies and extents previous
works on the development of IPDG methods with weighted averages [23, 3, 25],
which have been, in turn, motivated by developments in domain decomposition [21, 8]; see also the recent developments in [33].

This work is concerned with the proof of residual-type a posteriori error bounds in the energy norm for the hp-version robust IPDG method. The analysis is based on a variant of the ‘orthogonal’ H¹-reconstruction from [20, 12, 17], along with a new localised estimation of the diffusion tensor. The analysis tracks carefully the local variations in meshsize, polynomial degree and diffusion contrast. Moreover, in an important deviation from available a posteriori error analyses for IPDG methods, we also discuss the case of ‘rough’, i.e., essentially bounded from above and below, diffusion tensors. Indeed, to the best of our knowledge, the combined treatment of, possibly highly, locally variable rough diffusion coefficients, rendering natural the employment of the robust version of IPDG, has not been addressed in the literature.

For simplicity, we focus on the basic second order elliptic boundary value problem on an open bounded polyhedral Lipschitz domain $\Omega \subset \mathbb{R}^d$, $d \in \mathbb{N}$: find $u \in H^1(\Omega)$ such that

$$ (1.1) \qquad \begin{aligned} -\nabla \cdot (a\nabla u) + cu &= f && \text{in } \Omega; \\ u &= g && \text{on } \partial\Omega, \end{aligned} $$

with data $f \in L_2(\Omega)$, $g \in H^{\frac{1}{2}}(\partial\Omega)$ and a positive definite diffusion tensor $a \in [L_\infty(\Omega)]^{d \times d}$, with $c \in L_\infty(\Omega)$ non-negative function. The problem (1.1) is understood in the weak sense. Considerable extensions to the above problem are possible, but we refrain from exhausting full generality in the interest of clarity of the presentation.

The remainder of this work is structured as follows. The numerical methods are presented in Section 2, along with proofs of coercivity and continuity. Section 3 is concerned with the proof of a posteriori error bounds for the reconstructed/conforming part of the error. The analysis is completed in Section 4 with the proof of computable upper bound for the non-conforming remainder of the error. Finally in Section 5 we present some numerical investigations.

## 2. INTERIOR PENALTY DISCONTINUOUS GALERKIN METHODS

We consider a family of shape-regular meshes. Each mesh $\mathcal{T}$ consists of mutually disjoint open simplicial elements $K \in \mathcal{T}$, so that $\cup_{K \in \mathcal{T}} \bar{K} = \bar{\Omega}$. Let also $h_K := \text{diam}(K)$ the diameter of $K \in \mathcal{T}$, and define the mesh-function $h: \cup_{K \in \mathcal{T}} K \to \mathbb{R}_+$ by $h|_K = h_K$, $K \in \mathcal{T}$. We let $\omega_1(K)$ denote the collection of elements $\bar{K} \in \mathcal{T}$ sharing at least one vertex with the element $K$. Similarly, $\omega_k(K)$ denotes the collection of cells $\bar{K} \in \mathcal{T}$ sharing at least one vertex with $\omega_{k-1}(K)$.

The set of mesh vertices is denoted by $\mathcal{V}_h$ and is decomposed into the subset of interior vertices, $\mathcal{V}_h^i$, and the subset of boundary vertices, $\mathcal{V}_h^b$. For all $a \in \mathcal{V}_h$, $\mathcal{T}_a$ denotes the collection of mesh cells which share $a$ and $\omega_a$ the corresponding open subdomain. The set of mesh faces is denoted by collection of faces in $\mathcal{F}_h$. In addition, we define $\mathcal{F}_a$ as the collection of faces in $\mathcal{F}_h$ which share $a$.

Also, let $\Gamma := \cup_{K \in \mathcal{T}} \partial K$ denote the mesh skeleton and let $\Gamma_{\text{int}} := \Gamma \setminus \partial \Omega$. The mesh skeleton $\Gamma$ is decomposed into $(d-1)$-dimensional simplicial faces $F$, shared by at most two elements. Let also $\nabla_h$ denote the broken gradient defined as $\nabla_h v|_K := \nabla(v|_K)$, $K \in \mathcal{T}$.

The hp-discontinuous Galerkin space $S_{\mathcal{T}}^{\mathbf{P}}$ with respect to $\mathcal{T}$ is defined as

$$ S_{\mathcal{T}}^{\mathbf{P}} := \{v \in L_2(\Omega) : v|_K \in \mathcal{P}_{p_K}(K), K \in \mathcal{T}\}, $$

with $\mathcal{P}_{p_K}(K)$ denoting the space of $d$-variate polynomials of total degree up to $p_K$ on $K$; we collect all polynomial degrees in a vector $\mathbf{p} := \{p_K : K \in \mathcal{T}\}$.
Let $K_+$ and $K_-$ be two adjacent elements of $\mathcal{T}$ sharing a face $F \subset \partial K_+ \cap \partial K_- \subset \Gamma_{\text{int}}$. For $v$ and $q$ element-wise continuous scalar- and vector-valued functions, respectively, we define the weighted average across $F$ by

$$ \{v\}_{w_F}|_F := w_+v_+|_F + w_-v_-|_F, \quad \{q\}_{w_F}|_F := w_+q_+|_F + w_-q_-|_F, $$

for $w_F := (w_+, w_-) \in \mathbb{R}_+^2$ with $w_+ + w_- = 1$, respectively, with $v_{\pm}|_F$ denoting the trace of $v$ from the element $K_{\pm}$ on $F$, and correspondingly for $q$. Also, we define the jump across $F$ by

$$ [v]|_F := v_+ - v_-, \quad \{q\}|_F := q_+ - q_-, $$

using the convention that element numbering of $v_+$ is greater than that of $v_-$ to determine the sign. In addition, $n_F$ denoting the unit outward normal vector on $F$ pointing from $v_+$ to $v_-$. We collect all weight pairs in a weight function $\mathbf{w} : \Gamma_{\text{int}} \to \mathbb{R}_+^2$, with $\mathbf{w}|_F = w_F$, for each face $F \subset \Gamma_{\text{int}}$. On a boundary face $F \subset \partial\Omega \cap \partial K$, we set $\{v\} := v$, $\{q\} := q$, $[v] := v$, and $\{q\} := q$, respectively. In view of the latter, we extend $\mathbf{w}$ to $\Gamma$ by setting $\mathbf{w}|_{\partial\Omega} = (1,0)$ with $K_- = \emptyset$ there.

Assuming $a \in [W^{1,\infty}(\Omega, \mathcal{T})]^{d \times d}$, with $W^{1,\infty}(\Omega, \mathcal{T})$ denoting the 'broken' space, containing functions which are element-wise in $W^{1,\infty}(K)$, $K \in \mathcal{T}$, the (weighted) interior penalty discontinuous Galerkin (IPDG) method reads: find $u_h \in S_{\mathcal{T}}^{\mathcal{P}}$, such that

$$ (2.1) \qquad B(u_h, v_h) = \ell(v_h), \quad \text{for all } v_h \in S_{\mathcal{T}}^{\mathcal{P}}, $$

with

$$ B(u_h, v_h) := \int_{\Omega} (a \nabla_h u_h \cdot \nabla_h v_h + c u_h v_h) \, dx + \int_{\Gamma} \sigma \{[u_h][v_h]\} \, ds \\ - \int_{\Gamma} (\{\{a \nabla u_h\} \cdot \mathbf{n}\}_{\mathbf{w}}[v_h] + \theta \{\{a \nabla v_h\} \cdot \mathbf{n}\}_{\mathbf{w}}[u_h]\} \, ds, $$

and

$$ \ell(v_h) := \int_{\Omega} f v_h \, dx + \int_{\partial\Omega} g(\sigma v_h - \theta(a \nabla v_h) \cdot \mathbf{n}) \, ds, $$

with $\theta \in [-1, 1]$ and $\sigma : \Gamma \to \mathbb{R}$ the, so-called, discontinuity-penalization or penalty function, whose precise definition, along with that of the weights $\mathbf{w}$ is a central concern in this work.

The choice $\theta = 1$ in (2.1) yields the symmetric version, while the choice $\theta = -1$ the non-symmetric version of the weighted IPDG method.

The selection $w_F = (1/2, 1/2)$, for all $F \subset \Gamma_{\text{int}}$, corresponds to the classical IPDG method of Wheeler [36] and of Arnold [4]. The case $w_F = (1, 0)$ or $w_F = (0, 1)$ for all $F \subset \Gamma_{\text{int}}$ gives the original IPDG method of Baker [5]. Also, in [21, 8], the weighted average choice $w_F = (\delta_+, \delta_-)$ with

$$ \delta_{\pm} := \frac{\alpha_{\mp}}{\alpha_{+} + \alpha_{-}}, $$

with $\alpha_{\pm} := a|_{K_{\pm}}$, was proposed to ensure robustness with respect to variation of element-wise constant scalar diffusion coefficients, while in [23] the selection $\alpha_{\pm} := n_F^T a|_{K_{\pm}} n_F$ for element-wise constant diffusion tensors was provided, upon combining the selections in [21, 8] and [24]. Also, weighted averages involving also local mesh-sizes have been proposed for treating interface problems [3, 25] on non-matching grids on either side of the interface; see also [33]. Finally, in [19], the, so-called, robust interior penalty DG methods is presented admitting a choice of $\sigma$ yielding robust dependence not only with respect to the local variation of the diffusion, but also with respect to the local direction-wise variations in mesh-size and polynomial degree.
**Remark 2.1.** The above IPDG method (2.1) is well defined for $a \in [W^{1,\infty}(\Omega, \mathcal{T})]^{d \times d}$.
The definition of the method is unclear if $a \in [L_{\infty}(\Omega)]^{d \times d}$, as the face integrals involving the normal fluxes, $\{\!\{a\nabla w_h\}\!\}$, are not necessarily integrable on the skeleton $\Gamma$. Typical remedies involve the local projection of the diffusion tensor to element-wise constant functions or, polynomial, in general functions at the cost of introducing an additional consistency error.

To ensure the bilinear form $B: S_{\mathcal{T}}^{\mathbf{P}} \times S_{\mathcal{T}}^{\mathbf{P}} \to \mathbb{R}$ accepts functions in $H^1(\Omega)$ as arguments, we define its extension, still denoted by $B: S \times S \to \mathbb{R}$, for $S := H^1(\Omega) + S_{\mathcal{T}}^{\mathbf{P}}$, given by

$$B(w, v) := \int_{\Omega} (a\nabla_h w \cdot \nabla_h v + cwv) \, dx + \int_{\Gamma} \sigma\!\left[\!\left[w\right]\!\right]\!\left[\!\left[v\right]\!\right] \, ds \\ - \int_{\Gamma} \left(\{\!\{a\Pi_{\mathbf{p}-1}(\nabla w) \cdot \mathbf{n}\}\!\}_{\mathbf{w}}\!\left[\!\left[v\right]\!\right] + \theta\{\!\{a\Pi_{\mathbf{p}-1}(\nabla v) \cdot \mathbf{n}\}\!\}_{\mathbf{w}}\!\left[\!\left[w\right]\!\right]\right) \, ds,$$

for $w, v \in S$, with $\Pi_{\mathbf{p}-1}: L_2(\Omega) \to S_{\mathcal{T}}^{\mathbf{P}-1}$ denoting the component-wise orthogonal $L_2$-projection onto $S_{\mathcal{T}}^{\mathbf{P}-1}$, i.e., the space of one degree less elementwise than the original approximation space. Note that $B$ collapses to the original bilinear form in (2.1) when restricted to $S_{\mathcal{T}}^{\mathbf{P}} \times S_{\mathcal{T}}^{\mathbf{P}}$, i.e., it is an extension. For completeness, we also set

$$\ell(v) := \int_{\Omega} fv \, dx + \int_{\partial\Omega} g(\sigma v - \theta a \Pi_{\mathbf{p}-1}(\nabla v) \cdot \mathbf{n}) \, ds,$$

noting that $B(u, v) = \ell(v)$, for all $v \in H^1(\Omega)$, with $u \in H^1(\Omega)$ denoting the exact solution to the problem (1.1).

We now discuss the coercivity and continuity of the bilinear form $B(\cdot, \cdot)$ with respect to the norm

$$\|w\| := \left( \|\sqrt{a\nabla_h w}\|_{\Omega}^2 + \|\sqrt{cw}\|_{\Omega}^2 + \|\sqrt{\sigma}\!\left[\!\left[w\right]\!\right]\!\|_{\Gamma}^2 \right)^{\frac{1}{2}}.$$

We have

$$B(v, v) = \|\sqrt{a\nabla_h v}\|_{\Omega}^2 + \|\sqrt{cv}\|_{\Omega}^2 + \|\sqrt{\sigma}\!\left[\!\left[v\right]\!\right]\!\|_{\Gamma}^2 \\ - (1 + \theta) \int_{\Gamma} \{\!\{a\Pi_{\mathbf{p}-1}(\nabla v) \cdot \mathbf{n}\}\!\}_{\mathbf{w}}\!\left[\!\left[v\right]\!\right] \, ds.$$

Now, for a generic face $F$ shared by two elements $K_+$ and $K_-$, we have

$$- \int_{\Gamma} \{\!\{a\Pi_{\mathbf{p}-1}(\nabla v) \cdot \mathbf{n}\}\!\}_{\mathbf{w}}\!\left[\!\left[v\right]\!\right] \, ds \\ \geq - \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} w_* \|a|_{K_*} \mathbf{n}\|_{L_{\infty}(F)} \|\Pi_{p_{K_*}-1}(\nabla v)|_{K_*}\|_F \right) \|\!\left[\!\left[v\right]\!\right]\!\|_F \\ \geq - \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} w_* \|a|_{K_*} \mathbf{n}\|_{L_{\infty}(F)} \sqrt{\frac{p_{K_*}(p_{K_*} + d - 1)|F|}{d|K_*|}} \|\nabla v\|_{K_*} \right) \|\!\left[\!\left[v\right]\!\right]\!\|_F,$

using the trace inverse estimate $\|v\|_F^2 \leq (p+1)(p+d)|F|/(d|K|)\|v\|_K^2$ for any $v \in \mathcal{P}_p(K)$ [35]. Setting, now,

$$\zeta_* := \left( (d+1)\|a|_{K_*} \mathbf{n}\|_{L_{\infty}(F)} \|a^{-1/2}\|_{L_{\infty}(K_*)} \sqrt{\frac{p_{K_*}(p_{K_*} + d - 1)|F|}{d|K_*|}} \right)^{-1},$$
with $* \in \{+, -\}$, for brevity, the above gives

$$
\begin{align*}
& \left| \int_{\Gamma} \{ \|a\Pi_{\mathbf{p}-1}(\nabla v) \cdot \mathbf{n}\|_{\mathbf{w}} \|v\| \} ds \right| \\
& \le \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} \frac{w_*}{(d+1)\zeta_*} \|\sqrt{a\nabla v}\|_{K_*} \right) \|\|v\|\|_F \\
& \le \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} \frac{1}{d+1} \|\sqrt{a\nabla v}\|_{K_*} \right) \max_{* \in \{+, -\}} \frac{w_*}{\zeta_*} \|\|v\|\|_F \\
(2.2) \qquad & \le \left( \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} \frac{1}{d+1} \|\sqrt{a\nabla v}\|_{K_*} \right)^2 \right)^{\frac{1}{2}} \left( \sum_{F \subset \Gamma} \max_{* \in \{+, -\}} \frac{w_*^2}{\zeta_*^2} \|\|v\|\|_F \right)^{\frac{1}{2}} \\
& \le \left( \sum_{F \subset \Gamma_*} \sum_{* \in \{+, -\}} \frac{2}{(d+1)^2} \|\sqrt{a\nabla v}\|_{K_*}^2 \right)^{\frac{1}{2}} \left( \sum_{F \subset \Gamma} \max_{* \in \{+, -\}} \frac{w_*^2}{\zeta_*^2} \|\|v\|\|_F \right)^{\frac{1}{2}} \\
& \le 2 \|\sqrt{a\nabla v}\|_{\Omega} \left( \sum_{F \subset \Gamma} \max_{* \in \{+, -\}} \frac{w_*^2}{\zeta_*^2} \|\|v\|\|_F \right)^{\frac{1}{2}} \\
& \le \frac{1}{3+\theta} \|\sqrt{a\nabla v}\|_{\Omega}^2 + (3+\theta) \left( \sum_{F \subset \Gamma} \max_{* \in \{+, -\}} \frac{w_*^2}{\zeta_*^2} \|\|v\|\|_F \right)^{\frac{1}{2}},
\end{align*}
$$

since each element has exactly $d+1$ ($d-1$)-dimensional faces.
Selecting, now,

$$
\tau|_{F} = (1 + \theta)(3 + \theta) \max_{* \in \{+, -\}} \frac{w_{*}^{2}}{\zeta_{*}^{2}},
$$

on every face $F \subset \Gamma$, we deduce

$$
B(v, v) \geq \frac{2}{3+\theta} \|\sqrt{a\nabla_h v}\|_{\Omega}^2 + \|\sqrt{cv}\|_{\Omega}^2 + \|\sqrt{\sigma - \tau}\|_{v}\|_{\Gamma}^2.
$$

A good choice for $\sigma$, as we shall see below, satisfies the equation $(1+\theta)\sigma = (3+\theta)\tau$, thereby dictating the selection

$$
\sigma|_F := (3 + \theta)^2 \max_{* \in \{+, -\}} \frac{w_*^2}{\zeta_*^2},
$$

for each $F \subset \Gamma$, whenever $\theta \in (-1, 1]$, and can be selected as an arbitrary non-negative number $\sigma \ge 0$, whenever $\theta = -1$. With the above choice of $\sigma$, we deduce the coercivity of the bilinear form:

$$
(2.3) \quad B(v,v) \ge \frac{2}{3+\theta} \left( \|\sqrt{a\nabla_h v}\|_{\Omega}^2 + \|\sqrt{\sigma\|v\|}_{\Gamma}\|^2 + \|\sqrt{cv}\|_{\Omega}^2 \ge \frac{2}{3+\theta} \|v\|^2 \right).
$$

The continuity of the bilinear form $B(\cdot, \cdot)$ follows using similar arguments. To this end, working as in (2.2), we have, for $\theta \in (-1, 1]$,

$$
\begin{align}
B(w, v) &\le \|w\|\|v\| + \frac{2}{3+\theta} \left( \|\sqrt{a\nabla w}\|_{\Omega} \|\sqrt{\sigma}\|_F + \|\sqrt{a\nabla v}\|_{\Omega} \|\sqrt{\sigma}\|_F \right) \nonumber \\
(2.4) \qquad &\le \frac{5+\theta}{3+\theta} \|w\|\|v\|. \nonumber
\end{align}
$$

For $\theta = -1$, we have, instead

$$
(2.5) \qquad B(w, v) \le (1 + C_{\text{NIP}}) \|w\| \|v\|,
$$

with $C_{\text{NIP}} := 2 \| \max_{* \in \{+, -\}} w_*(\zeta_* \sigma)^{-1} \|_{L_\infty(\Gamma)}$, for any $\sigma : \Gamma \to \mathbb{R}$ non-negative function.
We now comment on particular cases. For the classical IP(DG) for $\theta \in (-1, 1]$ [36, 4, 28], we have $w_* = \frac{1}{2}$, thereby we deduce

$$ \sigma_{\text{IP}} = \left(\frac{3+\theta}{4}\right)^2 \max_{* \in \{+, -\}} \zeta_*^{-2}, $$

while for the Robust IP(DG) [19], we select $w_* = \frac{\zeta_*}{\zeta_+ + \zeta_-}$, thereby, deducing

$$ \sigma_{\text{RIP}} = (3 + \theta)^2 (\zeta_+ + \zeta_-)^{-2}. $$

The terminology ‘robust’ stems from the fact that $\sigma_{\text{RIP}}$ remains bounded under extreme local variations of diffusion coefficient, mesh-size, polynomial degree; we refer to [19] for a detailed discussion. We note that the above analysis extends slightly (and corrects a missing factor of 4) the respective developments in [19].

We now return to the question of definition of an IPDG method for elliptic problems for which the diffusion tensor is just assumed to be essentially bounded, viz., $a \in [L_\infty(\Omega)]^{d \times d}$, and positive definite. Instead of projecting $a$ into a piecewise constant or polynomial function, possibly subordinate to the mesh, we instead consider the weighted IPDG variant reading: find $u_h \in S_T^p$, such that

$$ (2.6) \qquad \mathcal{B}(u_h, v_h) = \iota(v_h), \quad \text{for all } v_h \in S_T^p, $$

with

$$ \mathcal{B}(u_h, v_h) := \int_{\Omega} (a \nabla_h u_h \cdot \nabla_h v_h + c u_h v_h) \, dx + \int_{\Gamma} \sigma \llbracket u_h \rrbracket \llbracket v_h \rrbracket \, ds \\ - \int_{\Gamma} (\{\Pi_{\mathbf{p}-1}(a \nabla u_h) \cdot \mathbf{n}\} \mathbf{w} \llbracket v_h \rrbracket + \theta \{\Pi_{\mathbf{p}-1}(a \nabla v_h) \cdot \mathbf{n}\} \mathbf{w} \llbracket u_h \rrbracket) \, ds, $$

and

$$ \iota(v_h) := \int_{\Omega} f v_h \, dx + \int_{\partial \Omega} g(\sigma v_h - \theta \Pi_{\mathbf{p}-1}(a \nabla v_h) \cdot \mathbf{n}) \, ds. $$

The key difference between (2.6) and the (R)IPDG method (2.1) is the application of orthogonal $L_2$-projection onto the finite element space for the whole numerical fluxes. In practice, this requires to compute the orthogonal $L_2$ projection of $a \nabla \phi_i$ for each local basis function $\phi_i$. This is a local operation, made easier when using elemental $L_2$-orthogonal bases, affecting only the assembly step, which is fully parallelisable, e.g., on a GPU chip. We note that the complexity of performing the projections is entirely comparable to first projecting $a$ and then applying the (R)IPDG method (2.1).

To show coercivity and continuity properties, we work similarly to above. In particular, we have

$$ \mathcal{B}(v, v) = \|\sqrt{a} \nabla_h v\|_{\Omega}^2 + \|\sqrt{c} v\|_{\Omega}^2 + \|\sqrt{\sigma} \llbracket v \rrbracket\|_{\Gamma}^2 - (1+\theta) \int_{\Gamma} \{\Pi_{\mathbf{p}-1}(a \nabla v) \cdot \mathbf{n}\} \mathbf{w} \llbracket v \rrbracket \, ds, $$

so that, for a generic face $F$ shared by two elements $K_+$ and $K_-$, we get

$$ \begin{aligned} & - \int_{\Gamma} \{\Pi_{\mathbf{p}-1}(a \nabla v) \cdot \mathbf{n}\} \mathbf{w} \llbracket v \rrbracket \, ds \\ & \geq - \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} w_* \|\Pi_{p_{K_*}-1}(a \nabla v)|_{K_*} \|_F \right) \|\llbracket v \rrbracket\|_F \\ & \geq - \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} w_* \sqrt{\frac{p_{K_*} (p_{K_*} + d - 1) |F|}{d |K_*|}} \|a \nabla v\|_{K_*} \right) \|\llbracket v \rrbracket\|_F, \end{aligned} $$

using the trace inverse estimate and $\|\mathbf{n}\|_{L_\infty(F)} \le 1$. Setting, now,

$$ \tilde{\zeta}_* := \left( (d+1) \|\sqrt{a}\|_{L_\infty(K_*)} \sqrt{\frac{p_{K_*} (p_{K_*} + d - 1) |F|}{d |K_*|}} \right)^{-1}, $$
for $* \in \{+, -\}$, the above gives

$$ (2.7) \qquad \begin{aligned} & \left| \int_{\Gamma} \{ \Pi_{\mathbf{p}-1}(a\nabla v) \cdot \mathbf{n} \} \mathbf{w} \, [v] \, \mathrm{d}s \right| \\ & \le \sum_{F \subset \Gamma} \left( \sum_{* \in \{+, -\}} \frac{w_*}{(d+1)\zeta_*} \|\sqrt{a}\nabla v\|_{K_*} \right) \| [v] \|_F \\ & \le \frac{1}{3+\theta} \|\sqrt{a}\nabla v\|_{\Omega}^2 + (3+\theta) \left( \sum_{F \subset \Gamma} \max_{* \in \{+, -\}} \frac{w_*^2}{\zeta_*^2} \| [v] \|_F^2 \right), \end{aligned} $$

as before. Selecting, now, $\bar{\tau}|_F = (1+\theta)(3+\theta)\max_{* \in \{+, -\}} w_*^2 \bar{\zeta}_*^{-2}$, on every face $F \subset \Gamma$, we deduce

$$ B(v, v) \ge \frac{2}{3+\theta} \|\sqrt{a}\nabla_h v\|_{\Omega}^2 + \|\sqrt{c}v\|_{\Omega}^2 + \|\sqrt{\sigma - \bar{\tau}}[v]\|_{\Gamma}^2. $$

We select $\sigma$ satisfying $(1+\theta)\sigma = (3+\theta)\bar{\tau}$, thereby

$$ \sigma|_F := (3+\theta)^2 \max_{* \in \{+, -\}} \frac{w_*^2}{\bar{\zeta}_*^2}, $$

for each $F \subset \Gamma$, for $\theta \in (-1, 1]$, arbitrary $\sigma \ge 0$, for $\theta = -1$. With the above choice of $\sigma$, we deduce both the coercivity (2.3) and continuity (2.4) of the bilinear form, the latter for $\theta \in (-1, 1]$, while for $\theta = -1$, we have, instead (2.5) with the constant $C_{\mathrm{NIP}}$ replaced by $\bar{C}_{\mathrm{NIP}} := 2\|\max_{* \in \{+, -\}} w_*(\bar{\zeta}_*\sigma)^{-1}\|_{L_{\infty}(\Gamma)}$, for any $\sigma : \Gamma \to \mathbb{R}$ non-negative function.

**Remark 2.2.** The variant (2.6) can also be defined when the diffusion tensor vanishes locally upon adopting the convention that $\bar{\zeta}_*^{-1} = 0$ in regions where $a = 0$ also. For an alternative IPDG method that is robust with respect to diffusion degeneration we refer to [24] and its ‘robust’ extension in [19, Section 6].

### 3. A POSTERIORI ERROR BOUNDS FOR A CONFORMING RECONSTRUCTION

A key ingredient of a typical proof of a *posteriori* error bounds in the energy norm for non-conforming methods is a recovery of the approximate solution $u_h$ into the solution space. Indeed, employing the stability of local nodal averaging operators, Karakashian & Pascal proved the first rigorous *a posteriori* error bounds in the energy norm for the IPDG method discretising elliptic problems [31]. An alternative recovery/reconstruction employing the PDE problem has been recently used [12] for the proof of *a posteriori* error bounds in the energy norm for IPDG methods on polygonal/polyhedral meshes on simply connected domains. More recently, has been extended on hybrid high-order methods on simplicial meshes on the polytopic Lipschitz domain [17].

It is useful to define some contrast factors related to the diffusion coefficient $a$. For all $K \in \mathcal{T}$ and all $a \in \mathcal{V}_h$, we set

$$ (3.1) \qquad \mu_{\omega_k(K)}(a) := a_K^\sharp (a_{\omega_k(K)}^\flat)^{-1}, $$

with $a_K^\sharp := \sup_{\mathbf{x} \in K} |a(\mathbf{x})|$, $a_K^\flat := \inf_{\mathbf{x} \in K} |a(\mathbf{x})|$ and $a_{\omega_k(K)}^\flat := \inf_{\mathbf{x} \in \omega_k(K)} |a(\mathbf{x})|$, where $\omega_k(K)$ denotes the patch contains $k$-th layer elements. In addition,

$$ (3.2) \qquad \mu_a(a) := a_a^\sharp (a_a^\flat)^{-1}, $$

with $a_a^\sharp := \sup_{\mathbf{x} \in \omega_a} |a(\mathbf{x})|$ and $a_a^\flat := \inf_{\mathbf{x} \in \omega_a} |a(\mathbf{x})|$. Moreover, it is convenient to set

$$ (3.3a) \qquad a_F^\flat := \min(a_{K_+}^\flat, a_{K_-}^\flat), \quad \forall F \subset \partial K_+ \cap \partial K_- $$

$$ (3.3b) \qquad a_F^\flat := a_K^\flat, \quad \forall F \subset \partial K \subset \partial \Omega. $$
Following [17], we define the potential reconstruction $u_c \in H_g^1(\Omega)$ constructed from the dG solution $u_h \in S_T^p$, with $H_g^1(\Omega) := \{v \in H^1(\Omega) | v|_{\partial\Omega} = g\}$.

For the a posteriori error analysis, we begin by decomposing the error as

$$ (3.4) \qquad e := u - u_h = (u - u_c) + (u_c - u_h) =: e_c + e_d. $$

Observe the identity $B(e_d, v) = 0$ for any $v \in H^1(\Omega)$. Then, we have, respectively, (3.5)

$$ \begin{aligned} \frac{2}{3+\theta} \|e_c\|^2 &\le B(e_c, e_c) = B(e, e_c) - B(e_d, e_c) = \ell(e_c) - B(u_h, e_c) - B(e_d, e_c) \\ &= (\ell(e_c - v_h) - B(u_h, e_c - v_h)) - B(e_d, e_c) = T_1 + T_2, \end{aligned} $$

for any $v_h \in S_T^p$. Setting $\eta := e_c - \Pi_p e_c$ for brevity, we perform element-wise integration by parts to deduce

$$ \begin{aligned} T_1 &= \sum_{K \in \mathcal{T}} \int_K \left( (f + \nabla \cdot (a \nabla u_h) - cu_h) \eta \,dx - \int_{\partial K} a \nabla u_h \cdot \mathbf{n} \eta \,ds \right) \\ &\quad + \int_{\Gamma} \left( \{\langle (a \nabla u_h) \cdot \mathbf{n} \rangle_w \rangle[\eta] + \theta \{\langle a \Pi_{\mathbf{p}-1}(\nabla \eta) \cdot \mathbf{n} \rangle_w \langle u_h \rangle - \sigma \langle u_h \rangle \rangle[\eta] \right) \,ds \\ &\quad + \int_{\partial \Omega} g(\sigma\eta - \theta a \Pi_{\mathbf{p}-1}(\nabla \eta) \cdot \mathbf{n}) \,ds \\ (3.6) \qquad &= \sum_{K \in \mathcal{T}} \int_K \left( (f + \nabla \cdot (a \nabla u_h) - cu_h) \eta \,dx - \int_{\Gamma_{\text{int}}} \langle (a \nabla u_h) \cdot \mathbf{n} \rangle \{\langle \eta \rangle \}_{\mathbf{1}-w} \,ds \right) \\ &\quad + \int_{\Gamma} \left( \theta \{\langle a \Pi_{\mathbf{p}-1}(\nabla \eta) \cdot \mathbf{n} \rangle_w \langle u_h \rangle - \sigma \langle u_h \rangle \rangle[\eta] \right) \,ds \\ &\quad + \int_{\partial \Omega} g(\sigma\eta - \theta a \Pi_{\mathbf{p}-1}(\nabla \eta) \cdot \mathbf{n}) \,ds, \end{aligned} $$

since $\Pi_{\mathbf{p}-1}(\nabla u_h) = \nabla u_h$ and, upon employing the identity [33, eq. (3.15)]:

$$ \begin{aligned} & \sum_{K \in \mathcal{T}} \int_{\partial K} a \nabla u_h \cdot \mathbf{n} \eta \,ds \\ &= \int_{\Gamma} \{\langle (a \nabla u_h) \cdot \mathbf{n} \rangle_w \rangle[\eta] \,ds + \int_{\Gamma_{\text{int}}} \langle (a \nabla u_h) \cdot \mathbf{n} \rangle \{\langle \eta \rangle \}_{\mathbf{1}-w} \,ds, \end{aligned} $$

respectively, with $(\mathbf{1} - \mathbf{w})|_F := 1 - w_F$, for each $F \subset \Gamma_{\text{int}}$.

We are still to select suitable $v_h \in S_T^p$ in (3.13). Selecting $v_h := \Pi_p e_c$, we have

$$ (3.7) \qquad \|\eta\|_K \le C_{L_2}^{\text{app}} h_K p_K^{-1} \|\nabla e_c\|_K, $$

and

$$ (3.8) \qquad \|\eta\|_F \le C_{L_2}^{\text{app}} \sqrt{h_K p_K^{-1}} \|\nabla e_c\|_K, $$

for a suitable constant $C_{L_2}^{\text{app}} > 0$, independent of $h_K, p_K$ and of $e_c$, but may be dependent on the element shape. A computable upper bound to $C_{L_2}^{\text{app}}$ can be found in [16, 34]. Moreover, we have

$$ (3.9) \qquad \|\nabla \eta\|_K \le C_{L_2,K}^{\text{st}} p_K^{\frac{1}{2}} \|\nabla e_c\|_K. $$

with $C_{L_2,K}^{\text{st}} := C p_K^{\frac{1}{2}}$, see [15, Appendix C].

Defining, now, the residuals $\mathcal{R}_{\text{el}} : \Omega \setminus \Gamma \to \mathbb{R}$, $\mathcal{R}_{\text{fa}} : \Gamma_{\text{int}} \to \mathbb{R}$, and $\mathcal{R}_{\text{pe}} : \Gamma \to \mathbb{R}$, with

$$ \mathcal{R}_{\text{el}|K} := (a_K^{\flat})^{-\frac{1}{2}} h_K p_K^{-1} r_{\text{el},K}, $$
whereby $r_{\text{el},K} := (I - \Pi_{p_K})(f + \nabla \cdot (a\nabla u_h) - cu_h)|_K$, and

$$ \mathcal{R}_{\text{fa}}|_F := (a_F^b)^{-\frac{1}{2}} \max_{*\{+, -\}} (1 - w_*) \sqrt{2(d+1)h_{K_*} p_{K_*}^{-1}} \left\| (a\nabla u_h) \cdot \mathbf{n} \right\| $$

$$ \mathcal{R}_{\text{pe}}|_F := \begin{cases} \sqrt{2(d+1)\max_{*\{+, -\}} (a_{K_*}^b)^{-1}h_{K_*}p_{K_*}^{-1}\sigma\left\|u_h\right\|}, & \text{if } F \subset \Gamma_{\text{int}}; \\ \sqrt{2(d+1)(a_K^b)^{-1}h_Kp_K^{-1}\sigma(u_h-g)}, & \text{if } F \subset \partial K \cap \partial \Omega. \end{cases} $$

Returning to (3.6), we further estimate each term on the right-hand side: for the first term we have

$$ \sum_{K \in \mathcal{T}} \int_K r_{\text{el},K} \eta \, dx \le \sum_{K \in \mathcal{T}} C_{L_2}^{\text{app}} \|(h_K p_K^{-1}) r_{\text{el},K} \|_K \|\nabla e_c\|_K \\ \le \sum_{K \in \mathcal{T}} C_{L_2}^{\text{app}} \|\mathcal{R}_{\text{el}}\|_K \|\sqrt{a} \nabla e_c\|_K \le C_{L_2}^{\text{app}} \|\mathcal{R}_{\text{el}}\|_{\Omega} \|\sqrt{a} \nabla e_c\|_{\Omega}. $$

For the second term on the right-hand side of (3.6), we have

$$ \begin{aligned} & \int_{\Gamma_{\text{int}}} \left\| \left\| (a \nabla u_h) \cdot \mathbf{n} \right\| \right\|_{\{\eta\}} \mathbf{1}_{-w} \, ds \\ & \le C_{L_2}^{\text{app}} \sum_{F \subset \Gamma_{\text{int}}} \left\| \left\| (a \nabla u_h) \cdot \mathbf{n} \right\| \right\|_F \sum_{*\{+, -\}} (1 - w_*) \sqrt{h_{K_*} p_{K_*}^{-1}} \|\nabla e_c\|_{K_*} \\ & \le C_{L_2}^{\text{app}} \sum_{F \subset \Gamma_{\text{int}}} \max_{*\{+, -\}} (1 - w_*) \sqrt{(a_{K_*}^b)^{-1} h_{K_*} p_{K_*}^{-1}} \left\| \left\| (a \nabla u_h) \cdot \mathbf{n} \right\| \right\|_F \sum_{*\{+, -\}} \|\sqrt{a} \nabla e_c\|_{K_*} \\ & \le C_{L_2}^{\text{app}} \left( \sum_{F \subset \Gamma_{\text{int}}} \frac{1}{2(d+1)} \|\mathcal{R}_{\text{fa}}^{L_2}\|_F^2 \right)^{\frac{1}{2}} \left( \sum_{F \subset \Gamma_{\text{int}}} \left( \sum_{*\{+, -\}} \|\sqrt{a} \nabla e_c\|_{K_*} \right)^2 \right)^{\frac{1}{2}} \\ & \le C_{L_2}^{\text{app}} \left( \sum_{F \subset \Gamma_{\text{int}}} (d+1)^{-1} \|\mathcal{R}_{\text{fa}}^{L_2}\|_F^2 \right)^{\frac{1}{2}} \left( \sum_{F \subset \Gamma_{\text{int}}} \sum_{*\{+, -\}} \|\sqrt{a} \nabla e_c\|_{K_*}^2 \right)^{\frac{1}{2}} \\ & \le C_{L_2}^{\text{app}} \|\mathcal{R}_{\text{fa}}\|_{\Gamma_{\text{int}}} \|\sqrt{a} \nabla e_c\|_{\Omega}. \end{aligned} $$

To estimate the third term on right-hand side of (3.6), we set

$$ \begin{aligned} \xi_* &:= \zeta_* \|a^{-1/2}\|_{L_\infty(K_*)} \\ &= \left( (d+1) \|a|_{K_*} \mathbf{n}\|_{L_\infty(F)} \sqrt{\frac{p_{K_*} (p_{K_*} + d-1) |F|}{d|K_*|}} \right)^{-1}, \end{aligned} \quad (3.10) $$
for brevity, we split the interior skeleton $\Gamma_{\text{int}}$ and the boundary, we work as for
(2.7) and, additionally, we use (3.9), getting, respectively,

$$
\begin{align*}
& \left| \int_{\Gamma_{\text{int}}} \left\{ a \Pi_{\mathbf{p}-1} (\nabla \eta) \cdot \mathbf{n} \right\} \mathbf{w} \left[ \left[ u_h \right] \right] \mathrm{d}s \right| \\
& \le \sum_{F \subset \Gamma_{\text{int}}} \left( \sum_{* \in \{+, -\}} \frac{w_*}{(d+1)\xi_*} \| \nabla \eta \|_{K_*} \right) \left\| \left[ u_h \right] \right\|_F \\
& \le \sum_{F \subset \Gamma_{\text{int}}} \left( \sum_{* \in \{+, -\}} \frac{w_* C_{L_2, K^*}^{\text{st}}}{(d+1)\xi_*} \| \nabla e_c \|_{K_*} \right) \left\| \left[ u_h \right] \right\|_F \\
& \le \sum_{F \subset \Gamma_{\text{int}}} \left( \sum_{* \in \{+, -\}} \frac{w_* C_{L_2, K^*}^{\text{st}}}{(d+1)\xi_*(a_{K_*}^b)^{\frac{1}{2}}} \| \sqrt{a} \nabla e_c \|_{K_*} \right) \left\| \left[ u_h \right] \right\|_F \\
& \le \sum_{F \subset \Gamma_{\text{int}}} \left( \sum_{* \in \{+, -\}} (d+1)^{-1} \| \sqrt{a} \nabla e_c \|_{K_*} \right) \\
& \qquad \times \max_{* \in \{+, -\}} \frac{w_*}{\xi_*} C_{L_2, K^*}^{\text{st}} (a_{K_*}^b)^{-\frac{1}{2}} \left\| \left[ u_h \right] \right\|_F \\
& \le \left( \sum_{F \subset \Gamma_{\text{int}}} \sum_{* \in \{+, -\}} \frac{2}{(d+1)^2} \| \sqrt{a} \nabla e_c \|_{K_*}^2 \right)^{\frac{1}{2}} \\
& \qquad \times \left( \sum_{F \subset \Gamma_{\text{int}}} \max_{* \in \{+, -\}} \left( \frac{w_*}{\xi_*} \frac{C_{L_2, K^*}^{\text{st}}}{(a_{K_*}^b)^{\frac{1}{2}}} \right)^2 \left\| \left[ u_h \right] \right\|_F^2 \right)^{\frac{1}{2}} \\
& \le \sqrt{\frac{2}{d+1}} \| \sqrt{a} \nabla e_c \|_{\Omega} \left( \sum_{F \subset \Gamma_{\text{int}}} \max_{* \in \{+, -\}} \left( \frac{w_*}{\xi_*} \frac{C_{L_2, K^*}^{\text{st}}}{(a_{K_*}^b)^{\frac{1}{2}}} \right)^2 \left\| \left[ u_h \right] \right\|_F^2 \right)^{\frac{1}{2}} \\
& = \| \chi \left[ u_h \right] \|_{\Gamma_{\text{int}}} \| \sqrt{a} \nabla e_c \|_{\Omega},
\end{align*}
$$

for $\theta \in [-1, 1]$, with $\chi : \Gamma \to \mathbb{R}$, defined face-wise by

$$
\chi|_F := \sqrt{2}(d+1)^{-\frac{1}{2}} \max_{* \in \{+, -\}} \frac{w_* C_{L_2, K^*}^{\text{st}}}{\xi_* (a_{K_*}^b)^{\frac{1}{2}}},
$$

for $F \subset \Gamma_{\text{int}}$ and $\chi|_F := \sqrt{2}(d+1)^{-\frac{1}{2}}\xi_+^{-1}C_{L_2,K^+}^{\text{st}}(a_{K_+}^b)^{-\frac{1}{2}}$ for $F \subset \partial\Omega$ (the latter to be required below). Also, we have

$$
\begin{align*}
& \int_{\Gamma_{\text{int}}} \sigma \llbracket u_h \rrbracket \llbracket \eta \rrbracket \mathrm{d}s \\
&\le \sum_{F \subset \Gamma_{\text{int}}} \| \sigma \llbracket u_h \rrbracket \|_F \| \llbracket \eta \rrbracket \|_F \\
&\le C_{L_2}^{\text{app}} \sum_{F \subset \Gamma_{\text{int}}} \| \sigma \llbracket u_h \rrbracket \|_F \left( \sum_{* \in \{+, -\}} \sqrt{h_{K_*} p_{K_*}^{-1}} \| \nabla e_c \|_{K_*} \right) \\
&\le C_{L_2}^{\text{app}} \sum_{F \subset \Gamma_{\text{int}}} \max_{* \in \{+, -\}} \sqrt{(a_{K_*}^b)^{-1} h_{K_*} p_{K_*}^{-1}} \| \sigma \llbracket u_h \rrbracket \|_F \left( \sum_{* \in \{+, -\}} \| \sqrt{a} \nabla e_c \|_{K_*} \right) \\
&\le \sqrt{2} C_{L_2}^{\text{app}} \left( \sum_{F \subset \Gamma_{\text{int}}} \max_{* \in \{+, -\}} \frac{h_{K_*}}{a_{K_*}^b p_{K_*}} \| \sigma \llbracket u_h \rrbracket \|_F^2 \right)^{\frac{1}{2}} \left( \sum_{F \subset \Gamma_{\text{int}}} \sum_{* \in \{+, -\}} \| \sqrt{a} \nabla e_c \|_{K_*}^2 \right)^{\frac{1}{2}} \\
&\le C_{L_2}^{\text{app}} \left( 2(d+1) \sum_{F \subset \Gamma_{\text{int}}} \max_{* \in \{+, -\}} \frac{h_{K_*}}{a_{K_*}^b p_{K_*}} \| \sigma \llbracket u_h \rrbracket \|_F^2 \right)^{\frac{1}{2}} \| \sqrt{a} \nabla e_c \|_{\Omega} \\
&= C_{L_2}^{\text{app}} \| \mathcal{R}_{\text{pe}} \|_{\Gamma_{\text{int}}} \| \sqrt{a} \nabla e_c \|_{\Omega}.
\end{align*}
$$
For the respective boundary terms, working completely analogously to the respective interior terms, we deduce

$$ \int_{\partial\Omega} (u_h - g)(\theta a \Pi_{\mathbf{p}-1}(\nabla\eta) \cdot \mathbf{n} u_h - \sigma u_h \eta) \, ds \\ \leq (\|\chi(u_h - g)\|_{\partial\Omega} + C_{L_2}^{\text{app}} \|\mathcal{R}_{\text{pe}}^{L_2}\|_{\partial\Omega}) \|\sqrt{a} \nabla e_c\|_{\Omega}. $$

Combining the above, and performing straightforward further estimation, we have the following result.

**Theorem 3.1.** Assume $a \in [W^{1,\infty}(\Omega, \mathcal{T})]^{d \times d}$, let $u \in H^1(\Omega)$ be the solution to (1.1) and let $u_h \in S_{\mathcal{T}}^{\mathbf{p}}$ its approximation satisfying (2.1). Then, we have the error bound

$$ (3.12) \qquad \|e_c\| \le C_1 (\|\mathcal{R}_{\text{el}}\|_{\Omega} + \|\mathcal{R}_{\text{fa}}\|_{\Gamma_{\text{int}}} + \|\mathcal{R}_{\text{pe}}\|_{\Gamma_{\text{int}}} + \|\mathcal{R}_{\text{pe}}\|_{\partial\Omega}) \\ \qquad \qquad + C_2 (\|\chi\|_{u_h}\|_{\Gamma_{\text{int}}} + \|\chi(u_h - g)\|_{\partial\Omega}) + C_3 \|e_d\|, $$

whereby $C_2 := \frac{3+\theta}{2}$, $C_1 := C_{L_2}^{\text{app}} C_2$ and $C_3 := \frac{5+\theta}{2}$, for $\theta \in [-1, 1]$.

*Proof.* Combining the above estimates, we get

$$ \frac{2}{3+\theta} \|e_c\|^2 \le C_{L_2}^{\text{app}} \left( \|\mathcal{R}_{\text{el}}\|_{\Omega} + \|\mathcal{R}_{\text{fa}}\|_{\Gamma_{\text{int}}} + \|\mathcal{R}_{\text{pe}}\|_{\Gamma_{\text{int}}} + \|\mathcal{R}_{\text{pe}}\|_{\partial\Omega} \right) \|\sqrt{a} \nabla e_c\|_{\Omega} \\ + \left( \|\chi\|_{u_h}\|_{\Gamma_{\text{int}}} + \|\chi(u_h - g)\|_{\partial\Omega} \right) \|\sqrt{a} \nabla e_c\|_{\Omega} + \frac{5+\theta}{3+\theta} \|e_d\| \|e_c\|. $$

for $\theta \in [-1, 1]$. Since $\|\sqrt{a} \nabla e_c\|_{\Omega} \le \|e_c\|$, elementary calculations give the result.

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

We note that (3.12) is not an a posteriori error bound, as it contains $e_d$ on the right-hand side. We will arrive to a fully computable bound will be completed below by showing an a posteriori error bound for $\|e_d\|$.

Next, we also show a bound for the variant of the classical weighted IPDG method (2.6), which is valid for $a \in [L_{\infty}(\Omega)]^{d \times d}$ also. In this case, noting that $\Pi_{\mathbf{p}-1}(a\nabla u_h)|_K \in L_2(K)$, for all $K \in \mathcal{T}$, we add and subtract the term

$$ \int_{\Omega} \Pi_{\mathbf{p}-1}(a \nabla_h u_h) \cdot \nabla_h \eta \, dx $$

on the right-hand side of the error equation, to deduce

$$ (3.13) \qquad \begin{aligned} l(\eta) - \mathcal{B}(u_h, \eta) = & \sum_{K \in \mathcal{T}} \int_K (f + \nabla \cdot \Pi_{\mathbf{p}-1}(a \nabla u_h) - cu_h) \hat{\eta} \, dx \\ & - \int_{\Gamma_{\text{int}}} \left[ \Pi_{\mathbf{p}-1}(a \nabla u_h) \cdot \mathbf{n} \right] \{\eta\}_{1-\mathbf{w}} \, ds \\ & - \int_{\Gamma_{\text{int}}} \left( \theta \left\{ \Pi_{\mathbf{p}-1}(a \nabla \eta) \cdot \mathbf{n} \right\}_{\mathbf{w}} \left[ u_h \right] - \sigma \left[ u_h \right] \left[ \eta \right] \right) \, ds \\ & + \int_{\partial\Omega} (g - u_h)(\sigma\eta - \theta\Pi_{\mathbf{p}-1}(a \nabla \eta) \cdot \mathbf{n}) \, ds \\ & + \int_{\Omega} (I - \Pi_{\mathbf{p}-1})(a \nabla_h u_h) \cdot \nabla_h \eta \, dx, \end{aligned} $$

for $\eta = e_c - v_h$, $v_h \in S_{\mathcal{T}}^{\mathbf{p}}$.

We select $v_h := \Pi_{\mathbf{p}}e_c$ and define the residuals $\mathfrak{R}_{\text{el}} : \Omega \setminus \Gamma \to \mathbb{R}$, and $\mathfrak{R}_{\text{fa}} : \Gamma_{\text{int}} \to \mathbb{R}$, with

$$ \mathfrak{R}_{\text{el}}|_K := (a_K^b)^{-\frac{1}{2}} h_K p_K^{-1} (I - \Pi_{p_K}) (f - cu_h), $$

$$ \mathcal{R}_{\text{fa}}|_F := (a_F^b)^{-\frac{1}{2}} \max_{*\{+, -\}} (1 - w_*) \sqrt{s(d+1)} h_K p_K^{-1} \left[ \Pi_{\mathbf{p}-1}(a \nabla u_h) \cdot \mathbf{n} \right], $$
as well as the non-standard “oscillation” term $\mathfrak{R}_a : \Omega \setminus \Gamma \to \mathbb{R}$, given by

$$ \mathfrak{R}_a|_K := (a_K^b)^{-\frac{1}{2}} C_{L_2,K}^{\text{st}} (I - \Pi_{\mathbf{p}-1})(a \nabla u_h)|_K $$

for all $K \in \mathcal{T}$.

Returning to (3.13), we estimate each term on the right-hand side: for the first term we have

$$ \sum_{K \in \mathcal{T}} \int_K (f + \nabla \cdot \Pi_{\mathbf{p}-1}(a \nabla u_h) - cu_h) \eta \, dx = \int_{\Omega} (f - cu_h) \eta \, dx \\ \leq C_{L_2}^{\text{app}} \|\mathfrak{R}_{\text{el}}\|_{\Omega} \|\sqrt{a} \nabla e_c\|_{\Omega}. $$

For the third term, we let

$$ (3.14) \qquad \bar{\xi}_* := \left( (d+1)\|a\|_{L_{\infty}(K_*)} \sqrt{\frac{p_{K_*}(p_{K_*}+d-1)|F|}{d|K_*|}} \right)^{-1} $$

and, following similar steps as in (3.11), we deduce

$$ \left| \int_{\Gamma_{\text{int}}} \left\{ \Pi_{\mathbf{p}-1}(a \nabla \eta) \cdot \mathbf{n} \right\} \mathbf{w} \left[ u_h \right] \, ds \right| \leq \| \bar{\chi} \left[ u_h \right] \|_{\Gamma_{\text{int}}} \| \sqrt{a} \nabla e_c \|_{\Omega}, $$

for $\bar{\chi} : \Gamma \to \mathbb{R}$, with

$$ \bar{\chi}|_F := \sqrt{2(d+1)^{-\frac{1}{2}}} \max_{* \in \{+, -\}} (w_* \bar{\xi}_*^{-1}) C_{L_2, K_*}^{\text{st}} (a_{K_*}^b)^{-\frac{1}{2}}, $$

for every face $F \subset \Gamma_{\text{int}}$ and $\bar{\chi}|_F := \sqrt{2(d+1)^{-\frac{1}{2}}} \bar{\xi}_* + C_{L_2, K_*}^{\text{st}} (a_{K_*}^b)^{-\frac{1}{2}}$ for $F \subset \partial\Omega$.

The bound for the discontinuity-penalisation component from (3.11) is still valid.

For the boundary terms, we have

$$ \int_{\partial\Omega} (u_h - g)(\theta\Pi_{\mathbf{p}-1}(a\nabla\eta) \cdot \mathbf{n}u_h - \sigma u_h\eta) \, ds \\ \leq (\|\bar{\chi}(u_h - g)\|_{\partial\Omega} + C_{L_2}^{\text{app}} \|\mathcal{R}_{\text{pe}}\|_{\partial\Omega}) \|\sqrt{a}\nabla e_c\|_{\Omega}. $$

Finally, for the last term on the right-hand side of (3.13), we have

$$ \int_{\Omega} (I - \Pi_{\mathbf{p}-1})(a \nabla_h u_h) \cdot \nabla_h \eta \, dx \leq \|\mathfrak{R}_a\|_{\Omega} \|\sqrt{a} \nabla e_c\|_{\Omega}. $$

Combining the above, along with straightforward arguments, we arrive at the following result.

**Theorem 3.2.** For $a \in [L_{\infty}(\Omega)]^{d \times d}$, let $u \in H^1(\Omega)$ be the solution to (1.1) and let $u_h \in S_{\mathcal{T}}^p$ its approximation satisfying (2.6). Then, we have the bound

$$ \|\|e_c\|\| \leq C_1 \left( \|\mathfrak{R}_{\text{el}}\|_{\Omega} + \|\mathfrak{R}_{\text{fa}}\|_{\Gamma_{\text{int}}} + \|\mathcal{R}_{\text{pe}}\|_{\Gamma_{\text{int}}} + \|\mathcal{R}_{\text{pe}}\|_{\partial\Omega} \right) \\ + C_2 \left( \|\bar{\chi}\left[ u_h \right]\|_{\Gamma_{\text{int}}} + \|\bar{\chi}(u_h - g)\|_{\partial\Omega} + \|\mathfrak{R}_a\|_{\Omega} \right) + C_3 \|e_d\|, \quad (3.15) $$

whereby $C_2 := \frac{3+\theta}{2} \|a^{-\frac{1}{2}}\|_{L_{\infty}(\Omega)}$ and $C_1 := C_{L_2}^{\text{app}} C_2$, and $C_3 := \frac{5+\theta}{2}$, for $\theta \in [-1, 1]$.

### 4. COMPLETING THE A POSTERIORI ERROR BOUNDS

To arrive at a completely computable upper bound, we need to further estimate the non-conforming part of the error $\|\|e_d\|\|$.

We start by defining the function $u_c \in H_g^1(\Omega)$ introduced in (3.4). We solve local minimization problems in $H^1(\omega_a)$ with suitable boundary conditions for every mesh vertex $a \in V_h$ (recall that $\omega_a$ denotes the star associated with the vertex $a$). Let $\psi_a$ be the hat basis function equal to 1 at $a$ and having a support in the vertex
star $\omega_a$. Recall that the hat basis functions satisfy the following partition-of-unity property:

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

Before introducing the potential reconstruction, we introduce an auxiliary SIPDG bilinear form $A: H^1(\Omega, \mathcal{T}) \times H^1(\Omega, \mathcal{T}) \to \mathbb{R}$, given by

$$ (4.2) \qquad A(w, v) := \int_{\Omega} a \nabla_h w \cdot \nabla_h v \, dx + \int_{\Gamma} \sigma[\|w\|\|v\|] \, ds \\ - \int_{\Gamma} (\{a \Pi_{\mathbf{p}-1}(\nabla w) \cdot \mathbf{n}\}_{\mathbf{w}}[\|v\|] + \{a \Pi_{\mathbf{p}-1}(\nabla v) \cdot \mathbf{n}\}_{\mathbf{w}}[\|w\|]) \, ds, $$

for $w, v \in H^1(\Omega, \mathcal{T})$. Following the analysis in Section 2, we have

$$ (4.3) \qquad A(w, w) \ge \frac{1}{2} \|w\|_d^2 \quad \text{and} \quad A(v, w) \le 2\|v\|_d \|w\|_d, $$

$$ \text{with } \|v\|_d := (\|\sqrt{a} \nabla_h w\|_{\Omega}^2 + \|\sqrt{\sigma}[w]\|_{\Gamma}^2)^{\frac{1}{2}}. $$

**Definition 4.1** (Patchwise and global potential reconstruction). For all $a \in \mathcal{V}_h$, let $u_c^a \in H_g^1(\omega_a)$ be the solution of the following well-posed problem:

$$ (4.4) \qquad A(u_c^a, v_a) = A(\psi_a u_h, v_a), \quad \forall v_a \in H_0^1(\omega_a), $$

with

$$ (4.5) \qquad H_g^1(\omega_a) := \{v \in H^1(\omega_a) \mid v|_{\partial\omega_a \cap \partial\Omega} = \psi_a g \text{ and } v|_{\partial\omega_a \cap \partial\Omega} = 0\}. $$

An equivalent definition is

$$ (4.6) \qquad u_c^a := \arg \min_{\rho_a \in H_g^1(\omega_a)} A(\rho_a - (\psi_a u_h), \rho_a - (\psi_a u_h)). $$

Then, extending $u_c^a$ by zero to $\Omega$, we set

$$ (4.7) \qquad u_c := \sum_{a \in \mathcal{V}_h} v_c^a. $$

Notice that we indeed have $u_c \in H_g^1(\Omega)$ as required in (3.4); this follows from the partition-of-unity property (4.1) and the definition (4.5) of $H_g^1(\omega_a)$.

**Lemma 4.2** (Local nonconformity error). Let $u_c^a$ solve (4.4) and set $e_d^a := u_c^a - \psi_a u_h$ for all $a \in \mathcal{V}_h$. The following holds:

$$ (4.8) \qquad \|\sqrt{a} \nabla_h e_d^a\|_{\omega_a} \le C_d^a \left( \|\tilde{\chi}[u_h]\|_{\mathcal{F}_a \cap \Gamma_{\text{int}}}^2 + \|\tilde{\chi}(u_h - g)\|_{\mathcal{F}_a \cap \partial\Omega}^2 + \|\mathfrak{R}_{\text{ta}}^{\text{QI}}\|_{\mathcal{F}_a}^2 \right)^{\frac{1}{2}}, $$

with

$$ (4.9) \qquad \mathcal{R}_{\text{ta}}^{\text{QI}}|_F := \begin{cases} \sqrt{(a_F^b)^{-1} \max_{\{+, -\}} h_{K_*} p_{K_*}^{-1}} [\nabla u_h \times \mathbf{n}], & \text{if } F \subset \mathcal{F}_a \cap \Gamma_{\text{int}}; \\ \sqrt{(a_F^b)^{-1} h_{K_+} p_{K_+}^{-1}} (\nabla (u_h - g_D) \times \mathbf{n}_\Omega), & \text{if } F \subset \mathcal{F}_a \cap \partial\Omega, \end{cases} $$

where the constant $C_d^a$ depends on the mesh shape-regularity, on $\gamma$ and on $\mu_a(A)^{\frac{1}{2}}$, but is independent of $h$ and $p$.

*Proof.* The proof primarily follows the approach detailed in [12, Section 3.1], applied to each local vertex patch. The key ingredients include the Helmholtz decomposition for simply connected domains [26, Theorem 3] and the hp-optimal quasi-interpolation operator introduced in [32]. The proof can be derived following techniques in the proof of [17, Lemma 5.4] for the $H^1$-reconstruction. $\square$

Next, we will derive the $L^2$-norm error bound for $e_d^a$.
**Lemma 4.3** (Local nonconformity $L_2$-norm error). Let $u_c^a$ solve (4.4) and set $e_d^a := u_c^a - \psi_a u_h$ for all $a \in V_h$. The following holds:

$$ (4.10) \qquad \|\sqrt{c}e_d^a\|_{\omega_a} \le \tilde{C}_d^a h_{\omega_a} \left( \|\tilde{\chi}_*[|u_h]\|_{\mathcal{F}_a \cap \Gamma_{\text{int}}}^2 + \|\tilde{\chi}_*(u_h - g)\|_{\mathcal{F}_a \cap \partial\Omega}^2 + \|\mathfrak{R}_{\text{ta},*}^{\text{QI}}\|_{\mathcal{F}_a}^2 \right)^{\frac{1}{2}}, $$

with

$$ (4.11) \qquad \tilde{\chi}_*|_F := \begin{cases} \sqrt{2c_a^\sharp}(d+1)^{-\frac{1}{2}} C_{\text{QI}}^{\text{st}} \max_{* \in \{+, -\}} \frac{w_*}{\xi_*}, & \text{if } F \subset \mathcal{F}_a \cap \Gamma_{\text{int}}; \\ \sqrt{2c_a^\sharp}(d+1)^{-\frac{1}{2}} C_{\text{QI}}^{\text{st}} \frac{w_+}{\xi_+}, & \text{if } F \subset \mathcal{F}_a \cap \partial\Omega, \end{cases} $$

and

$$ (4.12) \qquad \mathcal{R}_{\text{ta},*}^{\text{QI}}|_F := \begin{cases} \sqrt{c_a^\sharp} \max_{* \in \{+, -\}} h_{K_*} p_{K_*}^{-1} [[\nabla u_h \times \mathbf{n}]], & \text{if } F \subset \mathcal{F}_a \cap \Gamma_{\text{int}}; \\ \sqrt{c_a^\sharp} h_{K_+} p_{K_+}^{-1} (\nabla (u_h - g_D) \times \mathbf{n}_\Omega) & \text{if } F \subset \mathcal{F}_a \cap \partial\Omega, \end{cases} $$

where the constant $\tilde{C}_d^a$ depends on the mesh shape-regularity and on $\gamma$, but is independent of $h$ and $p$.

*Proof.* First, recalling the definition of $u_c^a$, we have $\|e_d^a\|_{\mathcal{F}_a \cap \Gamma_{\text{int}}} = -[\psi_a u_h]_{\mathcal{F}_a \cap \Gamma_{\text{int}}}$ and $e_d^a|_{\partial\omega_a \setminus \mathcal{F}_a} = 0$. Next, applying the broken Poincaré inequality from [7], we have

$$ \|e_d^a\|_{\omega_a}^2 \le C_{\text{bPF}}^2 h_{\omega_a}^2 \left( \|\nabla_h e_d^a\|_{\omega_a}^2 + \|h^{-\frac{1}{2}} \Pi_0^F [\psi_a u_h]\|_{\mathcal{F}_a \cap \Gamma_{\text{int}}}^2 \right. \\ \left. + \|h^{-\frac{1}{2}} \Pi_0^F (\psi_a(u_h - g))\|_{\mathcal{F}_a \cap \partial\Omega}^2 \right), $$

with $C_{\text{bPF}} > 0$ only depends on the mesh shape regularity. Invoking the stability of the $L^2$-orthogonal projection, $\|\psi_a\|_{L_\infty(\omega_a)} = 1$, mesh shape regularity and the result (4.8) by choosing diffusion tensor $a$ to be identity matrix, we infer

$$ \begin{aligned} \|\sqrt{c}e_d^a\|_{\omega_a}^2 &\le c_a^\sharp (C_{\text{bPF}} h_{\omega_a})^2 \left( \|\nabla_h e_d^a\|_{\omega_a}^2 + \|h^{-\frac{1}{2}}[u_h]\|_{\mathcal{F}_a \cap \Gamma_{\text{int}}}^2 + \|h^{-\frac{1}{2}}(u_h - g)\|_{\mathcal{F}_a \cap \partial\Omega}^2 \right) \\ &\le (\tilde{C}_d^a h_{\omega_a})^2 \left( \|\tilde{\chi}_*[|u_h]\|_{\mathcal{F}_a \cap \Gamma_{\text{int}}}^2 + \|\tilde{\chi}_*(u_h - g)\|_{\mathcal{F}_a \cap \partial\Omega}^2 + \|\mathfrak{R}_{\text{ta},*}^{\text{QI}}\|_{\mathcal{F}_a}^2 \right). \end{aligned} \quad \square $$

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

*Remark 4.4* (L$_2$-norm bound). The above $L_2$-norm error bound of $e_d^a$ is $h$-optimal but $p$-subopitmal by one order in the DG norm. Moreover, as we do not assume reaction coefficient $c$ to be strictly positive, the above bound does not depend on the contrast of $c$.

**Theorem 4.5** (Global bound on nonconformity error). The following holds:

$$ \begin{aligned} & \|\sqrt{a}\nabla_h e_d\|_{\Omega}^2 + \|\sqrt{c}e_d\|_{\Omega}^2 \\ & \le C_d^2 \sum_{K \in \mathcal{T}} \left( \|\tilde{\chi}[|u_h]\|_{\partial K \cap \Gamma_{\text{int}}}^2 + \|\tilde{\chi}(u_h - g)\|_{\partial K \cap \partial\Omega}^2 \right. \\ & \qquad \left. + \|\mathfrak{R}_{\text{ta}}^{\text{QI}}\|_{\partial K}^2 + h_K^2 \|\tilde{\chi}_*[|u_h]\|_{\partial K \cap \Gamma_{\text{int}}}^2 \right. \\ & \qquad \left. + h_K^2 \|\tilde{\chi}_*(u_h - g)\|_{\partial K \cap \partial\Omega}^2 + h_K^2 \|\mathfrak{R}_{\text{ta},*}^{\text{QI}}\|_{\partial K}^2 \right), \end{aligned} \quad (4.13) $$

where the constant $C_d$ depends on the mesh shape-regularity, on $\gamma$, on $\max_{a \in V_h} \chi_a(A)^{\frac{1}{2}}$ and $\max_{a \in V_h} (c_a^\sharp)^{\frac{1}{2}}$, but is independent of $h$ and $k$.
*Proof.* Using that $u_c|_K = \sum_{a \in \mathcal{V}_K} u_c^a|_K$ for all $K \in \mathcal{T}$, where $\mathcal{V}_K$ is the set of vertices of $K$, the partition-of-unity property (4.1), and the triangle inequality, we infer that

$$
\begin{align*}
\|\sqrt{a}\nabla_h e_d\|_{\Omega}^2 &= \sum_{K \in \mathcal{T}} \|\sqrt{a}\nabla(u_c - u_h|_K)\|_K^2 \\
&= \sum_{K \in \mathcal{T}} \left\| \sum_{a \in \mathcal{V}_K} \sqrt{a}\nabla(u_c^a - \psi_a u_h|_K) \right\|_K^2 \\
&\le \sum_{K \in \mathcal{T}} \sum_{a \in \mathcal{V}_K} (d+1) \|\sqrt{a}\nabla(u_c^a - \psi_a u_h|_K)\|_K^2 \\
&= \sum_{a \in \mathcal{V}_h} (d+1) \|\sqrt{a}\nabla_h e_d^a\|_{\omega_a}^2.
\end{align*}
$$

Similar, we infer

$$
\begin{align*}
\|\sqrt{c}e_d\|_{\Omega}^2 &= \sum_{K \in \mathcal{T}} \|\sqrt{c}(u_c - u_h|_K)\|_K^2 \\
&= \sum_{K \in \mathcal{T}} \left\| \sum_{a \in \mathcal{V}_K} \sqrt{c}(u_c^a - \psi_a u_h|_K) \right\|_K^2 \\
&\le \sum_{K \in \mathcal{T}} \sum_{a \in \mathcal{V}_K} (d+1) \|\sqrt{c}(u_c^a - \psi_a u_h|_K)\|_K^2 \\
&= \sum_{a \in \mathcal{V}_h} (d+1) \|\sqrt{c}e_d^a\|_{\omega_a}^2.
\end{align*}
$$

Invoking Lemma 4.2, Lemma 4.3 and using that $h_{\omega_a} \approx h_K$ for all $K \in \omega_a$ and all $a \in \mathcal{V}_h$ completes the proof. □

Upon completing the upper bound, some comments on lower bounds are also in order. The individual a posteriori error estimators of the previous theorem are familiar from known bounds in the literature [29, 20, 12], albeit with different constants, we prefer not to include these in the discussion here in the interest of brevity, as there are no essential improvements to be provided in this context, compared to the results in the literature.

## 5. NUMERICAL EXAMPLES

In this section, we present numerical examples to illustrate our theoretical results.

**5.1. Example 1: Adaptive algorithm for singular solution.** We select $f, a := I_{2 \times 2}$ and Dirichlet boundary conditions on the L-shaped domain $\Omega := (-1, 1)^2 \setminus \{(0, 1) \times (-1, 0)\}$, so that the exact solution is in polar coordinates

$$ (5.1) \qquad u = r^{\frac{2}{3}} \sin(2\theta/3) + \exp(-1000((x-0.5)^2 + (y-0.25)^2)) \\ \qquad \qquad \qquad + \exp(-1000((x-0.5)^2 + (y-0.75)^2)) $$

We test an adaptive algorithm driven by the a posteriori error estimator from Section 4. The adaptive algorithm starts from a coarse mesh and uses the estimator to mark mesh cells for refinement through a bulk-chasing criterion (also known as Dörfler's marking). The adaptive algorithm can be classically described as

SOLVE → ESTIMATE → MARK → REFINE.

We first test the convergence rate of the above adaptive algorithm with $p \in \{1, 2, 3, 4\}$ and a bulk-chasing criterion set to 30%. The energy error and the a posteriori error estimator are reported in Figures 1 and 2 with convergence rates computed as a
<figure><img src="image_3.png" /><figcaption>FIGURE 1. Example 1. DG-norm error and a posteriori error estimator as a function of DoFs for $p \in \{1, 2\}$.</figcaption></figure>

function of numerical degrees of freedom (DoFs). We observe that the energy error and the a posteriori error estimator converge at the optimal rate $O(DoFs^{-p/2})$. Moreover, we observe in Figure 3 that the effectivity index remains well behaved as a function of DoFs and that it slightly increases with the polynomial degree $p \ge 1$, with values between 1.5 and 2.8 asymptotically. In addition, we notice that there
<figure><img src="image_4.png" /><figcaption>FIGURE 2. Example 1. DG-norm error and a posteriori error estimator as a function of DoFs for $p \in \{3, 4\}$.</figcaption></figure>

is a pre-asymptotic regimes that the effectivity index grow to 3.5. We point out that although the solution exhibits a corner singularity, it appears that the two smooth functions require more degrees of freedom to resolve than the singularity. This behaviour can be observed in the sequence of meshes produced by the adaptive algorithm with $p = 4$, shown in Figure 4.
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Example 1. Effectivity index as a function of DoFs for $p \in \{1, 2, 3, 4\}$.</figcaption></figure>

REFERENCES

[1] M. AINSWORTH, *A synthesis of a posteriori error estimation techniques for conforming, non-conforming and discontinuous Galerkin finite element methods*, in Recent advances in adaptive computation, vol. 383 of Contemp. Math., Amer. Math. Soc., Providence, RI, 2005, pp. 1-14.

[2] M. AINSWORTH AND R. RANKIN, *Fully computable bounds for the error in nonconforming finite element approximations of arbitrary order on triangular elements*, SIAM J. Numer. Anal., 46 (2008), pp. 3207-3232.

[3] C. ANNAVARAPU, M. HAUTEFEUILLE, AND J. E. DOLBOW, *A robust Nitsche's formulation for interface problems*, Comput. Methods Appl. Mech. Engrg., 225/228 (2012), pp. 44-54.

[4] D. ARNOLD, *An interior penalty finite element method with discontinuous elements*, SIAM J. Numer. Anal., 19 (1982), pp. 742-760.

[5] G. BAKER, *Finite element methods for elliptic equations using nonconforming elements*, Math. Comp., 31 (1977), pp. 45-59.

[6] 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.

[7] S. C. BRENNER, *Poincaré-Friedrichs inequalities for piecewise H¹ functions*, SIAM J. Numer. Anal., 41 (2003), pp. 306-324.

[8] E. BURMAN AND P. ZUNINO, *A domain decomposition method based on weighted interior penalties for advection-diffusion-reaction problems*, SIAM J. Numer. Anal., 44 (2006), pp. 1612-1638.

[9] A. CANGIANI, Z. DONG, E. GEORGOULIS, AND P. HOUSTON, *hp-Version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes*, ESAIM: M2AN, 50 (2016), pp. 699-725.

[10] 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.

[11] ——, *hp-version discontinuous Galerkin methods on essentially arbitrarily-shaped elements*, Math. Comp., 91 (2022), pp. 1-35.

[12] ——, *A posteriori error estimates for discontinuous Galerkin methods on polygonal and polyhedral meshes*, SIAM J. Numer. Anal., 61 (2023), pp. 2352-2380.

[13] A. CANGIANI, Z. DONG, E. H. GEORGOULIS, AND P. HOUSTON, *hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes*, SpringerBriefs in Mathematics, Springer, Cham, 2017.
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Example 1. The sequence of meshes generated by the adaptive dG algorithm with $p = 4$.</figcaption></figure>

[14] A. CANGIANI, E. GEORGOULIS, AND P. HOUSTON, hp-Version discontinuous Galerkin methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci., 24 (2014), pp. 2009–2041.

[15] C. CARSTENSEN, B. GRÄSS LE, AND N. T. TRAN, *Adaptive hybrid high-order method for guaranteed lower eigenvalue bounds*, Numer. Math., 156 (2024), pp. 813–851.

[16] A. CHERNOV, *Optimal convergence estimates for the trace of the polynomial $L^2$-projection operator on a simplex*, Math. Comp., 81 (2012), pp. 765–787.

[17] Z. DONG AND A. ERN, hp-error analysis of mixed-order hybrid high-order methods for elliptic problems on simplicial meshes. working paper or preprint, Oct. 2024.

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

[19] ——, *Robust interior penalty discontinuous Galerkin methods*, J. Sci. Comput., 92 (2022), pp. Paper No. 57, 23.

[20] Z. DONG, L. MASCOTTO, AND O. J. SUTTON, *Residual-based a posteriori error estimates for hp-discontinuous Galerkin discretizations of the biharmonic problem*, SIAM J. Numer. Anal., 59 (2021), pp. 1273–1298.

[21] M. DRYJA, *On discontinuous Galerkin methods for elliptic problems with discontinuous coefficients*, vol. 3, 2003, pp. 76–85. Dedicated to Raytcho Lazarov.

[22] A. ERN AND A. STEPHANSEN, *A posteriori energy-norm error estimates for advection-diffusion equations approximated by weighted interior penalty methods*, J. Comp. Math., 26 (2008), pp. 488–510.

[23] A. ERN, A. F. STEPHANSEN, AND P. ZUNINO, *A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity*, IMA J. Numer. Anal., 29 (2009), pp. 235–256.

[24] E. 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.
[25] T. GUSTAFSSON, R. STENBERG, AND J. VIDEMAN, *On Nitsche's method for elastic contact problems*, SIAM J. Sci. Comput., 42 (2020), pp. B425-B446.

[26] R. HIPTMAIR AND C. PECHSTEIN, *Regular decompositions of vector fields - continuous, discrete, and structure-preserving*, Tech. Rep. 2019-18, Seminar for Applied Mathematics, ETH Zürich, Switzerland, 2019.

[27] P. HOUSTON, D. SCHÖTZAU, AND T. 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.

[28] P. HOUSTON, C. SCHWAB, AND E. SÜLI, *Discontinuous hp-finite element methods for advection-diffusion-reaction problems*, SIAM J. Numer. Anal., 39 (2002), pp. 2133-2163.

[29] O. KARAKASHIAN AND F. PASCAL, *A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems*, SIAM J. Numer. Anal., 41 (2003), pp. 2374-2399.

[30] ——, *Convergence of adaptive discontinuous Galerkin approximations of second-order elliptic problems*, SIAM J. Numer. Anal., 45 (2007), pp. 641-665.

[31] O. A. KARAKASHIAN AND F. PASCAL, *A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems*, SIAM J. Numer. Anal., 41 (2003), pp. 2374-2399.

[32] M. KARKULIK AND J. M. MELENK, *Local high-order regularization and applications to hp-methods*, Comput. Math. Appl., 70 (2015), pp. 1606-1639.

[33] N. KOPTEVA AND R. RANKIN, *Pointwise a posteriori error estimates for discontinuous Galerkin methods for singularly perturbed reaction-diffusion equations*, SIAM J. Numer. Anal., 61 (2023), pp. 1938-1961.

[34] J. MELENK AND T. WURZER, *On the stability of the boundary trace of the polynomial L2-projection on triangles and tetrahedra*, Comput. Math. Appl., 67 (2014), pp. 944-965.

[35] T. WARBURTON AND J. S. HESTHAVEN, *On the constants in hp-finite element trace inverse inequalities*, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 2765-2773.

[36] M. WHEELER, *An elliptic collocation-finite element method with interior penalties*, SIAM J. Numer. Anal., 15 (1978), pp. 152-161.

(Z.DONG) 1) INRIA, 48 RUE BARRAULT, 75647 PARIS, FRANCE AND 2) CERMICS, ECOLE DES PONTS, 77455 MARNE-LA-VALLÉE 2, FRANCE

Email address: zhaonan.dong@inria.fr

(E.H.GEORGOULIS) THE MAXWELL INSTITUTE FOR MATHEMATICAL SCIENCES & DEPARTMENT OF MATHEMATICS, SCHOOL OF MATHEMATICAL AND COMPUTER SCIENCES, HERIOT-WATT UNIVERSITY, EDINBURGH EH14 4AS, UNITED KINGDOM, AND DEPARTMENT OF MATHEMATICS, SCHOOL OF APPLIED MATHEMATICAL AND PHYSICAL SCIENCES, NATIONAL TECHNICAL UNIVERSITY OF ATHENS, ZOGRAFOU 15780, GREECE, AND IACM-FORTH, CRETE, GREECE

Email address: e.georgoulis@hw.ac.uk

(N.KOPTEVA) DEPARTMENT OF MATHEMATICS AND STATISTICS UNIVERSITY OF LIMERICK LIMERICK, V94 T9PX, IRELAND

Email address: Natalia.Kopteva@ul.ie