# Inexpensive polynomial-degree-robust equilibrated flux a posteriori estimates for isogeometric analysis

Gregor Gantner, Martin Vohralík

► To cite this version:

Gregor Gantner, Martin Vohralík. Inexpensive polynomial-degree-robust equilibrated flux a posteriori estimates for isogeometric analysis. Mathematical Models and Methods in Applied Sciences, In press, 34 (3), pp.477-522. 10.1142/S0218202524500076 . hal-03819048v3

HAL Id: hal-03819048

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

Submitted on 19 Nov 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
INEXPENSIVE POLYNOMIAL-DEGREE-ROBUST
EQUILIBRATED FLUX A POSTERIORI ESTIMATES
FOR ISOGEOMETRIC ANALYSIS

GREGOR GANTNER AND MARTIN VOHRALÍK

ABSTRACT. We consider isogeometric discretizations of the Poisson model problem, focusing on high polynomial degrees and strong hierarchical refinements. We derive a posteriori error estimates by equilibrated fluxes, i.e., vector-valued mapped piecewise polynomials lying in the $H(\text{div})$ space which appropriately approximate the desired divergence constraint. Our estimates are constant-free in the leading term, locally efficient, and robust with respect to the polynomial degree. They are also robust with respect to the number of hanging nodes arising in adaptive mesh refinement employing hierarchical B-splines. Two partitions of unity are designed, one with larger supports corresponding to the mapped splines, and one with small supports corresponding to mapped piecewise multilinear finite element hat basis functions. The equilibration is only performed on the small supports, avoiding the higher computational price of equilibration on the large supports or even the solution of a global system. Thus, the derived estimates are also as inexpensive as possible. An abstract framework for such a setting is developed, whose application to a specific situation only requests a verification of a few clearly identified assumptions. Numerical experiments illustrate the theoretical developments.

CONTENTS

<table><tr><td>1. Introduction</td><td>2</td></tr><tr><td>1.1. Available results</td><td>2</td></tr><tr><td>1.2. Outline</td><td>3</td></tr><tr><td>2. Model problem and its Galerkin discretization</td><td>3</td></tr><tr><td>3. Motivation and the new approach in a simplified spline setting</td><td>4</td></tr><tr><td>3.1. Standard equilibration in finite elements</td><td>4</td></tr><tr><td>3.2. A straightforward (expensive) generalization to IGA</td><td>5</td></tr><tr><td>3.3. Inexpensive equilibration in IGA</td><td>6</td></tr><tr><td>4. Hierarchical splines</td><td>8</td></tr><tr><td>4.1. Multivariate splines in the parameter domain (0, 1)<sup>d</sup></td><td>8</td></tr><tr><td>4.2. Hierarchical splines in the parameter domain (0, 1)<sup>d</sup></td><td>10</td></tr><tr><td>4.3. Hierarchical splines in the physical domain Ω</td><td>13</td></tr><tr><td>5. Partitions of unity and patchwise spaces</td><td>13</td></tr><tr><td>5.1. Partitions of unity based on hierarchical splines</td><td>13</td></tr><tr><td>5.2. Patchwise Sobolev spaces</td><td>15</td></tr><tr><td>5.3. Patchwise discrete subspaces</td><td>16</td></tr><tr><td>6. Abstract assumptions</td><td>18</td></tr><tr><td>7. Inexpensive equilibration</td><td>21</td></tr><tr><td>7.1. Inexpensive (lowest-order) residual lifting on the large patches ω<sub>a</sub></td><td>21</td></tr></table>

<footnote>This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 647134 GATIPOR). In addition, the first author has been supported by the Austrian Science Fund (FWF) under grant P29096 and J4379-N.</footnote>
<table>
  <tbody>
    <tr>
      <td>7.2. Projection operators for general bi-Lipschitz mappings F</td>
      <td>22</td>
    </tr>
    <tr>
      <td>7.3. Equilibrated flux on the small patches ω<sub>b</sub></td>
      <td>22</td>
    </tr>
    <tr>
      <td>7.4. Equilibrated flux on the large patches and the final equilibrated flux</td>
      <td>23</td>
    </tr>
    <tr>
      <td>8. A posteriori error estimates</td>
      <td>24</td>
    </tr>
    <tr>
      <td>8.1. Reliability</td>
      <td>25</td>
    </tr>
    <tr>
      <td>8.2. Efficiency</td>
      <td>26</td>
    </tr>
    <tr>
      <td>9. Numerical experiments</td>
      <td>31</td>
    </tr>
    <tr>
      <td>Appendix A. Proof of the broken polynomial extension property</td>
      <td>33</td>
    </tr>
    <tr>
      <td>References</td>
      <td>37</td>
    </tr>
  </tbody>
</table>

# 1. INTRODUCTION

A posteriori error estimates for standard finite element methods (FEM) are nowadays well understood [AO00, Rep08, Ver13]. On the one hand, they allow to assess the quality of the computed approximation, and, on the other hand, they indicate where to refine the underlying mesh of the computational domain. Among the existing error estimators, those based on equilibrated fluxes [PS47, LL83, DM99, BS08] have the advantage that they provide a *guaranteed* upper bound for the approximation error which is constant-free in the leading term. The estimator from [DM99, BS08] even turns out to be robust with respect to the polynomial degree [BPS09], i.e., it provides also a lower bound with an efficiency constant that does *not* depend on the polynomial degree. This result has been recently generalized in several directions; see [EV15, EV20] and the references therein. We also mention [DEV16], which covers standard FEM on triangular/rectangular meshes with hanging nodes, and its generalization to arbitrarily many hanging nodes [ESV17a].

In this work, we aim to generalize this result to isogeometric analysis (IGA) [HCB05, BBdVC+06, CHB09]. The central idea of IGA is to use the same ansatz functions for the representation of the problem geometry in computer-aided design (CAD) and for the discretization of the partial differential equation (PDE). While the CAD standard for spline representation in a multivariate setting relies on tensor-product B-splines, several extensions of the B-spline model have emerged that allow for adaptive refinement, e.g., hierarchical splines [VGJS11, GJS12], (analysis-suitable) T-splines [SZBN03, SLSH12, BdVBSV13], or LR-splines [DLP13, JKD14]; see also [JRK15, HKMP17, BGG+22] for a comparison of these approaches.

1.1. **Available results.** To steer an adaptive refinement, rigorous a posteriori error estimators have been developed. Assuming a certain admissibility condition of the employed meshes, the works [BG16b, GP20] generalize the weighted residual error estimator from standard FEM to IGA with hierarchical splines and analysis-suitable T-splines, respectively. It has even been proved that the corresponding adaptive algorithms converge at optimal algebraic rate with respect to the number of mesh elements [BG17, GHP17, GP20], see also the recent overview article [BGG+22]. However, the reliability and efficiency constants depend on the employed polynomial degree, which is indeed witnessed in numerical experiments, see, e.g., [BGG+22] for hierarchical splines. Similarly, for the hierarchical spline space from [BG16a], the work [BG18] derives a weighted residual estimator being the sum of indicators associated to basis functions instead of elements. It is shown to be reliable for arbitrary hierarchical meshes, with an unknown reliability constant that again particularly depends on the used polynomial degree *p*.

In the spirit of [Rep99, Rep00a, Rep00b], the works [KT15a, KT15b, Mat18] present guaranteed fully computable upper bounds of the approximation error for tensor-product
splines and hierarchical splines, respectively. A second estimate is also derived, giving a lower bound of the error. However, to compute these so-called functional-type estimators, a global minimization problem may need to be solved or an $H(\text{div})$ flux not in the equilibrium with the load may be employed, and the efficiency of the upper bound, as well as the reliability of the lower bound, are theoretically unclear. In the recent work [TCHM19], a non-$H(\text{div})$ approximation of an equilibrated flux is constructed for tensor-product splines in extension of the concepts of [LL83]. It only requires to solve locally, on the knot spans, a pair of a low-order problem for the normal fluxes together with a high-order problem for the equilibrated flux approximation. A generalization to hierarchical splines is briefly sketched and a corresponding numerical example is provided. While this yields a fully computable and approximately guaranteed upper bound on the error, efficiency of the resulting estimator is unclear.

## 1.2. Outline.
The present manuscript is organized as follows. We introduce the model problem and its general Galerkin discretization in Section 2. We then motivate and introduce our new approach in the simplified spline setting in Section 3. Sections 4 and 5 subsequently contain a detailed description of the discretization space formed by hierarchical B-splines. This technical description might be skipped in a first reading, since all our results are only based on a couple of abstract assumptions collected in Section 6, where we also verify them in the hierarchical B-splines setting. Section 7 describes in detail our inexpensive flux equilibration in the general context. The a posteriori error estimates are then derived in Section 8, where we prove their reliability and (local) efficiency. Numerical experiments illustrating the theoretical findings are collected in Section 9. Finally, Appendix A proves the broken polynomial extension property that is central for the polynomial-degree robustness in the present setting.

# 2. MODEL PROBLEM AND ITS GALERKIN DISCRETIZATION

Let $\Omega$ be an open bounded connected Lipschitz domain in $\mathbb{R}^d$ with $d = 2, 3$. We consider the Poisson model problem with homogeneous Dirichlet data of finding $u : \Omega \to \mathbb{R}$ such that

$$ 
\begin{array}{ll}
-\Delta u = f & \text{in } \Omega, \\
u = 0 & \text{on } \partial\Omega,
\end{array}
\qquad (1) $$

where $f \in L^2(\Omega)$ is a given source term. On an arbitrary subset $\omega \subset \bar{\Omega}$, let $(\cdot, \cdot)_\omega := \int_\omega (\cdot)(\cdot) dx$ and $\| \cdot \|_\omega$ denote the $L^2$-scalar product and the corresponding norm, respectively; we also denote by $\| \cdot \|_{\infty, \omega}$ the $L^\infty$-norm. Then, the weak formulation of problem (1) consists in finding $u \in H_0^1(\Omega)$ such that

$$ (\nabla u, \nabla v)_\Omega = (f, v)_\Omega \quad \text{for all } v \in H_0^1(\Omega). \qquad (2) $$

Let $V_h$ be a finite-dimensional subspace of $H_0^1(\Omega)$. The corresponding Galerkin approximation is to find $u_h \in V_h$ with

$$ (\nabla u_h, \nabla v_h)_\Omega = (f, v_h)_\Omega \quad \text{for all } v_h \in V_h. \qquad (3) $$

In this work, we will choose $V_h$ as the space of mapped hierarchical splines. To this end, we assume that $\Omega$ can be parametrized over $\hat{\Omega} := (0, 1)^d$ via a bi-Lipschitz mapping $F: \hat{\Omega} \to \Omega$ with positive Jacobian determinant, i.e., $\Omega = F(\hat{\Omega})$, $F \in W^{1,\infty}(\hat{\Omega})$, $F^{-1} \in W^{1,\infty}(\Omega)$, and $\det(DF) > 0$. Note that $F$ and its inverse $F^{-1}$ can both be continuously extended to Lipschitz continuous functions on $\bar{\Omega}$ and $\bar{\Omega}$, respectively.
### 3. MOTIVATION AND THE NEW APPROACH IN A SIMPLIFIED SPLINE SETTING

In this section, we recall equilibration for standard $C^0$ finite elements, discuss difficulties for a generalization to IGA, and introduce our new approach. In order to explain the main ideas as clearly as possible, we will consider a strongly simplified setting. In particular, we suppose in this section that $F$ is just the identity, so that $\Omega = \hat{\Omega} = (0, 1)^d$, and that $\mathcal{T}_h$ is a uniform tensor-product mesh of $\Omega$ consisting of elements $K$ being all $d$-rectangular parallelepipeds. Let $\mathbb{Q}^p(\mathcal{T}_h)$ denote all $\mathcal{T}_h$-piecewise tensor-product polynomials of fixed degree $p := (p, \dots, p)$ in each component, with $p \ge 1$. We assume additionally that the right-hand side $f$ in (1) is a $\mathcal{T}_h$-piecewise tensor-product polynomial of degree $p-1 := p - (1, \dots, 1)$, i.e.,

$$f \in \mathbb{Q}^{p-1}(\mathcal{T}_h).$$

3.1. **Standard equilibration in finite elements.** In finite elements, the ansatz space $V_h$ from (3) is the subspace of $\mathbb{Q}^p(\mathcal{T}_h)$ formed by functions which are continuous over the mesh interfaces and zero on the boundary of $\Omega$, i.e.,

$$V_h := \{v_h \in \mathbb{Q}^p(\mathcal{T}_h) \cap C^0(\Omega) : v_h|_{\partial\Omega} = 0\} = \mathbb{Q}^p(\mathcal{T}_h) \cap H_0^1(\Omega). \quad (4)$$

Let $\mathcal{V}_h$ collect the vertices of the mesh $\mathcal{T}_h$. For each $\mathbf{a} \in \mathcal{V}_h$, we define the associated continuous piecewise bilinear hat function

$$\psi_{\mathbf{a}} \in \mathbb{Q}^1(\mathcal{T}_h) \cap C^0(\Omega) \quad (5)$$

taking value 1 in the vertex $\mathbf{a}$ and 0 in all other vertices from $\mathcal{V}_h$. We denote the interior of its support by

$$\omega_{\mathbf{a}} := \text{int}(\text{supp}(\psi_{\mathbf{a}})). \quad (6)$$

Note that it corresponds to the union of all (closed) elements in $\mathcal{T}_h$ containing $\mathbf{a}$, which we denote by

$$\mathcal{T}_{\mathbf{a}} := \{T \in \mathcal{T}_h : T \subseteq \text{supp}(\psi_{\mathbf{a}})\}.$$

Crucially, the $\psi_{\mathbf{a}}$ form a partition of unity in that

$$\sum_{\mathbf{a} \in \mathcal{V}_h} \psi_{\mathbf{a}} = 1 \quad \text{in } \Omega.$$

The hat functions $\psi_{\mathbf{a}}$ and their supports $\text{supp}(\psi_{\mathbf{a}})$ are illustrated for $d = 2$ in Figure 1, left.

For each vertex $\mathbf{a} \in \mathcal{V}_h$, we further introduce the space

$$H_0(\text{div}, \omega_{\mathbf{a}}) := \begin{cases} \{v \in \mathcal{H}(\text{div}, \omega_{\mathbf{a}}) : v \cdot n_{\omega_{\mathbf{a}}} = 0 \text{ on } \partial\omega_{\mathbf{a}}\} & \text{if } \psi_{\mathbf{a}} \in H_0^1(\Omega), \\ \{v \in \mathcal{H}(\text{div}, \omega_{\mathbf{a}}) : v \cdot n_{\omega_{\mathbf{a}}} = 0 \text{ on } \partial\omega_{\mathbf{a}} \cap (\psi_{\mathbf{a}})^{-1}(\{0\})\} & \text{else,} \end{cases}$$

where $n_{\omega_{\mathbf{a}}}$ denotes the outer normal vector on $\partial\omega_{\mathbf{a}}$ and $v \cdot n_{\omega_{\mathbf{a}}}$ is understood in the appropriate weak sense. With the usual broken (elementwise) Raviart-Thomas space of order $\tilde{p} := (\tilde{p}, \dots, \tilde{p}) := p + 1$

$$RT^{\tilde{p}}(\mathcal{T}_a) := \begin{cases} \mathbb{Q}^{\tilde{p}+(1,0)}(\mathcal{T}_a) \times \mathbb{Q}^{\tilde{p}+(0,1)}(\mathcal{T}_a) & \text{if } d = 2, \\ \mathbb{Q}^{\tilde{p}+(1,0,0)}(\mathcal{T}_a) \times \mathbb{Q}^{\tilde{p}+(0,1,0)}(\mathcal{T}_a) \times \mathbb{Q}^{\tilde{p}+(0,0,1)}(\mathcal{T}_a) & \text{if } d = 3, \end{cases} \quad (7)$$

we also define the discrete subspace

$$V_h^a := RT^{\tilde{p}}(\mathcal{T}_a) \cap H_0(\text{div}, \omega_{\mathbf{a}}). \quad (8)$$
Then, the local equilibrated fluxes

$$ \sigma_h^a := \underset{\substack{v_h \in V_h^a \\ \nabla \cdot v_h = f \psi_a - \nabla u_h \cdot \nabla \psi_a}}{\operatorname{argmin}} \|v_h + \psi_a \nabla u_h\|_{\omega_a} \quad (9) $$

give rise to a global equilibrated flux

$$ \sigma_h := \sum_{a \in \mathcal{V}_h} \sigma_h^a \quad \text{with} \quad \nabla \cdot \sigma_h = f. $$

Following [BPS09, EV15, EV20], which rely on the tools from [CDD08, CM10, DGS12], one can show that the corresponding a posteriori estimator $||\sigma_h + \nabla u_h||_\Omega$ constitutes a guaranteed upper bound as well as a p-robust lower bound for the discretization error, i.e.,

$$ ||\nabla(u_h - u)||_\Omega \le ||\sigma_h + \nabla u_h||_\Omega \le C_{\text{eff}} ||\nabla(u_h - u)||_\Omega, \quad (10) $$

with a constant $C_{\text{eff}} > 0$ depending only on the shape-regularity of the mesh $\mathcal{T}_h$ but not on the polynomial degree $p$. The efficiency bound holds even locally on all mesh elements. This result has been extended to general rectangular meshes with an arbitrary number of hanging nodes in [DEV16, ESV17a], see also Remark 7.6 below.

## 3.2. A straightforward (expensive) generalization to IGA.
Compared to standard finite elements (4), the IGA/spline space is the subspace of $\mathbb{Q}^p(\mathcal{T}_h)$ formed by functions which are continuous over the mesh interfaces including their derivatives up to order $p-m$, i.e.,

$$ V_h := \{v_h \in \mathbb{Q}^p(\mathcal{T}_h) \cap C^{p-m}(\Omega) : v_h|_{\partial\Omega} = 0\}. \quad (11) $$

Here $m$ takes values between 1 and $p$. Examples for $m=1$, leading to the maximal smoothness $p-1$, and $p=1,3$ are given in Figure 1; note that for $m=p$, (11) and (4) coincide. A support of $\psi_a$ for the maximal smoothness $p-1$ and $p=5$ is then shown in the left part of Figure 2. Alternatively, in place of (11), one can write $V_h = \mathbb{S}^p(\mathcal{K}_h) \cap H_0^1(\Omega)$, where $\mathbb{S}^p(\mathcal{K}_h)$ is the spline space with the smoothness $p-m$ encoded in the global knot vector $\mathcal{K}_h$, see Section 4.1 below.

Whenever the requested smoothness $p-m$ is greater than 0, the finite element hat functions $\psi_a$ from (5) do not belong to the IGA space $V_h$ given by (11), being merely continuous. A straightforward adaptation of the flux reconstruction of the previous section would be to employ instead appropriately scaled $\psi_a$ from the spline space $\mathbb{Q}^p(\mathcal{T}_h) \cap C^{p-m}(\Omega) = \mathbb{S}^p(\mathcal{K}_h)$, still leading to the partition of unity

$$ \sum_{a \in \mathcal{V}_h} \psi_a = 1 \quad \text{in } \Omega, $$

where now $\mathcal{V}_h$ is a suitable node set rather than the set of vertices of the mesh $\mathcal{T}_h$. Standard B-splines, see Section 4.1, constitute such a partition of unity. The major inconvenience of such a construction is that the supports $\omega_a$ of $\psi_a$ (cf. (6)), which appear in the minimization (9), are much larger for splines. In particular, for the maximal smoothness $p-1$, the support of the corresponding B-splines of degree $p$ consists of $(p+1)^d$ elements, see Figure 2. This would make the dimension of the space $V_h^a$ from (8) and consequently the size of the linear system in (9) grow steeply with $p$, making the equilibration very expensive, see Figure 2 for illustration. In addition, in the construction (9), one needs $\psi_a \nabla u_h \in \mathbf{R}T\tilde{p}(\mathcal{T}_a)$ and $f\psi_a - \nabla u_h \cdot \nabla \psi_a \in \mathbb{Q}\tilde{p}(\mathcal{T}_a)$ to avoid all oscillation terms. In the finite element context of Section 3.1, this is satisfied for the choice $\tilde{p} = p+1$, but it would request $\tilde{p} = 2p$ here. Moreover, it is not clear whether (10) would still hold with a constant $C_{\text{eff}}$ independent of the polynomial degree $p$ of the considered B-splines.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. B-splines of degree $p = 1, 3$ and maximal smoothness $p-1$ on a two-dimensional uniform tensor mesh $\mathcal{T}_h$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2. A straightforward (expensive) extension of the FEM equilibration to IGA with B-splines of degree $p = 5$ and maximal smoothness $p-1$ on a two-dimensional uniform tensor mesh $\mathcal{T}_h$.</figcaption></figure>

### 3.3. Inexpensive equilibration in IGA.
Our main idea is to replace the straightforward approach of Section 3.2 by a construction that does not request the expensive solution of (9) employing the large spaces $V_h^a$ of (8) defined over the large spline supports $\omega_a$. We design a two-stage procedure where: first, lowest-order $Q^1(\mathcal{T}_a) \cap C^0(\omega_a)$ scalar problems are solved on each large patch $\mathcal{T}_a$; second, systems of the form (9), but only over the finite element hat basis functions supports $\omega_b$ with $2^d$ mesh elements, are solved. This is schematically viewed in Figure 3. Similar $p$-robust constructions were designed in [SMPZ08] in the context of an additive Schwarz smoother/preconditiner (just one global lowest-order solve and then $p$-order local patch problems) or in [ESV17a, ESV17b, PV22] for a posteriori error estimates (allowing for fast equilibration with $H^{-1}(\Omega)$ right-hand sides/arbitrary coarsening in parabolic time stepping/inexpensive estimates of the algebraic error).

More precisely, we first, on each (large) patch subdomain $\omega_a$ with local mesh $\mathcal{T}_a$ solve the primal lowest-order problem: find

$$r_h^a \in V_h^a := \begin{cases} \{v_h \in Q^1(\mathcal{T}_a) \cap C^0(\omega_a) : (v_h, 1)_{\omega_a} = 0\} & \text{if } \psi_a \in H_0^1(\Omega), \\ \{v_h \in Q^1(\mathcal{T}_a) \cap C^0(\omega_a) : v_h = 0 \text{ on } \partial\omega_a \setminus \psi_a^{-1}(\{0\})\} & \text{else} \end{cases}$$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Inexpensive equilibration in IGA in the simplified setting of uniform tensor-product mesh, $p = 5$ and maximal smoothness $p - 1$.</figcaption></figure>

such that

$$ (\nabla r_h^a, \nabla v_h)_{\omega_a} = (f, v_h \psi_a)_{\omega_a} - (\nabla u_h, \nabla (v_h \psi_a))_{\omega_a} \quad \text{for all } v_h \in V_h^a, \quad (12) $$

see Figure 3, steps 1)–3). This yields the scalar-valued lifting $r_h^a$ of the $\psi_a$-weighted residual of (3) with respect to (2). Second, we consider the finite element hat basis functions $\psi_b$ of vertices $b$ of the local mesh $\mathcal{T}_a$ and solve the dual high-order problems

$$ \sigma_h^{a,b} := \underset{v_h \in V_h^{a,b}}{\operatorname{argmin}} \quad \|v_h + \psi_b(\psi_a \nabla u_h + \nabla r_h^a)\|_{\omega_b}, \quad (13) \\ \nabla \cdot v_h = f \psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b) - \nabla r_h^a \cdot \nabla \psi_b $$

where $V_h^{a,b}$ is defined similarly to (8), with

$$ H_0(\operatorname{div}, \omega_b) := \begin{cases} \{v \in \mathcal{H}(\operatorname{div}, \omega_b) : v \cdot n_{\omega_b} = 0 \text{ on } \partial\omega_b\} & \text{if } \psi_a \psi_b \in H_0^1(\Omega), \\ \{v \in \mathcal{H}(\operatorname{div}, \omega_b) : v \cdot n_{\omega_b} = 0 \text{ on } \partial\omega_b \\ \qquad \cap (\psi_a \psi_b)^{-1}(\{0\})\} & \text{else,} \end{cases} $$

see Figure 3, steps 4)–5). As in Section 3.2, we need to increase the equilibration polynomial degree $\tilde{p}$; actually, owing to the presence of the finite element hat basis function $\psi_b$, we will need $\tilde{p} = 2p + 1$ and not just $\tilde{p} = 2p$ to satisfy $\psi_b(\psi_a \nabla u_h + \nabla r_h^a) \in \mathcal{R}T^{\tilde{p}}(\mathcal{T}_a)$ and $f\psi_a\psi_b - \nabla u_h \cdot \nabla (\psi_a\psi_b) - \nabla r_h^a \cdot \nabla \psi_b \in \mathbb{Q}^{\tilde{p}}(\mathcal{T}_a)$ to avoid all oscillation terms. All these local problems can be solved independently one from each other, allowing for efficient parallelization. Moreover, for every two same patch geometries, one matrix assembly and factorization is sufficient.

To finish, we sum the contributions $\sigma_h^{a,b}$ from (13) as

$$ \sigma_h^a := \sum_{b \in V_h^a} \sigma_h^{a,b} \quad (14) $$

and form the final equilibrated flux as

$$ \sigma_h := \sum_{a \in V_h} \sigma_h^a, \quad (15) $$

which yields a fully computable guaranteed (constant-free in the leading term) upper bound on the unknown error similarly to (10), see Proposition 8.2 for details. This bound is also locally and globally efficient, see respectively Propositions 8.4 and 8.6. The involved constant $C_{\text{eff}}$ is robust with respect to the polynomial degree $p$ of the used hierarchical splines and the number of hanging nodes of the underlying hierarchical
meshes, but theoretically not with respect to the smoothness of the splines and the number of patches $\omega_a$ overlapping in a point (see (84a)) (in our numerical experiments, though, $C_{\text{eff}}$ shows robustness with respect to the smoothness and the number of overlapping patches). Full details of the equilibration are given in Definitions 7.1, 7.2, and 7.5 below; Figure 6 illustrates the procedure on hierarchical meshes.

## 4. HIERARCHICAL SPLINES

In this section, we first describe the piecewise polynomial space of multivariate splines in the parameter domain $\hat{\Omega}$. We then introduce its hierarchical extension, covering highly graded local mesh refinement. Finally, the space $V_h$ from (3) will be given by the transformation of the latter one by the mapping $\mathbf{F}$.

### 4.1. Multivariate splines in the parameter domain $(0, 1)^d$.
We recall here the standard definition of multivariate splines in the parameter domain $(0, 1)^d$; for a detailed introduction, we refer, e.g., to [dB86, dB01, Sch07].

Let the integer $p \ge 1$ be a fixed positive polynomial degree and let

$$K_h = (K_{1(h)}, \dots, K_{d(h)})$$

be a fixed $d$-dimensional vector of $p$-open knot vectors, i.e., for each spatial dimension $1 \le i \le d$,

$$K_{i(h)} = (t_{i(h),0}, \dots, t_{i(h),N_{i(h)}+p})$$

is a $p$-open knot vector in $[0, 1]$, which means that

$$0 = t_{i(h),0} = \dots = t_{i(h),p} < t_{i(h),p+1} \le t_{i(h),p+2} < \dots < t_{i(h),N_{i(h)}} = \dots = t_{i(h),N_{i(h)}+p} = 1.$$

Moreover, we assume that each of the interior knots $t_{i(h),j} \in (0, 1)$ appears at most with multiplicity $p$. By definition, the boundary knots 0 and 1 have multiplicity $p+1$. An example is given in Figure 4, with polynomial degree $p=2$, number of knots minus $p$ minus one $= N_{i(h)} = 9$ (this will later correspond to the dimension of the B-splines space), and a varying multiplicity.

We define the resulting one-dimensional meshes

$$\hat{T}_{i(h)} := \{[t_{i(h),j-1}, t_{i(h),j}] : j \in \{1, \dots, N_{i(h)} + p\} \wedge t_{i(h),j-1} < t_{i(h),j}\}, \quad i = 1, \dots, d,$$

as well as the resulting tensor mesh

$$\hat{T}_h := \{\hat{T}_1 \times \dots \times \hat{T}_d : \hat{T}_i \in \hat{T}_{i(h)} \text{ for all } i \in \{1, \dots, d\}\}.$$

By $\mathbb{S}^p(K_{i(h)})$, we denote the set of all corresponding univariate splines, i.e., the set of all $\hat{T}_{i(h)}$-piecewise univariate polynomials of degree $p$ that are $p - \#t_{i(h),j}$-times continuously differentiable at any interior knot $t_{i(h),j}$, where $\#t_{i(h),j}$ denotes the corresponding multiplicity of $t_{i(h),j}$ within $K_{i(h)}$. Assuming for example that all interior multiplicities are equal to $p$, $\mathbb{S}^p(K_{i(h)})$ is just the space of all $\hat{T}_{i(h)}$-piecewise univariate polynomials of degree $p$ that are merely continuous, i.e., $C^0(0, 1)$, cf. (4). At the other extreme, when all interior multiplicities are equal to 1, then $\mathbb{S}^p(K_{i(h)})$ is the space of $\hat{T}_{i(h)}$-piecewise univariate polynomials of degree $p$ with continuous derivatives up to order $p-1$, cf. (11) with $m=1$. We refer again to Figure 4 for an example.

A basis of $\mathbb{S}^p(K_{i(h)})$, of dimension $N_{i(h)}$, is given by the set of B-splines

$$\{B_{i(h),j,p} : j \in \{1, \dots, N_{i(h)}\}\},$$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. The B-splines of degree $p = 2$ corresponding to the knot vector $\mathcal{K} = (0, 0, 0, 1/6, 2/6, 3/6, 4/6, 4/6, 5/6, 1, 1, 1)$ are depicted. They are at least $C^1$ at the knots $1/6, 2/6, 3/6$, and $5/6$, and at least $C^0$ at $4/6$.</figcaption></figure>

where the B-splines $B_{i(h),j,p}$ are recursively defined for all points $t \in (0, 1)$ via

$$ B_{i(h),j,0}(t) := \begin{cases} 1 & \text{if } t_{i(h),j} \le t < t_{i(h),j+1}, \\ 0 & \text{else,} \end{cases} \qquad (16a) $$

and

$$ B_{i(h),j,p}(t) := \frac{t - t_{i(h),j}}{t_{i(h),j+p} - t_{i(h),j}} B_{i(h),j,p-1}(t) + \frac{t_{i(h),j+p+1} - t}{t_{i(h),j+p+1} - t_{i(h),j+1}} B_{i(h),j+1,p-1}(t) \quad (16b) $$

with the formal convention that $\cdot/0 := 0$. Figure 4 gives an illustrative example, where the $N_{i(h)} = 9$ basis functions are depicted. It is easy to see that the support is

$$ \text{supp}(B_{i(h),j,p}) = [t_{i(h),j}, t_{i(h),j+p+1}] \qquad (17) $$

and we also remark that

$$ 0 \le B_{i(h),j,p} \le 1. \qquad (18) $$

We abbreviate $\mathbf{p} := (p, \dots, p) \in \mathbb{N}^d$. The space $\mathbb{S}^p(\mathcal{K}_h)$ of multivariate splines is defined as tensor-product of the univariate spline spaces. Note that each function in $\mathbb{S}^p(\mathcal{K}_h)$ is a $\hat{T}_h$-piecewise multivariate polynomial of degree $\mathbf{p}$. Again, assuming for example that all interior knots have multiplicity $p$, $\mathbb{S}^p(\mathcal{K}_h)$ is just the space of all continuous $\hat{T}_h$-piecewise multivariate polynomials of degree $\mathbf{p}$, i.e., up to boundary conditions, (4), while multiplicity 1 of all interior knots yields the space of all $\hat{T}_h$-piecewise multivariate polynomials of degree $\mathbf{p}$ with continuous derivatives up to order $p-1$, i.e., up to boundary conditions, (11) for $m = 1$. Clearly, the set of tensor-products of the univariate B-splines

$$ \{B_{h,j,p} : j \in \Pi_{i=1}^d \{1, \dots, N_{i(h)}\}\} \quad \text{with} \quad B_{h,j,p}(t) := \Pi_{i=1}^d B_{i(h),j,i,p}(t_i), \qquad (19) $$

provides a basis of $\mathbb{S}^p(\mathcal{K}_h)$.
4.2. **Hierarchical splines in the parameter domain $(0, 1)^d$.** We now introduce hierarchical splines, which are defined on a hierarchical mesh and are essentially coarse splines on coarse mesh elements and fine splines on fine mesh elements. For a detailed introduction, we refer, e.g., to the seminal work [VGJS11].

Let $p$ be a fixed positive polynomial degree as above and let $\mathcal{K}_0 = (\mathcal{K}_{1(0)}, \dots, \mathcal{K}_{d(0)})$ be a fixed initial $d$-dimensional vector of $p$-open knot vectors in $[0, 1]$ with interior multiplicities less than or equal to $p$, as in Section 4.1. Recall that we consider $d = 2, 3$. We set $\mathcal{K}_{\text{uni}(0)} := \mathcal{K}_0$ and recursively define $\mathcal{K}_{\text{uni}(\ell+1)}$ for $\ell \in \mathbb{N}_0$ as the uniform $h$-refinement of $\mathcal{K}_{\text{uni}(\ell)}$ with fixed multiplicity $m \in \mathbb{N}$, i.e., obtained by inserting the knot $(t_{i(\text{uni}(\ell)),j-1} + t_{i(\text{uni}(\ell)),j})/2$ to the knots $\mathcal{K}_{i(\text{uni}(\ell))}$ with multiplicity $m$ whenever $t_{i(\text{uni}(\ell)),j-1} < t_{i(\text{uni}(\ell)),j}$. We use analogous notation as in Section 4.1, replacing the index $h$ by $\text{uni}(\ell)$, e.g., we write $\hat{\mathcal{T}}_{\text{uni}(\ell)}$ for the induced tensor mesh. We stress that these spline spaces are nested in the sense that

$$ \mathbb{S}^p(\mathcal{K}_{\text{uni}(\ell)}) \subset \mathbb{S}^p(\mathcal{K}_{\text{uni}(\ell+1)}) \subset C^0((0, 1)^d), \quad (20) $$

where the last relation follows from the assumption that multiplicities of interior knots is less than or equal to $p$.

**Remark 4.1.** As a matter of fact, the definition of the uniform refinements $\mathcal{K}_{\text{uni}(\ell)}$ allows for a certain flexibility without changing the results of the present manuscript. For instance, instead of the natural dyadic uniform refinement, one could use n-adic refinement, i.e., insert the knots $(t_{i(\text{uni}(\ell)),j-1} + t_{i(\text{uni}(\ell)),j})/n, \dots, (n-1)(t_{i(\text{uni}(\ell)),j-1} + t_{i(\text{uni}(\ell)),j})/n$ to the knots $\mathcal{K}_{i(\text{uni}(\ell))}$ with multiplicity $m$ whenever $t_{i(\text{uni}(\ell)),j-1} < t_{i(\text{uni}(\ell)),j}$.

We say that a set

$$ \hat{\mathcal{T}}_h \subseteq \bigcup_{\ell \in \mathbb{N}_0} \hat{\mathcal{T}}_{\text{uni}(\ell)} $$

is a hierarchical mesh if it is a partition of $[0, 1]^d$ in the sense that $\bigcup \hat{\mathcal{T}}_h = [0, 1]^d$, where the intersection of two different elements $\hat{T} \neq \hat{T}'$ with $\hat{T}, \hat{T}' \in \hat{\mathcal{T}}_h$ has ($d$-dimensional) measure zero. Since $\hat{\mathcal{T}}_{\text{uni}(\ell)} \cap \hat{\mathcal{T}}_{\text{uni}(\ell')} = \emptyset$ for $\ell, \ell' \in \mathbb{N}_0$ with $\ell \neq \ell'$, we can define for any element $\hat{T} \in \hat{\mathcal{T}}_h$,

$$ \text{level}(\hat{T}) := \ell \in \mathbb{N}_0 \text{ with } \hat{T} \in \hat{\mathcal{T}}_{\text{uni}(\ell)}. $$

For an illustrative example of a hierarchical mesh, see Figure 5. In particular, any uniformly refined tensor mesh $\hat{\mathcal{T}}_{\text{uni}(\ell)}$ with $\ell \in \mathbb{N}_0$ is a hierarchical mesh.

For a hierarchical mesh $\hat{\mathcal{T}}_h$, we define a corresponding nested sequence $(\hat{\Omega}_h^\ell)_{\ell \in \mathbb{N}_0}$ of closed subsets of $[0, 1]^d$ as

$$ \hat{\Omega}_h^\ell := \bigcup_{\ell' \ge \ell} (\hat{T}_h \cap \hat{\mathcal{T}}_{\text{uni}(\ell')}), $$

i.e., $\hat{\Omega}_h^\ell$ covers all elements with level greater than or equal to $\ell$. Note that there exists a minimal integer $L_h$ such that $\hat{\Omega}_h^\ell = \emptyset$ for all $\ell \ge L_h$. It holds that

$$ \hat{\mathcal{T}}_h = \bigcup_{\ell \in \mathbb{N}_0} \{\hat{T} \in \hat{\mathcal{T}}_{\text{uni}(\ell)} : \hat{T} \subseteq \hat{\Omega}_h^\ell \wedge \hat{T} \not\subseteq \hat{\Omega}_h^{\ell+1}\}. \quad (21) $$

In the literature, one usually assumes that the sequence $(\hat{\Omega}_h^\ell)_{\ell \in \mathbb{N}_0}$ is given and the corresponding hierarchical mesh is defined via (21).
<figure><img src="image_6.png" /><figcaption>FIGURE 5. A two-dimensional hierarchical mesh $\hat{T}_h$ with level of all elements less than 4. Levels 0 to 3 are respectively highlighted in white, light grey, grey, and dark grey, also denoting the corresponding domains $\hat{\Omega}_h^0 \setminus \hat{\Omega}_h^1$, $\hat{\Omega}_h^1 \setminus \hat{\Omega}_h^2$, $\hat{\Omega}_h^2 \setminus \hat{\Omega}_h^3$, $\hat{\Omega}_h^3 \setminus \hat{\Omega}_h^4$, $\hat{\Omega}_h^4 = \emptyset$.</figcaption></figure>

We introduce the hierarchical basis

$$ \{B_a : \mathbf{a} \in \mathcal{I}_h\} \quad \text{with} \quad \mathcal{I}_h := \bigcup_{\ell \in \mathbb{N}_0} \left\{ \begin{array}{l} (\text{uni}(\ell), \mathbf{j}, \mathbf{p}) : \mathbf{j} \in \Pi_{i=1}^d \{1, \dots, N_i(\text{uni}(\ell))\} \\ \qquad \wedge \text{supp}(B_{\text{uni}(\ell), \mathbf{j}, \mathbf{p}}) \subseteq \hat{\Omega}_h^\ell \wedge \text{supp}(B_{\text{uni}(\ell), \mathbf{j}, \mathbf{p}}) \not\subseteq \hat{\Omega}_h^{\ell+1} \end{array} \right\}, \quad (22) $$

where we recall the definition (19) of a multivariate B-spline. Figure 6 gives an illustrative example. Its elements are referred to as (multivariate) hierarchical B-splines. For $\mathbf{a} \in \mathcal{I}_h$, the level of the corresponding hierarchical B-spline is well defined

$$ \text{level}(B_a) := \ell \in \mathbb{N}_0 \quad \text{with} \quad \mathbf{a} = (\text{uni}(\ell), \mathbf{j}, \mathbf{p}). \quad (23) $$

It is easy to check that if $\hat{T}_h$ is a tensor mesh, and hence coincides with some $\hat{T}_{\text{uni}(\ell)}$, then the hierarchical basis and the standard tensor-product B-spline basis are the same. One can prove that the hierarchical B-splines are linearly independent; see, e.g., [VGJS11, Lemma 2]. They span the space of hierarchical splines

$$ S^p(\mathcal{K}_0, m, \hat{T}_h) := \operatorname{span}(\{B_a : \mathbf{a} \in \mathcal{I}_h\}). \quad (24) $$

The hierarchical basis and the mesh $\hat{T}_h$ are compatible in the following sense: For all $B_a$, $\mathbf{a} \in \mathcal{I}_h$, the corresponding support can be written as union of elements in $\hat{T}_{\text{uni(level}(B_a))}$, i.e.,

$$ \text{supp}(B_a) = \bigcup \{\hat{T} \in \hat{T}_{\text{uni(level}(B_a))} : \hat{T} \subseteq \text{supp}(B_a)\}. $$

Each such element $\hat{T} \in \hat{T}_{\text{uni(level}(B_a))}$ with $\hat{T} \subseteq \text{supp}(B_a) \subseteq \hat{\Omega}_h^{\text{level}(B_a)}$ satisfies that $\hat{T} \in \hat{T}_h$ or $\hat{T} \subseteq \hat{\Omega}_h^{\text{level}(B_a)+1}$. In either case, we see that $\hat{T}$ can be written as union of elements in
<figure><img src="image_7.png" /><figcaption>FIGURE 6. A two-dimensional hierarchical mesh $\hat{T}_h$ (of Figure 5) is depicted in black. Assuming that all interior knots have multiplicity 1, in (a), the support of three hierarchical B-splines of degree $p = 2$ is highlighted in blue; two hierarchical B-splines of level 0 (left and bottom), and one hierarchical B-spline of level 1 (right). The corresponding local tensor meshes $\hat{T}_a$ are indicated in grey. In (b), these three hierarchical B-splines are depicted. For $\bar{p} = 2$, these functions also belong (up to scaling) to the partition of unity $\{\psi_a \circ F : a \in V_h\}$ defined in Section 5.1 below. Finally, in (c), for each of them, three functions of the corresponding local partition of unity $\{\psi_b \circ F : b \in V_h^a\}$ defined in Section 5.1 below over the local tensor meshes are depicted.</figcaption></figure>

$\hat{T}_h$ with level greater or equal to level($B_a$). Altogether, we have that

$$ \mathrm{supp}(B_a) = \bigcup_{\ell \ge \mathrm{level}(B_a)} \{\hat{T} \in \hat{T}_h \cap \hat{T}_{\mathrm{uni}(\ell)} : \hat{T} \subseteq \mathrm{supp}(B_a)\}. \quad (25) $$

Moreover, $\mathrm{supp}(B_a)$ must contain at least one element of level($B_a$). Otherwise one would get the contradiction $\mathrm{supp}(B_a) \subseteq \hat{\Omega}_h^{\mathrm{level}(B_a)+1}$.

We now present the following characterization of hierarchical splines, which has been verified in [SM16, Section 3],

$$ \mathcal{S}^p(\mathcal{K}_0, m, \hat{T}_h) = \{\hat{v}_h : (0,1)^d \to \mathbb{R} : \hat{v}_h|_{(0,1)^d \setminus \hat{\Omega}_h^{\ell+1}} \in \mathcal{S}^p(\mathcal{K}_{\mathrm{uni}(\ell)})|_{(0,1)^d \setminus \hat{\Omega}_h^{\ell+1}} \text{ for all } \ell \in \mathbb{N}_0\}. \quad (26) $$
In particular, each hierarchical spline is a $\hat{T}_h$-piecewise tensor-product polynomial of degree $p$. Put into words, hierarchical splines are coarse splines on coarse mesh elements, and they are fine splines on fine mesh elements.

Finally, we say that a hierarchical mesh $\hat{T}_h$ is finer than another hierarchical mesh $\hat{T}_H$ if $\hat{T}_h$ is obtained from $\hat{T}_H$ via iterative dyadic bisection. Formally, this can be stated as $\hat{\Omega}_H^\ell \subseteq \hat{\Omega}_h^\ell$ for all $\ell \in \mathbb{N}_0$. In this case, (26) shows that the corresponding hierarchical spline spaces are nested, i.e.,

$$ \mathbb{S}^p(\mathcal{K}_0, m, \hat{T}_H) \subseteq \mathbb{S}^p(\mathcal{K}_0, m, \hat{T}_h). \qquad (27) $$

In particular, we see that

$$ \mathbb{S}^p(\mathcal{K}_{\text{uni}(0)}) \subseteq \mathbb{S}^p(\mathcal{K}_0, m, \hat{T}_h) \subseteq \mathbb{S}^p(\mathcal{K}_{\text{uni}(L_h-1)}). $$

4.3. **Hierarchical splines in the physical domain $\Omega$.** If $\hat{T}_h$ is a hierarchical mesh in the parameter domain $(0, 1)^d$, we set

$$ \mathcal{T}_h := \{ \mathbf{F}(\hat{T}) : \hat{T} \in \hat{T}_h \}, $$

where $\mathbf{F}$ is the bi-Lipschitz mapping from Section 2. In this context, $h \in L^\infty(\Omega)$ denotes the mesh size function defined by $h|_T := \text{diam}(T)$ for all $T \in \mathcal{T}_h$. Moreover, we set

$$ \mathbb{S}^p(\mathcal{K}_0, m, \mathcal{T}_h) := \{ \hat{v} \circ \mathbf{F}^{-1} : \hat{v} \in \mathbb{S}^p(\mathcal{K}_0, m, \hat{T}_h) \}. $$

The conforming ansatz space for the Galerkin discretization (3) is then defined as

$$ V_h := \mathbb{S}^p(\mathcal{K}_0, m, \mathcal{T}_h) \cap H_0^1(\Omega). \qquad (28) $$

**Remark 4.2.** In practice, the parametrization $\mathbf{F}$ is often given in terms of non-uniform rational B-splines (NURBS) along with corresponding polynomial degree $\mathbf{p}_F$ and global knot vector $\mathcal{K}_F$. To guarantee good approximation properties of $V_h$, in particular good a priori estimates, the functions in $\mathbb{S}^p(\mathcal{K}_0, m, \hat{T}_h)$ should have the same or lower smoothness as $\mathbf{F}$ along the knot lines on $\hat{\Omega}$ corresponding to $\mathcal{K}_F$. Such a requirement is not needed for the presented a posteriori error analysis to hold. Moreover, the form of the oscillation terms (77b) and (84c) for reliability and efficiency, respectively, suggests that the approximation property of $V_h$ has minor influence on the size of these terms, see also Remark 5.1.

## 5. PARTITIONS OF UNITY AND PATCHWISE SPACES

In this section, we prepare the necessary material for defining and analyzing our equilibrated fluxes in the IGA context later. We particularly design a partition of unity based on hierarchical B-splines and define continuous- and discrete-level local spaces.

5.1. **Partitions of unity based on hierarchical splines.** We now first construct a partition of unity on $\Omega$ consisting of hierarchical B-splines with the smallest-possible polynomial degree $\bar{p}$ but sufficient smoothness to be contained in $\mathbb{S}^p(\mathcal{K}_0, m, \mathcal{T}_h)$. Subsequently, on the local tensor meshes of the support of each of these hierarchical B-splines, we form a partition of unity by piecewise multilinear hat functions with merely $C^0(\Omega)$ continuity.

Let $\bar{p} \le p$ be a supplementary polynomial degree and let

$$ \bar{\mathcal{K}}_0 = (\bar{\mathcal{K}}_{1(0)}, \dots, \bar{\mathcal{K}}_{d(0)}) $$

be a fixed $d$-dimensional vector of $\bar{p}$-open knot vectors

$$ \bar{\mathcal{K}}_{i(0)} = (\bar{t}_{i(0),0}, \dots, \bar{t}_{i(0),\bar{N}_{i(0)}+\bar{p}}) $$
such that $(\bar{t}_{i(0),\bar{p}}, \dots, \bar{t}_{i(0),\overline{N}_{i(0)}})$ is a subsequence of $(t_{i(0),p}, \dots, t_{i(0),N_{i(0)}})$ which is obtained by reducing multiplicities of the latter knots to at least one and

$$0 = \bar{t}_{i(0),0} = \dots = \bar{t}_{i(0),\bar{p}} \quad \text{and} \quad \bar{t}_{i(0),\overline{N}_{i(0)}} = \dots = \bar{t}_{i(0),\overline{N}_{i(h)}+\bar{p}} = 1.$$

In particular, the tensor mesh corresponding to $\overline{\mathcal{K}}_0$ coincides with the initial tensor-mesh $\hat{T}_0$ corresponding to $\mathcal{K}_0$. To guarantee that

$$\mathbb{S}^{\bar{p}}(\overline{\mathcal{K}}_{i(0)}) \subseteq \mathbb{S}^{\mathcal{p}}(\mathcal{K}_{i(0)}), \qquad (29)$$

we further suppose that

$$p - \#t_{i(0),j} \le \bar{p} - \overline{\#}t_{i(0),j} \qquad (30)$$

for all interior knots $t_{i(0),j}$ in $(0,1)$ (which determines the smoothness of the considered splines), where $\overline{\#}$ denotes the multiplicity within $\overline{\mathcal{K}}_{i(0)}$. Next, we set $\overline{\mathcal{K}}_{\text{uni}(0)} := \overline{\mathcal{K}}_0$ and recursively define $\overline{\mathcal{K}}_{\text{uni}(\ell+1)}$ for $\ell \in \mathbb{N}_0$ as the uniform $h$-refinement of $\overline{\mathcal{K}}_{\text{uni}(\ell)}$ with fixed multiplicity $\overline{m} \in \mathbb{N}$ such that

$$p - m \le \bar{p} - \overline{m}. \qquad (31)$$

In words, the knots of $\overline{\mathcal{K}}_{\text{uni}(\ell+1)}$ in $[0,1]$ are the knots of $\overline{\mathcal{K}}_{\text{uni}(\ell)}$ in $[0,1]$ plus the points $(t_{i(\text{uni}(\ell)),j-1} + t_{i(\text{uni}(\ell)),j})/2$ with multiplicity $\overline{m}$ if $t_{i(\text{uni}(\ell)),j-1} < t_{i(\text{uni}(\ell)),j}$ for $j \in \{\overline{p} + 1, \dots, \overline{N}_{i(\text{uni}(\ell))}\}$. If, for example, all interior knots in $\mathcal{K}_0$ have the same multiplicity 1 and $m=1$, i.e., the corresponding splines are $C^{p-1}$ along initial as well as new lines, one can only choose $\overline{p} = p$ and $\overline{m} = 1$. If all interior knots in $\mathcal{K}_0$ have the same multiplicity $p$ and $m=p$, i.e., the corresponding splines are only $C^0$ along initial as well as new lines, one can choose $\overline{p} \le p$ and $\overline{m} \le \overline{p}$ arbitrarily, which leads us to $\overline{p} = 1$ and $\overline{m} = 1$. While the analysis below does not rely on this, we will always choose the smallest-possible polynomial degree $\overline{p}$ together with $\overline{m} := 1$. This is the most sensible choice from a practical as well as theoretical point of view, as the efficiency constant of our estimator depends on the degree $\overline{p}$; see Proposition 6.9 and Remark 6.10 as well as Proposition 8.6 and Remark 8.7 below. In Remark 8.8, we further discuss that the choice of $\overline{p}$ and $\overline{m}$ has little influence on the data oscillation terms arising in our error estimator.

Let again $\hat{T}_h$ be a hierarchical mesh with corresponding sets $(\hat{\Omega}_h^\ell)_{\ell \in \mathbb{N}_0}$ as in Section 4.2. Setting $\overline{p} := (\overline{p}, \dots, \overline{p})$, we define the B-splines $\overline{B}_{(\text{uni}(\ell),j,\overline{p})}$, the hierarchical basis $\{\overline{B}_a : a \in \overline{\mathcal{I}}_h\}$ with index set $\overline{\mathcal{I}}_h$, the level level($\overline{B}_a$) for $a \in \mathcal{I}_h$, and the spanned space of hierarchical splines $\mathbb{S}^p(\overline{\mathcal{K}}_0, \overline{m}, \hat{T}_h)$ as in (22)–(24). Again, [SM16, Section 3] gives an explicit characterization for the spanned space of hierarchical splines. Our assumptions on the knot multiplicities, which imply the nestedness $\mathbb{S}^{\overline{p}}(\overline{\mathcal{K}}_{\text{uni}(\ell)}) \subseteq \mathbb{S}^{\mathcal{p}}(\mathcal{K}_{\text{uni}(\ell)})$ for all $\ell \in \mathbb{N}_0$, thus give that

$$\mathbb{S}^{\overline{p}}(\overline{\mathcal{K}}_0, \overline{m}, \hat{T}_h) \subseteq \mathbb{S}^{\mathcal{p}}(\mathcal{K}_0, m, \hat{T}_h). \qquad (32)$$

Now $1 \in \mathbb{S}^{\overline{p}}(\overline{\mathcal{K}}_0, \overline{m}, \hat{T}_h)$ yields the existence of a partition of unity on the parameter domain

$$1 = \sum_{a \in \overline{\mathcal{I}}_h} c_a \overline{B}_a.$$

One can prove that the coefficients $c_a \in \mathbb{R}$ satisfy that

$$0 \le c_a \le 1 \qquad (33)$$
for all $\mathbf{a} \in \bar{\mathcal{I}}_h$; see, e.g., [BG16a, Lemma 3.2]. Consequently, we can define

$$ \psi_a := (c_a \bar{B}_a) \circ \mathbf{F}^{-1} \quad \text{for all } \mathbf{a} \in \mathcal{V}_h := \{\mathbf{a} \in \bar{\mathcal{I}}_h : c_a > 0\} \quad (34) $$

and observe that the $\psi_a$ form a partition of unity on the physical domain

$$ \sum_{a \in \mathcal{V}_h} \psi_a = 1 \quad \text{in } \Omega. $$

Henceforth, we call $\mathcal{V}_h$ the set of nodes and $\mathbf{a} \in \mathcal{V}_h$ a node. For further use, we abbreviate $\omega_a := \text{int}(\text{supp}(\psi_a))$ as well as $\hat{\omega}_a := \mathbf{F}^{-1}(\omega_a)$ for all $\mathbf{a} \in \mathcal{V}_h$; we will use the terminology large patches for $\omega_a$ or $\hat{\omega}_a$. Figure 6 gives an illustrative example.

Below, we will also crucially use a second partition of unity on each large patch $\omega_a$. Let $\hat{\mathcal{T}}_a$ be the smallest uniform tensor mesh refinement of $\{\hat{T} \in \hat{\mathcal{T}}_h : \hat{T} \subseteq \text{supp}(B_a)\}$, i.e.,

$$ \hat{\mathcal{T}}_a := \{\hat{T} \in \hat{\mathcal{T}}_{\text{uni}(\ell_a)} : \hat{T} \subseteq \text{supp}(B_a)\} $$

with $\ell_a := \max\{\text{level}(\hat{T}) : \hat{T} \in \hat{\mathcal{T}}_h \wedge \hat{T} \subseteq \text{supp}(B_a)\}$, and $\mathcal{T}_a$ the corresponding mesh of the large patch $\omega_a$ in the physical domain; see again Figure 6. Moreover, let $\mathbb{Q}^1(\hat{\mathcal{T}}_a)$ be the set of all $\hat{\mathcal{T}}_a$-piecewise tensor-product polynomials of degree $1 := (1, \dots, 1)$ and

$$ \mathbb{Q}^1(\mathcal{T}_a) := \{\hat{v} \circ \mathbf{F}^{-1}|_{\omega_a} : \hat{v} \in \mathbb{Q}^1(\hat{\mathcal{T}}_a)\}. \quad (35) $$

Finally, let $\mathcal{V}_h^a$ be the set of all vertices in the local mesh $\mathcal{T}_a$. We denote by $\psi_b$ the hat function associated with the vertex $\mathbf{b} \in \mathcal{V}_h^a$; this is the unique function in $\mathbb{Q}^1(\mathcal{T}_a) \cap C^0(\omega_a)$ taking value 1 in the vertex $\mathbf{b}$ and 0 in all other vertices from $\mathcal{V}_h^a$. Observe that the $\psi_b$ form a partition of unity on the large patches $\omega_a$

$$ \sum_{b \in \mathcal{V}_h^a} \psi_b = 1 \quad \text{in } \omega_a. $$

We abbreviate $\omega_b := \text{int}(\text{supp}(\psi_b))$ as well as $\hat{\omega}_b := \mathbf{F}^{-1}(\omega_b)$ for all $\mathbf{b} \in \mathcal{V}_h^a$, for which we use the name small patches. Figure 6 gives again an illustrative example.

### 5.2. Patchwise Sobolev spaces.
For a node $\mathbf{a} \in \mathcal{V}_h$, define a local Sobolev space on the large patch $\omega_a$ as

$$ H_*^1(\omega_a) := \begin{cases} \{v \in H^1(\omega_a) : (v, 1)_{\omega_a} = 0\} & \text{if } \psi_a \in H_0^1(\Omega), \\ \{v \in H^1(\omega_a) : v = 0 \text{ on } \partial\omega_a \setminus \psi_a^{-1}(\{0\})\} & \text{else.} \end{cases} \quad (36) $$

This is the mean-value-free subspace of $H^1(\omega_a)$ in the interior of $\Omega$, and the trace-free (on that part of $\partial\omega_a$ where $\psi_a$ is nonzero) subspace of $H^1(\omega_a)$ adjacent to the boundary of $\Omega$. For vector-valued functions, we will use

$$ H_0(\text{div}, \omega_a) := \begin{cases} \{v \in \mathbf{H}(\text{div}, \omega_a) : v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a\} & \text{if } \psi_a \in H_0^1(\Omega), \\ \{v \in \mathbf{H}(\text{div}, \omega_a) : v \cdot n_{\omega_a} = 0 \text{ on } \partial\omega_a \cap (\psi_a)^{-1}(\{0\})\} & \text{else,} \end{cases} \quad (37) $$

where $n_{\omega_a}$ denotes the outer normal vector on $\partial\omega_a$ and $v \cdot n_{\omega_a}$ is understood in the appropriate weak sense. These are the normal-component-free subspaces of $\mathbf{H}(\text{div}, \omega_a)$: everywhere on $\partial\omega_a$ in the interior of $\Omega$ and on that part of $\partial\omega_a$ where $\psi_a$ is zero adjacent to the boundary of $\Omega$.
For a node $\mathbf{a} \in \mathcal{V}_h$ and a vertex $\mathbf{b} \in \mathcal{V}_h^a$, we also define some spaces on the small patches $\omega_b$. In particular, we let

$$L_{*}^{2}(\omega_{b}) := \begin{cases} \{v \in L^{2}(\omega_{b}) : (v, 1)_{\omega_{b}} = 0\} & \text{if } \psi_{a}\psi_{b} \in H_{0}^{1}(\Omega), \\ L^{2}(\omega_{b}) & \text{else,} \end{cases} \qquad (38)$$

where $\psi_b$, being defined as function on $\omega_a$, is identified with its extension by zero onto $\Omega$ (which is in general not in $H^1(\Omega)$, cf. Figure 6 (c)). We also define, as in (36) and (37),

$$H_{*}^{1}(\omega_{b}) := \begin{cases} \{v \in H^{1}(\omega_{b}) : (v, 1)_{\omega_{b}} = 0\} & \text{if } \psi_{a}\psi_{b} \in H_{0}^{1}(\Omega), \\ \{v \in H^{1}(\omega_{b}) : v = 0 \text{ on } \partial\omega_{b} \setminus (\psi_{a}\psi_{b})^{-1}(\{0\})\} & \text{else} \end{cases} \qquad (39)$$

and

$$\mathbf{H}_0(\mathrm{div}, \omega_b) := \begin{cases} \{\mathbf{v} \in \mathbf{H}(\mathrm{div}, \omega_b) : \mathbf{v} \cdot \mathbf{n}_{\omega_b} = 0 \text{ on } \partial\omega_b\} & \text{if } \psi_a\psi_b \in H_0^1(\Omega), \\ \{\mathbf{v} \in \mathbf{H}(\mathrm{div}, \omega_b) : \mathbf{v} \cdot \mathbf{n}_{\omega_b} = 0 \text{ on } \partial\omega_b \cap (\psi_a\psi_b)^{-1}(\{0\})\} & \text{else.} \end{cases} \qquad (40)$$

Note that these spaces actually depend on both $\mathbf{a}$ and $\mathbf{b}$, where $\mathbf{a}$ is omitted in our notation.

Finally, for $\mathbf{a} \in \mathcal{V}_h$, $\mathbf{b} \in \mathcal{V}_h^a$, and $\omega \in \{\omega_a, \omega_b\}$, define the Poincaré–Friedrichs constant as the minimal constant $C_{\mathrm{PF}}(\omega) > 0$ such that

$$\|v\|_{\omega} \leq \mathrm{diam}(\omega) C_{\mathrm{PF}}(\omega) \|\nabla v\|_{\omega} \quad \text{for all } v \in H_{*}^{1}(\omega). \qquad (41)$$

Note that $C_{\mathrm{PF}}(\omega)$ only depends on the shape of $\omega$ and $\partial\omega$ where respectively $\psi_a$ or $\psi_a\psi_b$ is nonzero; if $\psi_a \in H_0^1(\Omega)$ or $\psi_a\psi_b \in H_0^1(\Omega)$ and for convex $\omega$, there in particular holds $C_{\mathrm{PF}}(\omega) \leq 1/\pi$ (“interior” cases). In the other, “boundary”, cases, $C_{\mathrm{PF}}(\omega) \leq 1$ when there exists a unit vector $\mathbf{m}$ such that the straight semi-line of direction $\mathbf{m}$ originating at (almost) any point in $\omega$ hits the boundary $\partial\omega$ there where respectively $\psi_a$ or $\psi_a\psi_b$ is nonzero, cf., e.g., [Voh05, VV12] and the references therein.

### 5.3. Patchwise discrete subspaces.
For a node $\mathbf{a} \in \mathcal{V}_h$, let us define the $H_*^1(\omega_a)$-conforming subspace of the mapped piecewise multilinear functions $\mathbb{Q}^1(\mathcal{T}_a)$ from (35) as

$$V_h^a := \mathbb{Q}^1(\mathcal{T}_a) \cap H_*^1(\omega_a). \qquad (42)$$

Define the vector-valued contravariant Piola transform

$$\mathbf{\Phi}(\cdot) := (\det(DF)^{-1}(DF)(\cdot)) \circ \mathbf{F}^{-1} \qquad (43)$$

and the scalar Piola transform by

$$\tilde{\mathbf{\Phi}}(\cdot) := (\det(DF)^{-1}(\cdot)) \circ \mathbf{F}^{-1}, \qquad (44)$$

which satisfy the identity

$$\tilde{\mathbf{\Phi}}(\nabla \cdot (\cdot)) = \nabla \cdot \mathbf{\Phi}(\cdot), \qquad (45)$$

see, e.g., [EG21, Chapter 9].

For a node $\mathbf{a} \in \mathcal{V}_h$ and a vertex $\mathbf{b} \in \mathcal{V}_h^a$, define the meshes of the small patches as $\mathcal{T}_b := \{T \in \mathcal{T}_a : T \subseteq \bar{\omega}_b\}$, $\tilde{\mathcal{T}}_b := \{\mathbf{F}^{-1}(T) : T \in \mathcal{T}_b\}$; note that the elements $\tilde{T} \in \tilde{\mathcal{T}}_b$ are rectangles for $d=2$ and rectangular cuboids for $d=3$. Let $\mathbb{Q}^{\tilde{p}}(\tilde{\mathcal{T}}_b)$ be the space of all $\tilde{\mathcal{T}}_b$-piecewise polynomials of some fixed degree $\tilde{p} = (\tilde{p}_1, \dots, \tilde{p}_d)$ in each coordinate, $\tilde{p} \geq 0$,
and let $L^2_*(\hat{\omega}_b)$ be defined as in (38) with $\omega_b$ replaced by $\hat{\omega}_b$. We then define the local spaces

$$ \hat{Q}_h^{a,b} := \mathbb{Q}^{\tilde{p}}(\hat{\mathcal{T}}_b) \cap L^2_*(\hat{\omega}_b), \quad (46a) $$

$$ Q_h^{a,b} := \{\tilde{\Phi}(\hat{q}_h) : \hat{q}_h \in \hat{Q}_h^{a,b}\} = \tilde{\Phi}(\hat{Q}_h^{a,b}). \quad (46b) $$

Note that since

$$ \tilde{\Phi}(L^2_*(\hat{\omega}_b)) = L^2_*(\omega_b), $$

$Q_h^{a,b}$ is contained in $L^2_*(\omega_b)$, i.e.,

$$ Q_h^{a,b} \subset L^2_*(\omega_b). \quad (47) $$

The mean-value constraint in (38) (when $\psi_a\psi_b \in H_0^1(\Omega)$) makes $Q_h^{a,b}$ a constrained subspace of mapped piecewise polynomials from $\mathbb{Q}^{\tilde{p}}(\hat{\mathcal{T}}_b)$ scaled by the factor $\det(DF)^{-1}$. With the set $\hat{Q}_h := \mathbb{Q}^{\tilde{p}}(\hat{\mathcal{T}}_h)$ of all $\hat{\mathcal{T}}_h$-piecewise polynomials of degree $\tilde{p}$, we also define the global (unconstrained) space via the mapping $\tilde{\Phi}$ as

$$ Q_h := \{\tilde{\Phi}(\hat{q}_h) : \hat{q}_h \in \hat{Q}_h\} = \tilde{\Phi}(\hat{Q}_h). $$

Let

$$ \mathbf{R}\mathbf{T}^{\tilde{p}}(\hat{\mathcal{T}}_b) := \begin{cases} \mathbb{Q}^{\tilde{p}+(1,0)}(\hat{\mathcal{T}}_b) \times \mathbb{Q}^{\tilde{p}+(0,1)}(\hat{\mathcal{T}}_b) & \text{if } d=2, \\ \mathbb{Q}^{\tilde{p}+(1,0,0)}(\hat{\mathcal{T}}_b) \times \mathbb{Q}^{\tilde{p}+(0,1,0)}(\hat{\mathcal{T}}_b) \times \mathbb{Q}^{\tilde{p}+(0,0,1)}(\hat{\mathcal{T}}_b) & \text{if } d=3 \end{cases} $$

be the usual broken (elementwise) Raviart-Thomas space on the rectangular/rectangular cuboid mesh $\hat{\mathcal{T}}_b$, see, e.g., [BBF13, Section 2.4.1]. We then set

$$ \hat{V}_h^{a,b} := \mathbf{R}\mathbf{T}^{\tilde{p}}(\hat{\mathcal{T}}_b) \cap \mathbf{H}_0(\mathrm{div}, \hat{\omega}_b), \quad (48a) $$

$$ V_h^{a,b} := \{\Phi(\hat{v}_h) : \hat{v}_h \in \hat{V}_h^{a,b}\} = \Phi(\hat{V}_h^{a,b}). \quad (48b) $$

Since

$$ \mathbf{H}_0(\mathrm{div}, \omega_b) = \Phi(\mathbf{H}_0(\mathrm{div}, \hat{\omega}_b)), $$

we also have

$$ V_h^{a,b} = \mathbf{R}\mathbf{T}^{\tilde{p}}(\mathcal{T}_b) \cap \mathbf{H}_0(\mathrm{div}, \omega_b), \quad (49) $$

where $\mathbf{R}\mathbf{T}^{\tilde{p}}(\mathcal{T}_b)$ is the space $\mathbf{R}\mathbf{T}^{\tilde{p}}(\hat{\mathcal{T}}_b)$ mapped by the Piola transform

$$ \mathbf{R}\mathbf{T}^{\tilde{p}}(\mathcal{T}_b) := \{\Phi(\hat{v}_h) : \hat{v}_h \in \mathbf{R}\mathbf{T}^{\tilde{p}}(\hat{\mathcal{T}}_b)\} = \Phi(\mathbf{R}\mathbf{T}^{\tilde{p}}(\hat{\mathcal{T}}_b)). \quad (50) $$

Crucially, by construction, see [BBF13, Section 2.4.1],

$$ \nabla \cdot \hat{V}_h^{a,b} = \hat{Q}_h^{a,b}, \quad (51a) $$

whereas by the Piola transform identity (45) and definition (46b), one also has

$$ \nabla \cdot V_h^{a,b} = Q_h^{a,b}. \quad (51b) $$

**Remark 5.1.** While it is in principle not relevant for the quality of our a posteriori estimator whether the knot lines corresponding to $V_h$ are aligned with those of a NURBS parametrization $\mathbf{F}$ (see Remark 4.2), this is essential for $Q_h^{a,b}$ and $\mathbf{R}\mathbf{T}^{\tilde{p}}(\mathcal{T}_b)$. Indeed, the latter two spaces must exhibit good approximation properties to obtain small oscillation terms (77b) and (84c). On the other hand, as $Q_h^{a,b}$ and $\mathbf{R}\mathbf{T}^{\tilde{p}}(\mathcal{T}_b)$ are discontinuous (piecewise with respect to the mesh $\mathcal{T}_b$), they do not see/need the regularity of $\mathbf{F}$ over the knot lines. We also refer to Remark 8.3 and Remark 8.5 for conditions under which the oscillations even vanish.
## 6. ABSTRACT ASSUMPTIONS

In this section, we attempt to describe as clearly as possible the underlying principles of a posteriori error analysis by equilibrated fluxes in the present context. For this purpose, we identify six abstract assumptions under which our subsequent analysis can be carried out. We then immediately verify these assumptions in the particular context of Sections 2–5.

As for the general setting, we merely need to assume:

**Assumption 6.1.** $\hat{\Omega}$ is an open bounded connected Lipschitz domain in $\mathbb{R}^d$, $\mathbf{F}$ is a bi-Lipschitz mapping, $\Omega := \mathbf{F}(\hat{\Omega})$, and $V_h$ is an arbitrary subspace of $H_0^1(\Omega)$.

As for the partitions of unity, the minimalist assumptions are (note that we do not require that the partitions are non-negative, i.e., $\psi_a, \psi_b \ge 0$):

**Assumption 6.2.** There is a finite index set $\mathcal{V}_h$ and functions $\psi_a$ such that

$$ \{\psi_a : \mathbf{a} \in \mathcal{V}_h\} \subset W^{1,\infty}(\Omega) \quad (52a) $$

form a partition of unity over $\Omega$ in the sense that

$$ \sum_{\mathbf{a} \in \mathcal{V}_h} \psi_{\mathbf{a}} = 1 \quad \text{in } \Omega. \quad (52b) $$

The interior $\omega_a$ of the support of any $\psi_a$ is a connected Lipschitz domain with $|\omega_a| > 0$. Moreover,

$$ \{\psi_a : \mathbf{a} \in \mathcal{V}_h\} \cap H_0^1(\Omega) \subset V_h. \quad (52c) $$

**Assumption 6.3.** For any node $\mathbf{a} \in \mathcal{V}_h$, there is a finite set of vertices $\mathcal{V}_h^a$ and functions $\psi_b$ such that

$$ \{\psi_b : \mathbf{b} \in \mathcal{V}_h^a\} \subset W^{1,\infty}(\omega_a) \quad (53a) $$

form a partition of unity over $\omega_a$ in that

$$ \sum_{\mathbf{b} \in \mathcal{V}_h^a} \psi_{\mathbf{b}} = 1 \quad \text{in } \omega_a. \quad (53b) $$

The interior $\omega_b$ of the support of any $\psi_b$ is a connected Lipschitz domain with $|\omega_b| > 0$. Moreover, for $H_*^1(\omega_a)$ given by (36), all $\psi_b$ such that $\psi_a\psi_b \in H_0^1(\Omega)$, where $\psi_b$ is identified with its extension by zero onto $\Omega$, are contained up to additive constants in a finite-dimensional subspace $V_h^a \subset H_*^1(\omega_a)$, i.e.,

$$ \{\psi_b : \mathbf{b} \in \mathcal{V}_h^a \wedge \psi_a \psi_b \in H_0^1(\Omega)\} \subset \begin{cases} \{v_h + c : v_h \in V_h^a, c \in \mathbb{R}\} & \text{if } \psi_a \in H_0^1(\Omega) \\ V_h^a & \text{else.} \end{cases} \quad (53c) $$

**Remark 6.4.** While in theory, there is almost no connection between the functions $\psi_a$ and the functions $\psi_b$, in practice, $\psi_a$ are mapped piecewise polynomials with high smoothness on some local mesh and the corresponding $\psi_b$ are mapped piecewise polynomials with low smoothness on the same mesh, see Sections 3.3 and 5.1. We stress that the choice of these partitions of unity determines the quality of the equilibrated flux estimator through the oscillations terms in Propositions 8.2 and 8.6.

Next, we recall the space $H_*^1(\omega_b)$ from (39), the Poincaré–Friedrichs inequality (41), and assume the following set of estimates:
**Assumption 6.5.** There exist generic positive constants $C_1, \dots, C_6 > 0$ such that for all $\mathbf{a} \in \mathcal{V}_h$ and $\mathbf{b} \in \mathcal{V}_h^a$, it holds that

$$ \|\psi_{\mathbf{a}}\|_{\infty, \omega_{\mathbf{a}}} \leq C_1, \qquad (54a) $$

$$ \|\nabla \psi_{\mathbf{a}}\|_{\infty, \omega_{\mathbf{a}}} \mathrm{diam}(\omega_{\mathbf{a}}) C_{\mathrm{PF}}(\omega_{\mathbf{a}}) \leq C_2, \qquad (54b) $$

$$ \|\nabla \psi_{\mathbf{a}}\|_{\infty, \omega_{\mathbf{b}}} \mathrm{diam}(\omega_{\mathbf{b}}) C_{\mathrm{PF}}(\omega_{\mathbf{b}}) \leq C_3, \qquad (54c) $$

$$ \sup_{\mathbf{x} \in \omega_{\mathbf{a}}} \#\{\mathbf{b}' \in \mathcal{V}_h^a : \mathbf{x} \in \omega_{\mathbf{b}'}\} \leq C_4, \qquad (54d) $$

$$ \|\psi_{\mathbf{b}}\|_{\infty, \omega_{\mathbf{b}}} \leq C_5, \qquad (54e) $$

$$ \|\nabla \psi_{\mathbf{b}}\|_{\infty, \omega_{\mathbf{b}}} \mathrm{diam}(\omega_{\mathbf{b}}) C_{\mathrm{PF}}(\omega_{\mathbf{b}}) \leq C_6. \qquad (54f) $$

Recall the space $\mathbf{H}_0(\mathrm{div}, \omega_{\mathbf{b}})$ from (40), the subspace $L_*^2(\omega_{\mathbf{b}})$ of $L^2(\omega_{\mathbf{b}})$ containing functions with mean value zero if $\psi_{\mathbf{a}}\psi_{\mathbf{b}} \in H_0^1(\Omega)$ from (38), and the similar space $L_*^2(\hat{\omega}_{\mathbf{b}})$ in the parameter domain. Recall also the contravariant Piola transform $\mathbf{\Phi}$ from (43). For local flux equilibration, we will rely on discrete spaces $Q_h^{a,b}$ and $V_h^{a,b}$ satisfying the following:

**Assumption 6.6.** For any node $\mathbf{a} \in \mathcal{V}_h$ and any vertex $\mathbf{b} \in \mathcal{V}_h^a$, there are finite-dimensional subspaces

$$ Q_h^{a,b} \subset L_*^2(\omega_{\mathbf{b}}) \quad \text{and} \quad V_h^{a,b} \subset \mathbf{H}_0(\mathrm{div}, \omega_{\mathbf{b}}) \qquad (55a) $$

satisfying the compatibility condition

$$ \nabla \cdot \mathbf{V}_h^{a,b} = Q_h^{a,b}. \qquad (55b) $$

Moreover, we suppose the existence of a global space $Q_h \subset L^2(\Omega)$ such that

$$ Q_h|_{\omega_{\mathbf{b}}} \subseteq Q_{h,c}^{a,b} := \begin{cases} \{q_h + \tilde{\Phi}(c) : q_h \in Q_h^{a,b}, c \in \mathbb{R}\} & \text{if } \psi_{\mathbf{a}}\psi_{\mathbf{b}} \in H_0^1(\Omega), \\ Q_h^{a,b} & \text{else.} \end{cases} \qquad (55c) $$

**Remark 6.7.** Using the Piola transforms $\mathbf{\Phi}$ and $\tilde{\Phi}$ from (43)–(44), let

$$ \hat{Q}_h^{a,b} := \tilde{\Phi}^{-1}(Q_h^{a,b}) \subset L_*^2(\hat{\omega}_{\mathbf{b}}) = \tilde{\Phi}^{-1}(L_*^2(\omega_{\mathbf{b}})) \qquad (56) $$

and

$$ \hat{\mathbf{V}}_h^{a,b} := \mathbf{\Phi}^{-1}(\mathbf{V}_h^{a,b}) \subset \mathbf{H}_0(\mathrm{div}, \hat{\omega}_{\mathbf{b}}) = \mathbf{\Phi}^{-1}(\mathbf{H}_0(\mathrm{div}, \omega_{\mathbf{b}})). \qquad (57) $$

From (55b) and (45), we in particular have

$$ \nabla \cdot \hat{\mathbf{V}}_h^{a,b} = \hat{Q}_h^{a,b}, $$

as in Section 5.3. The role of the global space $Q_h \subset L^2(\Omega)$ will be prominent below: please note that it is related to neither the node $\mathbf{a}$, nor to the vertex $\mathbf{b}$; this forces the local spaces $Q_{h,c}^{a,b}$ to contain patch-independent “base-blocks”. In practice, $Q_h|_{\omega_{\mathbf{b}}} = Q_{h,c}^{a,b}$ only for uniform mesh refinement but not, for example, in the setting of Figure 6, where a strict inclusion holds, since the local spaces $Q_{h,c}^{a,b}$ stem from the local uniform fine meshes $\mathcal{T}_a$ (in grey in Figure 6), whereas $Q_h$ is related to the (mapped) hierarchical mesh $\mathcal{T}_h$ (in black in Figure 6).

Finally, we will essentially employ the following finite dimension to infinite dimension extension property.
**Assumption 6.8.** There exist a generic constant $C_{st} \ge 1$, as well as a superspace $RT_h^{a,b} \subset [L^2(\omega_b)]^d$ verifying $V_h^{a,b} = RT_h^{a,b} \cap H_0(\text{div}, \omega_b)$ for all nodes $a \in V_h$ and vertices $b \in V_h^a$, such that

$$ \min_{\substack{v_h \in V_h^{a,b} \\ \nabla \cdot v_h = g_h}} \|v_h + \tau_h\|_{\omega_b} \le C_{st} \min_{\substack{v \in H_0(\text{div}, \omega_b) \\ \nabla \cdot v = g_h}} \|v + \tau_h\|_{\omega_b} \quad (58) $$

for all $g_h \in \nabla \cdot V_h^{a,b} = Q_h^{a,b}$ and all $\tau_h \in RT_h^{a,b}$.

We now verify that the above assumptions are satisfied in our IGA context:

**Proposition 6.9.** Assumptions 6.1–6.8 are satisfied in the context of Sections 2–5 with the choice $RT_h^{a,b} = RT_{\tilde{\mathcal{T}}_b}^{\tilde{p}}$. In particular, the constants in Assumption 6.5 can be taken as $C_1 = 1$, $C_4 = 2^d$, $C_5 = 1$ and such that $C_2, C_3, C_6$, and $C_{st}$ only depend on the shapes in $\tilde{\mathcal{T}}_{\text{uni}(0)}$ as well as the mapping $\mathbf{F}$ via $\max\{\|DF\|_{\infty,\hat{\Omega}}, \|(DF)^{-1}\|_{\infty,\hat{\Omega}}\}; C_2$ and $C_3$ additionally depend on the supplementary polynomial degree $\bar{p}$ from (30)–(31).

*Proof.* Assumptions 6.1–6.3 are immediately satisfied with the choices made in Sections 2–5, in particular fixing the space $V_h^a$ following (42).

We now verify Assumption 6.5. Thanks to (18) and (33), inequality (54a) is satisfied with the constant $C_1 = 1$, whereas (54d) and (54e) hold easily with respectively $C_4 = 2^d$ and $C_5 = 1$. Next, (54b), (54c), and (54f) follow from the facts that $\|\nabla \psi_a\|_{\infty, \omega_a} \lesssim \bar{p} \text{diam}(\omega_a)^{-1}$, $\text{diam}(\omega_b) \le \text{diam}(\omega_a)$, and $\|\nabla \psi_b\|_{\infty, \omega_b} \lesssim \text{diam}(\omega_b)^{-1}$, see, e.g., [BdVBSV14, Equation (2.7)], where we also use that $C_{\text{PF}}(\omega_a), C_{\text{PF}}(\omega_b) \lesssim 1$ (see, e.g., [VV12, Corollary 2.2]). Here, $A \lesssim B$ means that $A \le CB$ for a hidden constant $C > 0$ depending only on the shapes of the elements in $\tilde{\mathcal{T}}_{\text{uni}(0)}$ as well as on $\max\{\|DF\|_{\infty, \hat{\Omega}}, \|(DF)^{-1}\|_{\infty, \hat{\Omega}}\}$.

We now turn to Assumption 6.6. Properties (55a) and (55c) are trivially satisfied for the spaces defined in Section 5.3, see in particular (47) and (49), whereas the compatibility condition (55b) follows from the construction requirement of the Raviart–Thomas space in the parameter domain (51a) by (45) and definition (46b), see (51b).

We finally address Assumption 6.8. For the spaces defined in Section 5.3, where we take $RT_h^{a,b} = RT_{\tilde{\mathcal{T}}_b}^{\tilde{p}}$, this result is proved in the parameter domain in [BPS09, Theorems 5 and 7] in two space dimensions; in three space dimensions, one needs to rely instead on [CDD08], cf. [EV20, Theorem 2.5 and Corollary 3.3] building on [CM10, DGS12]. We give a proof in the physical domain in Appendix A, where the resulting constant $C_{st}$ depends only on the shapes of the elements in $\tilde{\mathcal{T}}_{\text{uni}(0)}$ and $\max\{\|DF\|_{\infty, \hat{\Omega}}, \|(DF)^{-1}\|_{\infty, \hat{\Omega}}\}$. We stress that $C_{st}$ is independent of the polynomial degree $\tilde{p}$. $\square$

**Remark 6.10.** Recall that the supplementary polynomial degree $\bar{p}$ from Section 5.1 depends on the considered smoothness but not necessarily on the polynomial degree $p$. In this sense, $C_2, C_3$ in Proposition 6.9 are independent of the polynomial degree $p$. We recall in particular that $\bar{p} = 1$ can be taken for $C^0$ splines, and in general $\bar{p} = k+1$ for $C^k$ splines. We admit, though, that $\bar{p} = p$ for $C^{p-1}$ splines, see the discussion in Section 5.1.

**Remark 6.11.** We mention that the abstract Assumptions 6.1–6.6 are also satisfied in the setting of Sections 3.1 and 3.2. Indeed, one can simply choose one function $\psi_b := 1$ on each patch $\omega_a$ so that $\omega_b = \omega_a$, $V_h^a := \{0\}$, $V_h^{a,b} := RT_{\tilde{\mathcal{T}}_a}^{\tilde{p}}(\mathcal{T}_a) \cap H_0(\text{div}, \omega_a)$ with $RT_{\tilde{\mathcal{T}}_a}^{\tilde{p}}(\mathcal{T}_a)$ from (7), $Q_h^{a,b} := \mathbb{Q}^{\tilde{p}}(\mathcal{T}_a) \cap L_*^2(\omega_a)$ and $Q_h := \mathbb{Q}^{\tilde{p}}(\mathcal{T}_h)$. For Assumption 6.8, in turn, to our knowledge, there only is a rigorous proof in the setting of Section 3.1 with “small” vertex patches $\omega_a$, and not in Section 3.2 with “large” vertex patches $\omega_a$.
## 7. INEXPENSIVE EQUILIBRATION

Let the abstract assumptions of Section 6 be satisfied, let $u$ solve (2), and let $u_h$ solve (3). The partition of unity functions $\psi_a$ from Assumption 6.2, in our setting the hierarchical B-splines $\psi_a$ of polynomial degree $\bar{p} = (\bar{p}, \dots, \bar{p})$ and multiplicity $\bar{m}$ from (34), lead to

$$ (\nabla u_h, \nabla \psi_a)_{\omega_a} = (f, \psi_a)_{\omega_a} \quad \text{for all } a \in \mathcal{V}_h \text{ with } \psi_a \in H_0^1(\Omega), \quad (59) $$

which is an immediate consequence of (3) and (52c). Thus, in view of (2),

$$ (\nabla(u - u_h), \nabla \psi_a)_{\omega_a} = 0 \quad \text{for all } a \in \mathcal{V}_h \text{ with } \psi_a \in H_0^1(\Omega). \quad (60) $$

This orthogonality is sufficient to localize the error $\|\nabla(u-u_h)\|_{\Omega}$ (or, equivalently, the dual norm of the residual) over the large patches $\omega_a$, see [BMV20] and the references therein, and then a flux equilibration can be easily devised. The issue, however, is that this leads to an expensive equilibration with higher-order (related to the polynomial degree $p$) mixed finite spaces on the large patches $\omega_a$, see Section 3.2. Our goal below is to design a much less expensive equilibration, with some inexpensive (typically piecewise multilinear) solve on the large patches $\omega_a$ followed by higher-order (related to $p$) mixed finite element solves on the small patches $\omega_b$ only, extending the construction from Section 3.3 to the present general setting. In order to achieve it, we crucially rely on the discrete patchwise spaces from Section 6/Section 5.3.

### 7.1. Inexpensive (lowest-order) residual lifting on the large patches $\omega_a$.
Our first step is to construct a discrete residual function $r_h^a \in V_h^a$, where $V_h^a$ is the local space from Assumption 6.3. Recall that in the context of hierarchical B-splines, $V_h^a$ merely consists of mapped piecewise multilinear functions, see (35) and (42).

**Definition 7.1.** For all nodes $a \in \mathcal{V}_h$, let $r_h^a \in V_h^a$ be such that

$$ (\nabla r_h^a, \nabla v_h)_{\omega_a} = (f, v_h \psi_a)_{\omega_a} - (\nabla u_h, \nabla (v_h \psi_a))_{\omega_a} \quad \text{for all } v_h \in V_h^a. \quad (61) $$

Figure 3, steps 1)–3), gives an illustration for (61) in the spline setting, which is an inexpensive scalar-valued lowest-order local problem on the larger patches $\omega_a$, lifting the $\psi_a$-weighted residual, cf. [CF00, BPS09, EV15, BMV20]. We note that (61) is a finite-dimensional version of the problem: find $r^a \in H_*^1(\omega_a)$ such that

$$ (\nabla r^a, \nabla v)_{\omega_a} = (f, v \psi_a)_{\omega_a} - (\nabla u_h, \nabla (v \psi_a))_{\omega_a} \quad \text{for all } v \in H_*^1(\omega_a). \quad (62) $$

This can be equivalently written as: find $r^a \in H_*^1(\omega_a)$ such that

$$ (\nabla r^a + \psi_a \nabla u_h, \nabla v)_{\omega_a} = (f \psi_a - \nabla u_h \cdot \nabla \psi_a, v)_{\omega_a} \quad \text{for all } v \in H_*^1(\omega_a). \quad (63) $$

From (63), we see that $-(\nabla r^a + \psi_a \nabla u_h)$ lies in $H_0(\text{div}, \omega_a)$ with the divergence equal to $f \psi_a - \nabla u_h \cdot \nabla \psi_a$. It can be also characterized by

$$ \begin{array}{ll} \underset{v \in H_0(\text{div}, \omega_a)}{\operatorname{argmin}} & \|v + \psi_a \nabla u_h\|_{\omega_a}, \quad (64) \\ \nabla \cdot v = f \psi_a - \nabla u_h \cdot \nabla \psi_a & \end{array} $$

which is the infinite-dimensional paradigm of (9).
7.2. **Projection operators for general bi-Lipschitz mappings F.** In order to proceed, we need to define some projection operators. Recall the local space $Q_h^{a,b}$ from Assumption 6.6 as well as $\hat{Q}_h^{a,b}$ from (56). For a given $g \in L^2(\omega_b)$, we define $\Upsilon_{Q_h^{a,b}g} \in Q_h^{a,b}$ as the Petrov-Galerkin projection of $g$ into $Q_h^{a,b}$ with the test space $\hat{Q}_h^{a,b} \circ F^{-1} = \{\hat{q}_h \circ F^{-1} : \hat{q}_h \in \hat{Q}_h^{a,b}\}$, i.e.,

$$ (\Upsilon_{Q_h^{a,b}g, q_h})_{\omega_b} = (g, q_h)_{\omega_b} \quad \text{for all } q_h \in \hat{Q}_h^{a,b} \circ F^{-1}. \qquad (65) $$

When the mapping $F$ is affine, the scaling factor $\det(DF)^{-1}$ in (44) is constant, so that the trial space $Q_h^{a,b}$ and the test space $\hat{Q}_h^{a,b} \circ F^{-1}$ become the same. In such a case, $\Upsilon_{Q_h^{a,b}}$ is simply the $L^2(\omega_b)$-orthogonal (Galerkin) projection. Let

$$ \begin{aligned} \hat{Q}_{h,c}^{a,b} &:= \{\tilde{\Phi}^{-1}(q_h) : q_h \in Q_{h,c}^{a,b}\} = \tilde{\Phi}^{-1}(Q_{h,c}^{a,b}), \\ \hat{Q}_h &:= \{\tilde{\Phi}^{-1}(q_h) : q_h \in Q_h\} = \tilde{\Phi}^{-1}(Q_h). \end{aligned} $$

For $g \in L^2(\omega_b)$ and $g \in L^2(\Omega)$, respectively, we define $\Upsilon_{Q_h^{a,b}g} \in Q_h^{a,b}$ and $\Upsilon_{Q_h}g \in Q_h$ analogously, i.e.,

$$ (\Upsilon_{Q_h^{a,b}g, q_h})_{\omega_b} = (g, q_h)_{\omega_b} \quad \text{for all } q_h \in \hat{Q}_{h,c}^{a,b} \circ F^{-1}, \qquad (66a) $$

$$ (\Upsilon_{Q_h}g, q_h)_{\Omega} = (g, q_h)_{\Omega} \quad \text{for all } q_h \in \hat{Q}_h \circ F^{-1}. \qquad (66b) $$

7.3. **Equilibrated flux on the small patches $\omega_b$.** Let $r_h^a$ be given by Definition 7.1 and recall the projection operator $\Upsilon_{Q_h^{a,b}}$ from Section 7.2. Then our equilibrated flux on the small patches $\omega_b$ is given by:

**Definition 7.2.** For all nodes $a \in \mathcal{V}_h$ and all vertices $b \in \mathcal{V}_h^a$, let

$$ \begin{aligned} \sigma_h^{a,b} &:= \underset{v_h \in \mathcal{V}_h^{a,b}}{\operatorname{argmin}} & \|v_h + \psi_b(\psi_a \nabla u_h + \nabla r_h^a)\|_{\omega_b}. \qquad (67) \\ & \nabla \cdot v_h = \Upsilon_{Q_h^{a,b}}(f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b) - \nabla r_h^a \cdot \nabla \psi_b) \end{aligned} $$

Figure 3, steps 4)–5), gives an illustration in the spline setting. From (63), we know that $-(\nabla r^a + \psi_a \nabla u_h)$ lies in $H_0(\text{div}, \omega_a)$ with the divergence equal to $f\psi_a - \nabla u_h \cdot \nabla \psi_a$. Thus, its cut-off by $\psi_b$, $-\psi_b(\psi_a \nabla u_h + \nabla r^a)$, lies in $H_0(\text{div}, \omega_b)$ with the divergence equal to $f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b) - \nabla r^a \cdot \nabla \psi_b$. Similarly to (64), it can be characterized implicitly by

$$ \begin{aligned} & \underset{v \in H_0(\text{div}, \omega_b)}{\operatorname{argmin}} & \|v + \psi_b(\psi_a \nabla u_h + \nabla r^a)\|_{\omega_b}. \qquad (68) \\ & \nabla \cdot v = f \psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b) - \nabla r^a \cdot \nabla \psi_b \end{aligned} $$

The flux $\sigma_h^{a,b}$ from (67) is then its discrete approximation.

The following lemma guarantees the existence and uniqueness of the minimizer of Definition 7.2:

**Lemma 7.3.** There exists a unique minimizer $\sigma_h^{a,b}$ of (67).

*Proof.* Let us abbreviate

$$ g := f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b) - \nabla r_h^a \cdot \nabla \psi_b \qquad (69) $$

and

$$ \tau := \psi_b(\psi_a \nabla u_h + \nabla r_h^a). $$

By definition of $\Upsilon_{Q_h^{a,b}}$, $\Upsilon_{Q_h^{a,b}g} \in Q_h^{a,b} = \nabla \cdot V_h^{a,b}$, so that the minimization set is nonempty, and existence and uniqueness follow by standard convexity arguments, see, e.g., [ET99,
Proposition 1.2]. In more details, problem (67) is equivalent to finding $\sigma_h^{a,b}$ with $\nabla \cdot \sigma_h^{a,b} = \Upsilon_{Q_h^{a,b}} g$ such that

$$ (\sigma_h^{a,b}, v_h)_{\omega_b} = -(\tau, v_h)_{\omega_b} \quad \text{for all } v_h \in V_h^{a,b} \text{ with } \nabla \cdot v_h = 0, \quad (70) $$

which is a square linear system. Existence and uniqueness thus follow when $\sigma_h^{a,b} = 0$ for zero data. Let thus $g = 0$ and $\tau = 0$. Since $\Upsilon_{Q_h^{a,b}} 0 = 0$ follows from (65), we can take $v_h = \sigma_h^{a,b}$ in (70), which implies $\|\sigma_h^{a,b}\|_{\omega_b} = 0$ and thus $\sigma_h^{a,b} = 0$. $\square$

Recalling (65) and (66a), the following is an important extension of Lemma 7.3:

**Lemma 7.4.** The projection $\Upsilon_{Q_h^{a,b}}$ in (67) can be replaced by $\Upsilon_{Q_{h,c}^{a,b}}$, i.e.,

$$ (\Upsilon_{Q_h^{a,b}} - \Upsilon_{Q_{h,c}^{a,b}})(f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b) - \nabla r_h^a \cdot \nabla \psi_b) = 0. \quad (71) $$

*Proof.* From (55c), the spaces $Q_h^{a,b}$ and $Q_{h,c}^{a,b}$ only differ if $\psi_a\psi_b \in H_0^1(\Omega)$, the difference only being transformed constants that remove the constraint. Let thus $\psi_a\psi_b \in H_0^1(\Omega)$ and recall the notation (69). From $V_h^{a,b} \subset H_0(\text{div}, \omega_b)$ (recall (40)) and (67), $0 = (\nabla \cdot \sigma_h^{a,b}, 1)_{\omega_b} = (\Upsilon_{Q_h^{a,b}} g, 1)_{\omega_b}$. Taking into account (65) and (66a), we thus only need to verify the Neumann compatibility condition $(g, 1)_{\omega_b} = 0$. Hence, test $g$ with 1 and obtain that

$$ (g, 1)_{\omega_b} = (f\psi_a, \psi_b)_{\omega_b} - (\nabla u_h \cdot \nabla \psi_a, \psi_b)_{\omega_b} - (\psi_a \nabla u_h, \nabla \psi_b)_{\omega_b} - (\nabla r_h^a, \nabla \psi_b)_{\omega_b}. \quad (72) $$

We consider two cases.

First, let $\psi_a \in H_0^1(\Omega)$. Then, we use that $\text{supp}(\psi_b) \subseteq \text{supp}(\psi_a)$ from (53a), (61) with the test function $v_h = \psi_b - |\omega_a|^{-1}(\psi_b, 1)_{\omega_a}$ having zero mean value on $\omega_a$ (which is possible due to (53c)), and (59) to see for the last term in (72) that

$$ \begin{align*} (\nabla r_h^a, \nabla \psi_b)_{\omega_b} &= (\nabla r_h^a, \nabla (\psi_b - |\omega_a|^{-1}(\psi_b, 1)_{\omega_a}))_{\omega_a} \\ &\stackrel{(61)}{=} (f\psi_a - \nabla u_h \cdot \nabla \psi_a, \psi_b - |\omega_a|^{-1}(\psi_b, 1)_{\omega_a})_{\omega_a} - (\psi_a \nabla u_h, \nabla \psi_b)_{\omega_a} \\ &\stackrel{(59)}{=} (f\psi_a - \nabla u_h \cdot \nabla \psi_a, \psi_b)_{\omega_b} - (\psi_a \nabla u_h, \nabla \psi_b)_{\omega_b}. \end{align*} $$

With (72), we then see that $(g, 1)_{\omega_b} = 0$.

It remains to consider $\psi_a \notin H_0^1(\Omega)$. In this case, we obtain again with $\text{supp}(\psi_b) \subseteq \text{supp}(\psi_a)$ and (61) with the test function $v_h = \psi_b$ (which is possible due to (53c)) that

$$ \begin{align*} (\nabla r_h^a, \nabla \psi_b)_{\omega_b} &= (\nabla r_h^a, \nabla \psi_b)_{\omega_a} \stackrel{(61)}{=} (f\psi_a - \nabla u_h \cdot \nabla \psi_a, \psi_b)_{\omega_a} - (\psi_a \nabla u_h, \nabla \psi_b)_{\omega_a} \\ &= (f\psi_a - \nabla u_h \cdot \nabla \psi_a, \psi_b)_{\omega_b} - (\psi_a \nabla u_h, \nabla \psi_b)_{\omega_b}. \end{align*} $$

With (72), we again see that $(g, 1)_{\omega_b} = 0$, which concludes the proof. $\square$

**7.4. Equilibrated flux on the large patches and the final equilibrated flux.** As in the simplified setting of Section 3.3, our equilibrated flux on the large patches $\omega_a$ and the final equilibrated flux are given by:

**Definition 7.5.** From (67), define the patchwise fluxes, for all $a \in \mathcal{V}_h$,

$$ \sigma_h^a := \sum_{b \in \mathcal{V}_h^a} \sigma_h^{a,b} \quad (73a) $$

and the equilibrated flux

$$ \sigma_h := \sum_{a \in \mathcal{V}_h} \sigma_h^a. \quad (73b) $$
**Remark 7.6.** With the definitions of Remark 6.11, the discrete residual function $r_h^a$ from (61) is simply zero, so that, up to the projection $\Upsilon_{Q_h^{a,b}}$, the local equilibrated flux $\sigma_h^{a,b}$ from (67) coincides with $\sigma_h^a$ from (73a) and (9).

We first discuss the contributions $\sigma_h^a$, which can be seen as discrete approximations of $-(\nabla r^a + \psi_a \nabla u_h)$ from (64). Recalling (37), we have.

**Lemma 7.7.** For all $a \in \mathcal{V}_h$, it holds that $\sigma_h^a \in \mathbf{H}_0(\mathrm{div}, \omega_a)$ and

$$ (\nabla \cdot \boldsymbol{\sigma}_h^a, q_h)_{\omega_a} = (f \psi_a - \nabla u_h \cdot \nabla \psi_a, q_h)_{\omega_a} \quad \text{for all } q_h \in \hat{\mathcal{Q}}_h \circ \mathbf{F}^{-1}. \quad (74) $$

*Proof.* By Definition 7.2, we have $\sigma_h^{a,b} \in \mathbf{V}_h^{a,b} \subset \mathbf{H}_0(\mathrm{div}, \omega_b)$ for all $b \in \mathcal{V}_h^a$. If $\psi_a \in H_0^1(\Omega)$, then $\psi_a \psi_b \in H_0^1(\Omega)$, and $\sigma_h^{a,b} \cdot \mathbf{n}_{\omega_b} = 0$ on $\partial \omega_b$ from (40), so that $\sigma_h^{a,b} \cdot \mathbf{n}_{\omega_a} = 0$ on $\partial \omega_a$ by virtue of $\mathrm{supp}(\psi_b) \subseteq \mathrm{supp}(\psi_a)$, which we assume in (53a). If $\psi_a \notin H_0^1(\Omega)$, then $\psi_a \psi_b \in H_0^1(\Omega)$ may still hold, in which case the previous reasoning applies. If $\psi_a \notin H_0^1(\Omega)$ and $\psi_a \psi_b \notin H_0^1(\Omega)$, then the homogeneous Neumann (no flow) boundary condition only applies on $\partial \omega_b \cap (\psi_a \psi_b)^{-1}(\{0\})$ in (40), and, in the sum of contributions over all $b \in \mathcal{V}_h^a$, on $\partial \omega_a \cap (\psi_a)^{-1}(\{0\})$. Thus $\sigma_h^a \in \mathbf{H}_0(\mathrm{div}, \omega_a)$.

To see (74), we use again definitions (67) and (73a), along with the partition of unity property (53b) and (71). Let $q_h \in \hat{\mathcal{Q}}_h \circ \mathbf{F}^{-1}$ be fixed. Then

$$ \begin{align*} (\nabla \cdot \boldsymbol{\sigma}_h^a, q_h)_{\omega_a} &= \sum_{b \in \mathcal{V}_h^a} (\nabla \cdot \boldsymbol{\sigma}_h^{a,b}, q_h)_{\omega_b} \\ &\stackrel{(67)}{=} \sum_{b \in \mathcal{V}_h^a} (\Upsilon_{Q_{h,c}^{a,b}} (f \psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b) - \nabla r_h^a \cdot \nabla \psi_b), q_h)_{\omega_b} \\ &\stackrel{(66a)}{=} \sum_{b \in \mathcal{V}_h^a} (f \psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b) - \nabla r_h^a \cdot \nabla \psi_b, q_h)_{\omega_b} \\ &\stackrel{(53b)}{=} (f \psi_a - \nabla u_h \cdot \nabla \psi_a, q_h)_{\omega_a}, \end{align*} $$

where we have crucially used Lemma 7.4 (on the second line) and the inclusion property (55c) (on the third line). $\square$

We now show that $\sigma_h$ is indeed an equilibrated flux:

**Lemma 7.8.** It holds that $\sigma_h \in \mathbf{H}(\mathrm{div}, \Omega)$ with

$$ (\nabla \cdot \boldsymbol{\sigma}_h, q_h)_{\Omega} = (f, q_h)_{\Omega} \quad \text{for all } q_h \in \hat{\mathcal{Q}}_h \circ \mathbf{F}^{-1}, \quad (75) $$

or equivalently, by virtue of (66b), $\Upsilon_{Q_h}(f - \nabla \cdot \boldsymbol{\sigma}_h) = 0$.

*Proof.* Definition 7.5, (37), and Lemma 7.7 immediately imply that $\sigma_h \in \mathbf{H}(\mathrm{div}, \Omega)$. To see the second point, we use again Lemma 7.7 and the partition of unity property (52b), giving

$$ (\nabla \cdot \boldsymbol{\sigma}_h, q_h)_{\Omega} \stackrel{(73b)}{=} \sum_{a \in \mathcal{V}_h} (\nabla \cdot \boldsymbol{\sigma}_h^a, q_h)_{\omega_a} \stackrel{(74)}{=} \sum_{a \in \mathcal{V}_h} (f \psi_a - \nabla u_h \cdot \nabla \psi_a, q_h)_{\omega_a} \stackrel{(52b)}{=} (f, q_h)_{\Omega}, $$

which concludes the proof. $\square$

## 8. A POSTERIORI ERROR ESTIMATES

We are now ready to present our a posteriori error estimates.
8.1. **Reliability.** With the dual norm

$$ \| \cdot \|_{H^{-1}(\Omega)} := \sup_{\substack{v \in H_0^1(\Omega) \\ \|\nabla v\|_{\Omega} = 1}} (\cdot, v)_{\Omega}, $$

one immediately gets the following reliability result:

**Proposition 8.1.** Let the abstract assumptions of Section 6 be satisfied, let $u$ solve (2), and let $u_h$ solve (3). Let the equilibrated flux be given by Definitions 7.1, 7.2, and 7.5. Then

$$ \|\nabla(u - u_h)\|_{\Omega} \le \|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega} + \|(1 - \Upsilon_{Q_h})(f - \nabla \cdot \boldsymbol{\sigma}_h)\|_{H^{-1}(\Omega)}. \quad (76) $$

*Proof.* The weak solution definition (2), the fact that $\boldsymbol{\sigma}_h \in \mathbf{H}(\mathrm{div}, \Omega)$, the Green theorem, and Lemma 7.8 show

$$ \begin{aligned} \|\nabla(u - u_h)\|_{\Omega} &= \sup_{\substack{v \in H_0^1(\Omega) \\ \|\nabla v\|_{\Omega} = 1}} ((f, v)_{\Omega} - (\nabla u_h, \nabla v)_{\Omega}) \\ &= \sup_{\substack{v \in H_0^1(\Omega) \\ \|\nabla v\|_{\Omega} = 1}} ((\nabla \cdot \boldsymbol{\sigma}_h, v)_{\Omega} - (\nabla u_h, \nabla v)_{\Omega} + (f - \nabla \cdot \boldsymbol{\sigma}_h, v)_{\Omega}) \\ &= \sup_{\substack{v \in H_0^1(\Omega) \\ \|\nabla v\|_{\Omega} = 1}} \left( -(\boldsymbol{\sigma}_h + \nabla u_h, \nabla v)_{\Omega} + ((1 - \Upsilon_{Q_h})(f - \nabla \cdot \boldsymbol{\sigma}_h), v)_{\Omega} \right) \\ &\le \|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega} + \|(1 - \Upsilon_{Q_h})(f - \nabla \cdot \boldsymbol{\sigma}_h)\|_{H^{-1}(\Omega)}. \end{aligned} $$

□

Let $\| \cdot \|_2$ denote the spectral norm of a square matrix. In the particular situation of Section 5.3, the $\| \cdot \|_{H^{-1}(\Omega)}$ in (76) can be further estimated by a computable weighted $L^2$-norm, forming a data oscillation term.

**Proposition 8.2.** Suppose that $\hat{Q}_h \circ \mathbf{F}^{-1}$ contains the space of piecewise constants with respect to some mesh $\mathcal{T}_h$ of $\Omega$, as is the case in Section 5.3. Assume that for all $T \in \mathcal{T}_h$, $\hat{T} := \mathbf{F}^{-1}(T)$ is convex. Then, there holds that

$$ \|\nabla(u - u_h)\|_{\Omega} \le \|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega} + \mathrm{osc}_h^{\mathrm{rel}}, \quad (77\mathrm{a}) $$

where

$$ \mathrm{osc}_h^{\mathrm{rel}} := \left( \sum_{T \in \mathcal{T}_h} \mathrm{osc}_h^{\mathrm{rel}}(T)^2 \right)^{1/2} \quad \text{with } \mathrm{osc}_h^{\mathrm{rel}}(T) := \frac{C_{\mathrm{rel}}}{\pi} \mathrm{diam}(\hat{T}) \|(1 - \Upsilon_{Q_h})(f - \nabla \cdot \boldsymbol{\sigma}_h)\|_{T} \quad (77\mathrm{b}) $$

and

$$ C_{\mathrm{rel}} := \| \det(D\mathbf{F}) \|_{\infty, \hat{\Omega}}^{1/2} \| \det(D\mathbf{F})^{-1} \|_{\infty, \hat{\Omega}}^{1/2} \sup_{\hat{x} \in \hat{\Omega}} \| D\mathbf{F}(\hat{x}) \|_{2}. \quad (77\mathrm{c}) $$

*Proof.* Let $v \in H_0^1(\Omega)$ with $\|\nabla v\|_{\Omega} = 1$ and let $v_h$ be its $L^2(\Omega)$-orthogonal projection onto the space of $\mathcal{T}_h$-piecewise constants. Recall from Lemma 7.8 that $\Upsilon_{Q_h}(f - \nabla \cdot \boldsymbol{\sigma}_h) = 0$. Hence, the Cauchy–Schwarz inequality and the Poincaré inequality as in (41) (with
constant $C_P(T))$ show that

$$((1 - \Upsilon_{Q_h})(f - \nabla \cdot \boldsymbol{\sigma}_h), v)_\Omega = (f - \nabla \cdot \boldsymbol{\sigma}_h, v - v_h)_\Omega = \sum_{T \in \mathcal{T}_h} (f - \nabla \cdot \boldsymbol{\sigma}_h, v - v_h)_T \\ \le \sum_{T \in \mathcal{T}_h} \|f - \nabla \cdot \boldsymbol{\sigma}_h\|_T \|v - v_h\|_T \le \sum_{T \in \mathcal{T}_h} \|f - \nabla \cdot \boldsymbol{\sigma}_h\|_T \mathrm{diam}(T) C_P(T) \|\nabla v\|_T \\ \le \left( \sum_{T \in \mathcal{T}_h} \mathrm{diam}(T)^2 C_P(T)^2 \|f - \nabla \cdot \boldsymbol{\sigma}_h\|_T^2 \right)^{1/2}.$$

According to [VV12, Corollary 2.2], it further holds for $\hat{T} := \mathbf{F}^{-1}(T)$ that

$$\mathrm{diam}(T)C_P(T) \le \mathrm{diam}(\hat{T}) C_P(\hat{T}) \|\det(DF)\|_{\infty, \hat{T}}^{1/2} \|\det(DF)^{-1}\|_{\infty, \hat{T}}^{1/2} \sup_{\hat{x} \in \hat{T}} \|DF(\hat{x})\|_2,$$

where $C_P(\hat{T})$ denotes the Poincaré constant of $\hat{T}$, which is smaller or equal to $1/\pi$ as we assumed that $\hat{T}$ is convex. $\square$

**Remark 8.3.** If there also holds the converse inclusion in (55c), i.e., $Q_{h,c}^{a,b} \subseteq Q_h$ (where elements in $Q_{h,c}^{a,b}$ are extended by zero outside of $\omega_b$), then $\nabla \cdot \boldsymbol{\sigma}_h \in Q_h$ and thus $(1 - \Upsilon_{Q_h})\nabla \cdot \boldsymbol{\sigma}_h = 0$. For the spaces of Section 5.3, the converse inclusion is satisfied for uniform refinements $\mathcal{T}_h$ of $\mathcal{T}_0$, but not in general. In this case, the second term in (77a) is, for smooth $f$, of order $O(h^{\bar{p}+2})$, where $h$ denotes the maximal diameter in $\mathcal{T}_h$, cf. [DEV16, Equation (3.12b)] and [EV15, Remark 3.6].

8.2. **Efficiency.** To prove efficiency, we will crucially rely on the patchwise Sobolev spaces from Section 5.2 and Assumption 6.8. We will employ the residual function $r^a \in H_*^1(\omega_a)$ from (62). The next proposition states local efficiency of the equilibrated flux estimator from Proposition 8.2:

**Proposition 8.4.** Let $\mathcal{T}_h$ be a mesh of $\Omega$, as is the case in Section 5.3. For all $T \in \mathcal{T}_h$, it holds that

$$\|\boldsymbol{\sigma}_h + \nabla u_h\|_T \le \sum_{\substack{\boldsymbol{a} \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \left( 2\sqrt{1 + (C_1 + C_2)^2} C_4 C_{\mathrm{veff}} \|\nabla(u - u_h)\|_{\omega_a} \right. \\ \left. \qquad + \sqrt{2C_4} C_{\mathrm{veff}} \mathrm{osc}_h^{\mathrm{eff}}(\omega_a, T) \right), \quad (78a)$$

where

$$\mathrm{osc}_h^{\mathrm{eff}}(\omega_a, T) := \left( \sum_{\substack{\boldsymbol{b} \in \mathcal{V}_h^a \\ |\omega_b \cap T| > 0}} \left( \mathrm{diam}(\omega_b)^2 C_{\mathrm{PF}}(\omega_b)^2 \|(1 - \Upsilon_{Q_{h,c}^{a,b}})(f \psi_a \psi_b - \nabla u_h \cdot \nabla(\psi_a \psi_b))\|_{\omega_b}^2 \right. \right. \\ \left. \left. + \|(1 - \Pi_{RT_h^{a,b}})(\psi_a \psi_b \nabla u_h)\|_{\omega_b}^2 \right) \right)^{1/2} \quad (78b)$$

and $\Pi_{RT_h^{a,b}} : [L^2(\omega_b)]^d \to RT_h^{a,b}$ denotes the $L^2$-orthogonal projection. With $C_F := \|\det(DF)\|_{\infty, \hat{\Omega}} \|\det(DF)^{-1}\|_{\infty, \hat{\Omega}}$, the constant $C_{\mathrm{veff}}$ is explicitly given as

$$C_{\mathrm{veff}} := 2C_{\mathrm{st}} \max \{C_3 C_5 + C_1(C_5 + C_6), C_5 + (C_5 + C_6) + C_F C_6\}. \quad (78c)$$

In the setting of Sections 2–5, all involved constants depend themselves only on the space dimension $d$, the polynomial degree $\bar{p}$ from Section 5.1 (which itself depends only on the
considered smoothness, see Remark 6.10), and $\max\{\|DF\|_{\infty,\hat{\Omega}}, \|(DF)^{-1}\|_{\infty,\hat{\Omega}}\}$. They do not depend on the polynomial degrees $p$ and $\tilde{p}$.

*Proof.* We prove the assertion in seven steps.

**Step 1:** Definition (73b), the partition of unity property (52b), and the triangle inequality show that

$$ \|\sigma_h + \nabla u_h\|_T \stackrel{(73b)}{=} \stackrel{(52b)}{\left\|} \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} (\sigma_h^a + \psi_a \nabla u_h) \stackrel{T}{\left\|} \le \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \|\sigma_h^a + \psi_a \nabla u_h\|_T. $$

**Step 2:** Next, we bound each summand separately. Let $a \in \mathcal{V}_h$ with $|\omega_a \cap T| > 0$. Then, definition (73a) and the partition of unity property (53b) together with the Cauchy-Schwarz inequality give that

$$ \begin{aligned} \|\sigma_h^a + \psi_a \nabla u_h\|_T &\stackrel{(73a)}{=} \stackrel{(53b)}{\left\|} \sum_{\substack{b \in \mathcal{V}_h^a \\ |\omega_b \cap T| > 0}} (\sigma_h^{a,b} + \psi_a \psi_b \nabla u_h) \stackrel{T}{\left\|} \\ &\le \stackrel{(54d)}{\sqrt{C_4}} \left( \sum_{\substack{b \in \mathcal{V}_h^a \\ |\omega_b \cap T| > 0}} \|\sigma_h^{a,b} + \psi_a \psi_b \nabla u_h\|_T^2 \right)^{1/2} \\ &\le \stackrel{(54e)}{\sqrt{2C_4}} \left( \sum_{\substack{b \in \mathcal{V}_h^a \\ |\omega_b \cap T| > 0}} (\|\sigma_h^{a,b} + \psi_b(\psi_a \nabla u_h + \nabla r_h^a)\|_T^2 + C_5^2 \|\nabla r_h^a\|_{\omega_b}^2) \right)^{1/2}. \end{aligned} \quad (79) $$

**Step 3:** We estimate the first summand of (79). To this end, let $b \in \mathcal{V}_h^a$ with $|\omega_b \cap T| > 0$. As in Lemma 7.3, we abbreviate $g := f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b) - \nabla r_h^a \cdot \nabla \psi_b$ as well as $\tau := \psi_b(\psi_a\nabla u_h + \nabla r_h^a)$. Then, definition (67) (or its equivalent formulation (70), which allows to stick in the projector $\Pi_{RT_h^{a,b}}$) with $V_h^{a,b} \subseteq RT_h^{a,b}$ together with (58) give that

$$ \begin{aligned} \|\sigma_h^{a,b} + \psi_b(\psi_a \nabla u_h + \nabla r_h^a)\|_T &\le \|\sigma_h^{a,b} + \tau\|_{\omega_b} \\ &\stackrel{(67)}{=} \min_{\substack{v_h \in \mathcal{V}_h^{a,b} \\ \nabla \cdot v_h = \Upsilon_{Q_h^{a,b}(g)}}} \|v_h + \Pi_{RT_h^{a,b}}(\tau)\|_{\omega_b} \stackrel{(58)}{\le} C_{\text{st}} \min_{\substack{v \in \mathcal{H}_0(\text{div},\omega_b) \\ \nabla \cdot v = \Upsilon_{Q_h^{a,b}(g)}}} \|v + \Pi_{RT_h^{a,b}}(\tau)\|_{\omega_b} \\ &\stackrel{(54e)}{\le} C_{\text{st}} \left( \min_{\substack{v \in \mathcal{H}_0(\text{div},\omega_b) \\ \nabla \cdot v = \Upsilon_{Q_h^{a,b}(g)}}} \|v + \tau\|_{\omega_b} + \|(1 - \Pi_{RT_h^{a,b}})(\psi_a \psi_b \nabla u_h)\|_{\omega_b} + C_5 \|\nabla r_h^a\|_{\omega_b} \right). \end{aligned} $$

Let $r^{a,b} \in H_*(\omega_b)$ solve

$$ (\nabla r^{a,b}, \nabla v)_{\omega_b} = (\Upsilon_{Q_h^{a,b}(g)}, v)_{\omega_b} - (\tau, \nabla v)_{\omega_b} \quad \text{for all } v \in H_*(\omega_b). $$

Then a standard primal-dual equivalence gives, as in, e.g., [EV20, Corollary 3.6],

$$ \min_{\substack{v \in \mathcal{H}_0(\text{div}, \omega_b) \\ \nabla \cdot v = \Upsilon_{Q_h^{a,b}(g)}}} \|v + \tau\|_{\omega_b} = \|\nabla r^{a,b}\|_{\omega_b} = \sup_{\substack{v \in H_*(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} ((\Upsilon_{Q_h^{a,b}(g)}, v)_{\omega_b} - (\tau, \nabla v)_{\omega_b}). $$
Thus,

$$
\begin{align}
\min_{\substack{v \in H_0(\mathrm{div}, \omega_b) \\ \nabla \cdot v = \Upsilon_{Q_h^{a,b}}(g)}} \|v + \tau\|_{\omega_b} &= \sup_{\substack{v \in H_0^*(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \left( -(\tau, \nabla v)_{\omega_b} + (g, v)_{\omega_b} - ((1 - \Upsilon_{Q_h^{a,b}})g, v)_{\omega_b} \right) \nonumber \\
&= \sup_{\substack{v \in H_0^*(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \left( (f\psi_a - \nabla u_h \cdot \nabla \psi_a, v\psi_b)_{\omega_b} - (\psi_a \nabla u_h + \nabla r_h^a, \nabla(v\psi_b))_{\omega_b} \right. \tag{80} \\
&\qquad \left. - ((1 - \Upsilon_{Q_h^{a,b}})(f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b) - \nabla r_h^a \cdot \nabla \psi_b), v)_{\omega_b} \right). \tag{81}
\end{align}
$$

**Step 4:** Recalling (39), we note that for all $v \in H_0^1(\omega_b)$,

$$
\begin{align*}
(f\psi_a, v\psi_b)_{\omega_b} &= (f, \underbrace{v\psi_a\psi_b}_{\in H_0^1(\Omega)})_{\omega_b} \stackrel{(2)}{=} (\nabla u, \nabla(v\psi_a\psi_b))_{\omega_b} \\
&= (\nabla u \cdot \nabla \psi_a, v\psi_b)_{\omega_b} + (\psi_a \nabla u, \nabla(v\psi_b))_{\omega_b}.
\end{align*}
$$

With the Poincaré–Friedrichs inequality (41) and Assumption 6.5, the latter equality
allows us to further estimate the term in (80)

$$
\begin{align*}
\sup_{\substack{v \in H_0^1(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \left( (f\psi_a - \nabla u_h \cdot \nabla \psi_a, v\psi_b)_{\omega_b} - (\psi_a \nabla u_h + \nabla r_h^a, \nabla(v\psi_b))_{\omega_b} \right) \\
&= \sup_{\substack{v \in H_0^1(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \left( (\nabla(u-u_h) \cdot \nabla \psi_a, v\psi_b)_{\omega_b} + (\psi_a \nabla(u-u_h) - \nabla r_h^a, \nabla(v\psi_b))_{\omega_b} \right) \\
&\le \|\nabla(u-u_h)\|_{\omega_b} \|\nabla \psi_a\|_{\infty, \omega_b} \sup_{\substack{v \in H_0^1(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \|v\psi_b\|_{\omega_b} \\
&\quad + (\|\psi_a\|_{\infty, \omega_b} \|\nabla(u-u_h)\|_{\omega_b} + \|\nabla r_h^a\|_{\omega_b}) \sup_{\substack{v \in H_0^1(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \|\nabla(v\psi_b)\|_{\omega_b} \\
&\le \|\nabla(u-u_h)\|_{\omega_b} (\|\nabla \psi_a\|_{\infty, \omega_b} \|\psi_b\|_{\infty, \omega_b} \operatorname{diam}(\omega_b) C_{\mathrm{PF}}(\omega_b)) \\
&\quad + (\|\psi_a\|_{\infty, \omega_b} \|\nabla(u-u_h)\|_{\omega_b} + \|\nabla r_h^a\|_{\omega_b}) (\|\nabla \psi_b\|_{\infty, \omega_b} \operatorname{diam}(\omega_b) C_{\mathrm{PF}}(\omega_b) + \|\psi_b\|_{\infty, \omega_b}) \\
&\stackrel{6.5}{\le} (C_3 C_5 + C_1(C_5 + C_6)) \|\nabla(u-u_h)\|_{\omega_b} + (C_5 + C_6) \|\nabla r_h^a\|_{\omega_b}.
\end{align*}
$$

**Step 5:** To estimate the term in (81), we use that $\Upsilon_{Q_h^{a,b}}g = \Upsilon_{Q_{h,c}^{a,b}}g$ from Lemma 7.4,
(44), (66a), the Poincaré–Friedrichs inequality (41), and (54f) to infer that

$$
\begin{align*}
\sup_{\substack{v \in H_0^1(\omega_b) \\ \|\nabla v\|_{\omega_b} = 1}} \left( (1 - \Upsilon_{Q_h^{a,b}})(f\psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b) - \nabla r_h^a \cdot \nabla \psi_b), v \right)_{\omega_b} & \\
&\stackrel{(44)}{\le} \left( \| (1 - \Upsilon_{Q_{h,c}^{a,b}})(f\psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b)) \|_{\omega_b} \right. \\
&\qquad \left. + \| \det(D\mathbf{F}) \|_{\infty, \hat{\omega}_b} \| \det(D\mathbf{F})^{-1} \|_{\infty, \hat{\omega}_b} \| \nabla r_h^a \|_{\omega_b} \| \nabla \psi_b \|_{\infty, \omega_b} \right) \operatorname{diam}(\omega_b) C_{\mathrm{PF}}(\omega_b) \\
&\stackrel{(54\mathrm{f})}{\le} \operatorname{diam}(\omega_b) C_{\mathrm{PF}}(\omega_b) \| (1 - \Upsilon_{Q_{h,c}^{a,b}})(f\psi_a \psi_b - \nabla u_h \cdot \nabla (\psi_a \psi_b)) \|_{\omega_b} + C_F C_6 \| \nabla r_h^a \|_{\omega_b}.
\end{align*}
$$
Together with Steps 3–4 and the Cauchy–Schwarz inequality, we conclude that

$$
\begin{align}
& \|\sigma_h^{a,b} + \psi_b(\psi_a \nabla u_h + \nabla r_h^a)\|_T^2 \nonumber \\
& \le C_{\text{veff}}^2 \left( \|\nabla(u-u_h)\|_{\omega_b}^2 + \|\nabla r_h^a\|_{\omega_b}^2 + \|(1-\Pi_{RT_h^{a,b}})(\psi_a \psi_b \nabla u_h)\|_{\omega_b}^2 \right. \tag{82} \\
& \qquad \left. + \text{diam}(\omega_b)^2 C_{\text{PF}}(\omega_b)^2 \|(1-\Upsilon_{Q_{h,c}^{a,b}})(f\psi_a \psi_b - \nabla u_h \cdot \nabla(\psi_a \psi_b))\|_{\omega_b}^2 \right). \nonumber
\end{align}
$$

**Step 6:** As a final auxiliary step, we bound $\|\nabla r_h^a\|_{\omega_a}$. Since $r_h^a$ from (61) is the Galerkin approximation of $r^a$ from (62), and crucially employing (36), we see with the variational formulation (2) that

$$
\begin{align*}
\|\nabla r_h^a\|_{\omega_a} &\le \|\nabla r^a\|_{\omega_a} \stackrel{(62)}{=} \sup_{\substack{v \in H_0^1(\Omega) \\ \|\nabla v\|_{\omega_a} = 1}} \left( (f, \underbrace{v\psi_a}_{v \in H_0^1(\Omega)})_{\omega_a} - (\nabla u_h, \nabla(v\psi_a))_{\omega_a} \right) \\
&\stackrel{(2)}{=} \sup_{\substack{v \in H_0^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} (\nabla(u-u_h), \nabla(v\psi_a))_{\omega_a} \le \|\nabla(u-u_h)\|_{\omega_a} \sup_{\substack{v \in H_0^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} \|\nabla(v\psi_a)\|_{\omega_a}.
\end{align*}
$$

Relying on the Poincaré inequality (41) and with (54b) and (54a), we can bound the supremum

$$
\sup_{\substack{v \in H_0^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} \|\nabla(v\psi_a)\|_{\omega_a} \le \sup_{\substack{v \in H_0^1(\omega_a) \\ \|\nabla v\|_{\omega_a} = 1}} \left( \|\nabla\psi_a\|_{\infty, \omega_a} \|v\|_{\omega_a} + \|\psi_a\|_{\infty, \omega_a} \|\nabla v\|_{\omega_a} \right) \stackrel{(54a)}{\le} C_1 + C_2. \tag{83}
$$

**Step 7:** Putting all steps together, also using that $C_5 \le C_{\text{veff}}$, we obtain that

$$
\begin{align*}
& \|\sigma_h + \nabla u_h\|_T \\
&\le \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \sqrt{2C_4} \left( \sum_{\substack{b \in \mathcal{V}_h^a \\ |\omega_b \cap T| > 0}} \left( \|\sigma_h^{a,b} + \psi_b(\psi_a \nabla u_h + \nabla r_h^a)\|_T^2 + C_5^2 \|\nabla r_h^a\|_{\omega_b}^2 \right) \right)^{1/2} \\
&\stackrel{(82)}{\le} \sqrt{2C_4} \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \left( \sum_{\substack{b \in \mathcal{V}_h^a \\ |\omega_b \cap T| > 0}} \left( C_{\text{veff}}^2 \left( \|\nabla(u-u_h)\|_{\omega_b}^2 + \|\nabla r_h^a\|_{\omega_b}^2 \right. \right. \right. \\
&\qquad \left. \left. + \|(1-\Pi_{RT_h^{a,b}})(\psi_a \psi_b \nabla u_h)\|_{\omega_b}^2 \right. \right. \\
&\qquad \left. + \text{diam}(\omega_b)^2 C_{\text{PF}}(\omega_b)^2 \|(1-\Upsilon_{Q_{h,c}^{a,b}})(f\psi_a \psi_b - \nabla u_h \cdot \nabla(\psi_a \psi_b))\|_{\omega_b}^2 \right) + C_5^2 \|\nabla r_h^a\|_{\omega_b}^2 \right)^{1/2} \\
&\stackrel{(54d)}{\le} 2C_4 C_{\text{veff}} \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \left( \|\nabla(u-u_h)\|_{\omega_a}^2 + \|\nabla r_h^a\|_{\omega_a}^2 \right)^{1/2} + \sqrt{2C_4} C_{\text{veff}} \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \text{osc}_h^{\text{eff}}(\omega_a, T) \\
&\stackrel{(78b)}{\le} 2\sqrt{1+(C_1+C_2)^2} C_4 C_{\text{veff}} \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \|\nabla(u-u_h)\|_{\omega_a} + \sqrt{2C_4} C_{\text{veff}} \sum_{\substack{a \in \mathcal{V}_h \\ |\omega_a \cap T| > 0}} \text{osc}_h^{\text{eff}}(\omega_a, T),
\end{align*}
$$

which concludes the proof. $\square$
**Remark 8.5.** In the situation of Section 5, $\mathrm{osc}_h^{\mathrm{eff}}(\omega_a, T)$ of (78b) vanishes if $\mathbf{F}$ is affine, $\mathcal{T}_h$ is a uniform refinement of $\mathcal{T}_0$, $f$ is a $\mathcal{T}_h$-piecewise polynomial of some degree $q = (q, \dots, q)$ with $q \ge 0$, and $\tilde{p} \ge \max\{q + \bar{p} + 1, p + \bar{p} + 1\}$. Note that $\tilde{p} \ge p + \bar{p} + 1$ in particular implies that $\psi_a\psi_b\nabla u_h + \psi_b\nabla r_h^a \in \mathbf{RT}^{\tilde{p}}(\mathcal{T}_b)$ that we take for $\mathbf{RT}_h^{a,b}$.

We now turn towards the global efficiency. In order to achieve robustness with respect to the strength of the hierarchical refinement (the number of hanging nodes), we do not straightforwardly use the element-related result of Proposition 8.4 but rather resort to its patch-related variant. With an assumption on the maximal overlap by the patches $\omega_a$ (not limiting the strength of the hierarchical refinement, see Remark 8.7), our global efficiency result is:

**Proposition 8.6.** Let $\mathcal{T}_h$ be a mesh of $\Omega$, as is the case in Section 5.3. Let $C_{\mathrm{over}} > 0$ be a constant such that

$$ \sup_{x \in \Omega} \#\{a \in \mathcal{V}_h : x \in \omega_a\} \le C_{\mathrm{over}}. \quad (84a) $$

Then, there holds that

$$ \begin{aligned} \|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega}^2 \le 4C_{\mathrm{over}}^2 C_4^2 C_{\mathrm{veff}}^2 ((C_1 + C_2)^2 + 1) \|\nabla(u - u_h)\|_{\Omega}^2 \\ + 4C_{\mathrm{over}} C_4 C_{\mathrm{veff}}^2 \sum_{a \in \mathcal{V}_h} \mathrm{osc}_h^{\mathrm{eff}}(\omega_a)^2, \end{aligned} \quad (84b) $$

where

$$ \begin{aligned} \mathrm{osc}_h^{\mathrm{eff}}(\omega_a)^2 := & \sum_{b \in \mathcal{V}_h^a} (\mathrm{diam}(\omega_b)^2 C_{\mathrm{PF}}(\omega_b)^2 \|(1 - \Upsilon_{Q_{h,c}^{a,b}})(f\psi_a\psi_b - \nabla u_h \cdot \nabla(\psi_a\psi_b))\|_{\omega_b}^2 \\ & + \|(1 - \Pi_{RT_h^{a,b}})(\psi_a\psi_b\nabla u_h)\|_{\omega_a}^2). \end{aligned} \quad (84c) $$

In the setting of Sections 2–5, all involved constants except of $C_{\mathrm{over}}$ depend themselves only on the space dimension $d$, the polynomial degree $\bar{p}$ from Section 5.1 (which itself depends only on the considered smoothness, see Remark 6.10), and $\max\{\|DF\|_{\infty,\hat{\Omega}}, \|(DF)^{-1}\|_{\infty,\hat{\Omega}}\}$. They do not depend on the polynomial degrees $p$ and $\tilde{p}$. $C_{\mathrm{over}}$ is discussed in Remark 8.7 below.

**Proof.** Proceeding as in the proof of Proposition 8.4 while relying on Definition 7.5 and the partitions of unity (52b) and (53b) together with the finite overlap assumptions (54d) and (84a), we see

$$ \begin{aligned} \|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega}^2 &\stackrel{(52b)}{=} \left\| \sum_{a \in \mathcal{V}_h} (\boldsymbol{\sigma}_h^a + \psi_a \nabla u_h) \right\|_{\Omega}^2 \stackrel{(84a)}{\le} C_{\mathrm{over}} \sum_{a \in \mathcal{V}_h} \|\boldsymbol{\sigma}_h^a + \psi_a \nabla u_h\|_{\omega_a}^2 \\ &\stackrel{(53b)}{=} C_{\mathrm{over}} \sum_{a \in \mathcal{V}_h} \left\| \sum_{b \in \mathcal{V}_h^a} (\boldsymbol{\sigma}_h^{a,b} + \psi_a \psi_b \nabla u_h) \right\|_{\omega_a}^2 \\ &\stackrel{(54d)}{\le} C_{\mathrm{over}} C_4 \sum_{a \in \mathcal{V}_h} \sum_{b \in \mathcal{V}_h^a} \|\boldsymbol{\sigma}_h^{a,b} + \psi_a \psi_b \nabla u_h\|_{\omega_b}^2. \end{aligned} $$

We now employ (82) plus the triangle and the Cauchy–Schwarz inequalities. Also using (54e) and the fact that $C_5 \le C_{\mathrm{veff}}$, this leads to

$$ \|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega}^2 \le C_{\mathrm{over}} C_4 C_{\mathrm{veff}}^2 \sum_{a \in \mathcal{V}_h} \left( \sum_{b \in \mathcal{V}_h^a} (\|\nabla(u - u_h)\|_{\omega_b}^2 + \|\nabla r_h^a\|_{\omega_b}^2) + \mathrm{osc}_h^{\mathrm{eff}}(\omega_a)^2 \right), $$
and we are left to treat the first two terms. The finite overlap assumptions (54d) and (84a) again imply

$$ \sum_{a \in \mathcal{V}_h} \sum_{b \in \mathcal{V}_h^a} \|\nabla(u - u_h)\|_{\omega_b}^2 \stackrel{(54d)}{\le} C_4 \sum_{a \in \mathcal{V}_h} \|\nabla(u - u_h)\|_{\omega_a}^2 \stackrel{(84a)}{\le} C_{\text{over}} C_4 \|\nabla(u - u_h)\|_{\Omega}^2. $$

Reasoning similarly and also employing the two estimates from Step 6 of the proof of Proposition 8.4, we see

$$ \sum_{a \in \mathcal{V}_h} \sum_{b \in \mathcal{V}_h^a} \|\nabla r_h^a\|_{\omega_b}^2 \stackrel{(54d)}{\le} C_4 \sum_{a \in \mathcal{V}_h} \|\nabla r_h^a\|_{\omega_a}^2 \le C_4(C_1 + C_2)^2 \sum_{a \in \mathcal{V}_h} \|\nabla(u - u_h)\|_{\omega_a}^2 \\ \stackrel{(84a)}{\le} C_{\text{over}} C_4(C_1 + C_2)^2 \|\nabla(u - u_h)\|_{\Omega}^2, $$

which altogether gives (84b). □

**Remark 8.7.** For (scaled) hierarchical B-splines $\psi_a$ of degree $\bar{p}$ on $\overline{\mathcal{H}}$-admissible (i.e., $\mathcal{H}$-admissible with respect to $\bar{p}$, $\overline{\mathcal{K}}_0$, and $\overline{m}$) meshes $\hat{T}_h$ of class $\mu$ introduced along with suitable refinement algorithms in [GHP17] for $\mu = 2$ and in [BGV18] for $\mu \ge 2$, the upper bound $C_{\text{over}}$ from (84a) only depends on the polynomial degree $\bar{p}$ (which itself only depends on the considered smoothness), the initial knot vector $\overline{\mathcal{K}}_0$, the multiplicity $\overline{m}$, and the grading parameter $\mu$, allowing for meshes with arbitrarily many hanging nodes. Meshes satisfying assumption (84a) with the constant $C_{\text{over}}$ only depending on the polynomial degree $\bar{p}$ are considered in the numerics Section 9 below.

**Remark 8.8.** We further mention that the degree $\bar{p}$ and the smoothness $\bar{p}-\overline{m}$ for (scaled) hierarchical B-splines $\psi_a$, or more precisely the approximation power of the spanned space, do not directly affect the oscillations (84c) and thus the efficiency of the estimator. Similarly, the chosen $\bar{p}$ and $\overline{m}$ do not directly affect the oscillations (77b) and thus the reliability of the estimator. Indeed, as discussed in Remarks 5.1 and 8.5, the presence and size of data oscillations is rather induced by the approximation properties of the spaces $Q_h^{a,b}$ and $RT^{\tilde{p}}(T_b)$ which are discontinuous (piecewise with respect to $T_b$).

## 9. NUMERICAL EXPERIMENTS

We consider problem (1) on the quarter ring depicted in Figure 7,

$$ \Omega := \{r(\cos(\varphi), \sin(\varphi)) : r \in (1/2, 1) \wedge \varphi \in (0, \pi/2)\} $$

with NURBS parametrization $\mathbf{F}$ as in [GHP17, Section 6.3], and prescribe the exact solution

$$ u(x, y) = xy \sin(4\pi(x^2 + y^2)). \qquad (85) $$

For polynomial degrees $p \in \{1, \dots, 5\}$ and multiplicities $m \in \{1, p\}$, we define the initial knot vectors by

$$ \mathcal{K}_{1(0)} := \mathcal{K}_{2(0)} := \left( \underbrace{0, \dots, 0}_{(p+1)\text{-times}}, \underbrace{1/2, \dots, 1/2}_{m\text{-times}}, \underbrace{1, \dots, 1}_{(p+1)\text{-times}} \right), $$

leading to piecewise $p$-degree polynomials with $C^{p-m}$ smoothness. As the corresponding polynomial degree $\bar{p}$ for the partition of unity by the $\psi_a$ in Section 5.1, we choose $\bar{p} := p+1-m$ and the corresponding multiplicity $\overline{m} := \bar{p}-p+m = 1$ following (31), so that the $\psi_a$ are mapped piecewise $\bar{p}$-degree polynomials of class $C^{\bar{p}-1}$. The polynomial degree $\tilde{p}$ for the flux equilibration in Section 5.3 is chosen in $\{p+1, p+2\}$. Following Remark 8.5, ignoring temporarily the source term $f$, this would only imply $\text{osc}_h^{\text{eff}}(\omega_a, T) = \text{osc}_h^{\text{eff}}(\omega_a) = 0$ in
<figure><img src="image_8.png" /><figcaption>FIGURE 7. The domain $\Omega$ considered in Section 9 along with the initial mesh $T_0$ and the first three refinements: in the right half of $\Omega$ (leading to arbitrary many hanging nodes) (top) and towards the point (1, 0) (violating the bounded patch overlap condition (84a)) (bottom).</figcaption></figure>

Propositions 8.4 and 8.6 if $F$ was affine (which is not the case here) and $\tilde{p} \ge p + \bar{p} + 1$ (which is only the case here for $m = p$, i.e., $C^0$, but not higher-smoothness splines; note that in the case of maximum $C^{p-1}$ smoothness, we should use $\tilde{p} \ge p + \bar{p} + 1 = 2p + 1$ theoretically, but we merely employ $\tilde{p} = p + 1$ or $\tilde{p} = p + 2$ numerically). Nevertheless, both choices $\tilde{p} \in \{p+1, p+2\}$ seem to perform numerically well in the considered test case also for high-smoothness cases, up to $m = 1$, corresponding to $C^{p-1}$ splines.

We consider four different refinements of the initial mesh: 1) uniform refinement, where in each step, all elements in the parameter domain are bisected in both directions; 2) adaptive refinement, where in each step, a minimal set of elements $M_h \subseteq T_h$ is marked via the Dörfler marking

$$ \theta \sum_{T \in \mathcal{T}_h} (\|\boldsymbol{\sigma}_h + \nabla u_h\|_T^2 + \mathrm{osc}_h^{\mathrm{rel}}(T)^2) \le \sum_{T \in \mathcal{M}_h} (\|\boldsymbol{\sigma}_h + \nabla u_h\|_T^2 + \mathrm{osc}_h^{\mathrm{rel}}(T)^2) $$

with $\theta = 0.5$ and subsequently refined via the refinement strategy from [GHP17] (see also Remark 8.7); 3) artificial refinement enforcing an arbitrary number of hanging nodes, where in each step, all elements in the parameter domain that are contained in $[0, 1/2] \times [0, 1]$ are bisected in both directions (see Figure 7, top); 4) artificial refinement enforcing an arbitrary number of overlapping patches $\omega_a$, where in each step, the element in the parameter domain containing the point (0, 0), which is mapped onto (1, 0) under $F$, is bisected in both directions (see Figure 7, bottom). In each case, new knots have multiplicity $m$.

The resulting effectivity indices $\frac{\|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega}}{\|\nabla(u-u_h)\|_{\Omega}}$ and $\frac{\|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega} + \mathrm{osc}_h^{\mathrm{rel}}}{\|\nabla(u-u_h)\|_{\Omega}}$ as function of the number of mesh elements $N$ in $T_h$ are displayed in Figures 8–11. Recall from Propositions 8.4 and 8.6 that the efficiency constant in (84b) may theoretically depend on the space dimension $d$, the polynomial degree $\bar{p}$ (which itself depends on the considered smoothness, see Remark 6.10), $\max\{\|DF\|_{\infty, \hat{\Omega}}, \|(DF)^{-1}\|_{\infty, \hat{\Omega}}\}$, and the overlap constant $C_{\mathrm{over}}$ from (84a) (which itself only depends on $\bar{p}$ for the first three refinements types, see Remark 8.7, but grows unboundedly in the fourth case). At least in this example, though, the dependence on $\bar{p}$ is not observed, and the equilibrated flux estimator $\|\boldsymbol{\sigma}_h + \nabla u_h\|_{\Omega}$ (plus oscillation
<figure><img src="image_9.png" /><figcaption>FIGURE 8. Effectivity indices $\frac{\|\sigma_h + \nabla u_h\|_{\Omega}}{\|\nabla(u-u_h)\|_{\Omega}}$ (left) and $\frac{\|\sigma_h + \nabla u_h\|_{\Omega} + \text{osc}_h^{\text{rel}}}{\|\nabla(u-u_h)\|_{\Omega}}$ (right) corresponding to the problem of Section 9 with uniform mesh refinement, polynomial degrees $p \in \{1, \dots, 5\}$, multiplicities $m \in \{1, p\}$.</figcaption></figure>

terms) seems to be not only, as proven, robust with respect to the polynomial degrees $p$ and $\tilde{p}$, but also with respect to the smoothness incarnated in $\tilde{p}$. The increase of the effectivity indices on adaptively refined meshes with $\tilde{p} = p+1$ in the right part of Figure 9 is only because of data oscillation, as we discuss below. Moreover, as theoretically shown in Proposition 8.6, we also numerically observe in Figure 10 the robustness with respect to the strength of the hierarchical refinement (number of hanging nodes). For the fourth refinement, (84a) is not uniformly satisfied, which, however, is again not reflected in the resulting efficiency constants. Note however that the computation of the equilibrated flux becomes very expensive in this case, as there are nodes $a$ whose patches $\omega_a$ coincide with the entire computational domain $\Omega$, so that the corresponding meshes $T_a$ are uniform refinements of $T_0$ up to the same level as the element containing the point (1, 0). Additional numerical experiments were carried out for an exact solution $u(x, y) = \sin(2\pi x) \sin(2\pi y)$ on the square (not displayed), with similar results.

The “data oscillation” terms $\text{osc}_h^{\text{rel}}$ from (77b) are displayed separately in Figure 12, again as function of the number of mesh elements $N$ in $T_h$. For the present smooth solution (85), one expects the error $\|\nabla(u - u_h)\|_{\Omega}$ to decay as $O(h^p) \approx O(N^{-p/2})$ for uniform mesh refinement. Recall from Remark 8.3 that $\text{osc}_h^{\text{rel}}$ are expected to decay in this case as $O(h^{\tilde{p}+2})$, i.e., as $O(N^{-p/2-3/2})$ or $O(N^{-p/2-2})$ for respectively $\tilde{p} = p+1$ and $\tilde{p} = p+2$. Figure 12 only concerns adaptive mesh refinement, where $\|\nabla(u - u_h)\|_{\Omega}$ is still expected to decay as $O(N^{-p/2})$. We do not have here theoretical indications for $\text{osc}_h^{\text{rel}}$, but we observe at least $O(N^{-p/2})$ in most cases, even though the chosen $\tilde{p}$ is theoretically inappropriate for higher smoothness as discussed above. For uniform mesh refinement (not displayed), we indeed observe $O(h^{\tilde{p}+2})$, following Remark 8.3.

## APPENDIX A. PROOF OF THE BROKEN POLYNOMIAL EXTENSION PROPERTY

In this section, we verify Assumption 6.8 for the spaces defined in Section 5.3. Recall the Piola transformations from (43) and (44)

$$ \Phi(\cdot) := (\det(DF)^{-1}(DF)(\cdot)) \circ F^{-1} \quad \text{and} \quad \tilde{\Phi}(\cdot) := (\det(DF)^{-1}(\cdot)) \circ F^{-1} $$
<figure><img src="image_10.png" /><figcaption>FIGURE 9. Effectivity indices $\frac{\|\sigma_h + \nabla u_h\|_{\Omega}}{\|\nabla(u-u_h)\|_{\Omega}}$ (left) and $\frac{\|\sigma_h + \nabla u_h\|_{\Omega} + \text{osc}_h^{\text{rel}}}{\|\nabla(u-u_h)\|_{\Omega}}$ (right) corresponding to the problem of Section 9 with adaptive mesh refinement, polynomial degrees $p \in \{1, \dots, 5\}$, multiplicities $m \in \{1, p\}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIGURE 10. Effectivity indices $\frac{\|\sigma_h + \nabla u_h\|_{\Omega}}{\|\nabla(u-u_h)\|_{\Omega}}$ (left) and $\frac{\|\sigma_h + \nabla u_h\|_{\Omega} + \text{osc}_h^{\text{rel}}}{\|\nabla(u-u_h)\|_{\Omega}}$ (right) corresponding to the problem of Section 9 with artificial refinement enforcing a high number of hanging nodes (see Figure 7, top), polynomial degrees $p \in \{1, \dots, 5\}$, multiplicities $m \in \{1, p\}$.</figcaption></figure>

and the identity (45)

$$ \tilde{\Phi}(\nabla(\cdot)) = \nabla \cdot \Phi(\cdot); \qquad (86) $$

see, e.g., [EG21, Lemma 9.6]. Following the imposition of the Neumann boundary condition in the space $H_0(\text{div}, \omega_b)$ in (40), we denote by $\mathcal{F}_b^N$ the boundary faces of the local mesh $\mathcal{T}_b$ if $\psi_a\psi_b \in H_0^1(\Omega)$ and such boundary faces of $\mathcal{T}_b$ where $F \subset (\psi_a\psi_b)^{-1}(\{0\})$ for $\psi_a\psi_b \notin H_0^1(\Omega)$. By $\mathcal{F}_b^{\text{int}}$, we denote the interior faces of $\mathcal{T}_b$. Analogously, we write $\hat{\mathcal{F}}_b^N$ and $\hat{\mathcal{F}}_b^{\text{int}}$ for the corresponding faces on the parameter mesh $\hat{\mathcal{T}}_b$. Finally, we write $\nabla_b(\cdot)$ for the $\mathcal{T}_b$- or $\hat{\mathcal{T}}_b$-piecewise divergence operator.
<figure><img src="image_12.png" /><figcaption>FIGURE 11. Effectivity indices $\frac{\|\sigma_h + \nabla u_h\|_{\Omega}}{\|\nabla(u-u_h)\|_{\Omega}}$ (left) and $\frac{\|\sigma_h + \nabla u_h\|_{\Omega} + \text{osc}_{h}^{\text{rel}}}{\|\nabla(u-u_h)\|_{\Omega}}$ (right) corresponding to the problem of Section 9 with artificial refinement enforcing a high number of overlapping patches (see Figure 7, bottom), polynomial degrees $p \in \{1, \dots, 5\}$, multiplicities $m \in \{1, p\}$.</figcaption></figure>

Since $\tau_h \in RT_h^{a,b} = RT^{\tilde{p}}(\mathcal{T}_b)$ and $V_h^{a,b} = RT^{\tilde{p}}(\mathcal{T}_b) \cap H_0(\text{div}, \omega_b)$ by (49), we can substitute in (58) $v_h + \tau_h = w_h$ with $w_h = \Phi(\hat{w}_h)$ and $\hat{w}_h \in RT^{\tilde{p}}(\hat{\mathcal{T}}_b)$, also using (50). This shows that

$$ \min_{\substack{v_h \in V_h^{a,b} \\ \nabla \cdot v_h = g_h}} \|v_h + \tau_h\|_{\omega_b} = \min_{\substack{w_h \in RT^{\tilde{p}}(\mathcal{T}_b) \\ w_h \cdot n_{\omega_b} |_{F} = \tau_h \cdot n_{\omega_b} |_{F} = r_F \\ [w_h \cdot n_{\omega_b}]_F = [\tau_h \cdot n_{\omega_b}]_F = r_F \\ \nabla_b \cdot w_h|_{T} = (g_h + \nabla_b \cdot \tau_h)|_T = r_T \\ \nabla \cdot w_h|_{\partial T} = r_T}} \|w_h\|_{\omega_b} = \min_{\substack{\hat{w}_h \in RT^{\tilde{p}}(\hat{\mathcal{T}}_b) \\ \Phi(\hat{w}_h) \cdot n_{\omega_b} |_{F} = r_F \\ [\Phi(\hat{w}_h) \cdot n_{\omega_b}]_F = r_F \\ \nabla_b \cdot \Phi(\hat{w}_h)|_{T} = r_T \\ \nabla \cdot \Phi(\hat{w}_h)|_{\partial T} = r_T}} \|\Phi(\hat{w}_h)\|_{\omega_b}. \quad (87) $$

With $n_{\hat{\omega}_b}$ denoting the outer normal vector on $\partial\hat{\omega}_b$, elementary analysis, cf. [EG21, Lemma 9.11], provides the relation

$$ (DF^{\top} \circ F^{-1}) n_{\omega_b} = n_{\hat{\omega}_b} \circ F^{-1} | (DF^{\top} \circ F^{-1}) n_{\omega_b} | $$

and thus

$$ \Phi(\hat{w}_h) \cdot n_{\omega_b} = \left( (\det(DF)^{-1} \hat{w}_h^{\top} (DF)^{\top}) \circ F^{-1} \right) n_{\omega_b} \\ = \left( \hat{w}_h \cdot n_{\hat{\omega}_b} \frac{|(DF^{\top} \circ F^{-1}) n_{\omega_b}|}{\det(DF)} \right) \circ F^{-1}. $$

Hence, with $\hat{F} := F^{-1}(F)$ and $\hat{\tau}_h := \Phi^{-1}(\tau_h)$, the first equation in the last minimum of (87) is equivalent to $\hat{w}_h \cdot n_{\hat{\omega}_b}|_{\hat{F}} = \hat{\tau}_h \cdot n_{\hat{\omega}_b}|_{\hat{F}} =: r_{\hat{F}}$, and the second equation is equivalent to $[\hat{w}_h \cdot n_{\hat{\omega}_b}]|_{\hat{F}} = [\hat{\tau}_h \cdot n_{\hat{\omega}_b}]|_{\hat{F}} =: r_{\hat{T}}$. The identity (86) shows that the third equation is equivalent to $\nabla_b \cdot \hat{w}_h|_{\hat{T}} = \tilde{\Phi}^{-1}(g_h + \nabla_b \cdot \tau_h) =: r_{\hat{T}}$. As $\|\Phi(\hat{w}_h)\|_{\omega_b} \lesssim \|\hat{w}_h\|_{\hat{\omega}_b}$ (with a hidden constant depending only on $\|DF\|_{\infty, \hat{\Omega}}$), we can formulate (87) in the parameter domain

$$ \min_{\substack{v_h \in V_h^{a,b} \\ \nabla \cdot v_h = g_h}} \|v_h + \tau_h\|_{\omega_b} \lesssim \min_{\substack{\hat{w}_h \in RT^{\tilde{p}}(\hat{\mathcal{T}}_b) \\ \hat{w}_h \cdot n_{\hat{\omega}_b}|_{\hat{F}} = r_{\hat{F}} \\ [\hat{w}_h \cdot n_{\hat{\omega}_b}]_{\hat{F}} = r_{\hat{F}} \\ \nabla_b \cdot \hat{w}_h|_{\hat{T}} = r_{\hat{T}} \\ \nabla \cdot \hat{w}_h|_{\partial \hat{T}} = r_{\hat{T}}}} \|\hat{w}_h\|_{\hat{\omega}_b}. \quad (88) $$
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Oscillation terms osc<sub>h</sub><sup>rel</sup> from (77b) corresponding to the problem of Section 9 with adaptive mesh refinement, polynomial degrees p ∈ {1, ..., 5}, multiplicities m ∈ {1, p}.</figcaption></figure>

Finally, an application of [BPS09, Theorems 5 and 7] in two space dimensions or a similar
procedure relying on [CDD08] in three space dimensions, see also [EV20, Theorem 2.5
and Corollary 3.3] building on [CM10, DGS12], yields that

$$
\begin{array}{ll}
\min_{\substack{\hat{w}_h \in RT^P(\hat{T}_b) \\ \hat{w}_h \cdot n_{\hat{\omega}_b}|_{\hat{F}} = r_{\hat{F}} \\ [\hat{w}_h \cdot n_{\hat{\omega}_b}]_{\hat{F}} = r_{\hat{F}} \\ \nabla_b \cdot \hat{w}_h|_{\hat{T}} = r_{\hat{T}} \\ \forall \hat{F} \in \hat{\mathcal{F}}^N_b \\ \forall \hat{F} \in \hat{\mathcal{F}}^{\text{int}}_b \\ \forall \hat{T} \in \hat{T}_b}} & \|\hat{w}_h\|_{\hat{\omega}_b} \le C_{\text{st}} \\
\min_{\hat{w} \in \mathcal{H}(\text{div}, \hat{T}_b)} & \|\hat{w}\|_{\hat{\omega}_b}, \\
\hat{w} \cdot n_{\hat{\omega}_b}|_{\hat{F}} = r_{\hat{F}} & \forall \hat{F} \in \hat{\mathcal{F}}^N_b \\
[\hat{w} \cdot n_{\hat{\omega}_b}]_{\hat{F}} = r_{\hat{F}} & \forall \hat{F} \in \hat{\mathcal{F}}^{\text{int}}_b \\
\nabla_b \cdot \hat{w}|_{\hat{T}} = r_{\hat{T}} & \forall \hat{T} \in \hat{T}_b
\end{array}
$$

for a generic constant $C_{\text{st}}$ only depending on the shapes of the elements in $\hat{T}_{\text{uni}(0)}$. Here, $\mathbf{H}(\text{div}, \hat{T}_b)$ is the space of $\hat{T}_b$-piecewise $\mathbf{H}(\text{div})$-functions on $\hat{\omega}_b$, with the normal component understood as in [EV20, equation (2.6)]. The required compatibility condition

$$
\sum_{\hat{T} \in \hat{T}_b} (r_{\hat{T}}, 1)_{\hat{T}} - \sum_{\hat{F} \in \hat{\mathcal{F}}^N_b \cup \hat{\mathcal{F}}^{\text{int}}_b} (r_{\hat{F}}, 1)_{\hat{F}} = 0
$$

in case of $\psi_a \psi_b \in H_0^1(\Omega)$ follows from the assumption that $g_h$, and thus $\tilde{\Phi}^{-1}(g_h)$, has
integral mean zero then. Finally,

$$
\begin{array}{ll}
\min_{\substack{\hat{w} \in \mathcal{H}(\text{div}, \hat{T}_b) \\ \hat{w} \cdot n_{\hat{\omega}_b}|_{\hat{F}} = r_{\hat{F}} \\ [\hat{w} \cdot n_{\hat{\omega}_b}]_{\hat{F}} = r_{\hat{F}} \\ \nabla_b \cdot \hat{w}|_{\hat{T}} = r_{\hat{T}} \\ \forall \hat{F} \in \hat{\mathcal{F}}^N_b \\ \forall \hat{F} \in \hat{\mathcal{F}}^{\text{int}}_b \\ \forall \hat{T} \in \hat{T}_b}} & \|\hat{w}\|_{\hat{\omega}_b} \lesssim \min_{\substack{v \in H_0(\text{div}, \omega_b) \\ \nabla \cdot v = g_h}} \|v + \tau_h\|_{\omega_b} \\
\end{array}
$$

with a hidden constant depending only on $\|(DF)^{-1}\|_{\infty, \hat{\Omega}}$, as in (87) and (88). This
concludes the proof.

REFERENCES

[AO00] Mark Ainsworth and J. Tinsley Oden. *A posteriori error estimation in finite element analysis*. Pure and Applied Mathematics (New York). Wiley-Interscience [John Wiley & Sons], New York, 2000.

[BBdVC$^{+06}$] Yuri Bazilevs, Lourenco Beirão da Veiga, J. Austin Cottrell, Thomas J. R. Hughes, and Giancarlo Sangalli. Isogeometric analysis: approximation, stability and error estimates for *h*-refined meshes. *Math. Models Methods Appl. Sci.*, 16(07):1031–1090, 2006.

[BBF13] Daniele Boffi, Franco Brezzi, and Michel Fortin. *Mixed finite element methods and applications*. Springer, Heidelberg, 2013.

[BdVBSV13] Lourenco Beirão da Veiga, Annalisa Buffa, Giancarlo Sangalli, and Rafael Vázquez. Analysis-suitable T-splines of arbitrary degree: definition, linear independence and approximation properties. *Math. Models Methods Appl. Sci.*, 23(11):1979–2003, 2013.

[BdVBSV14] Lourenço Beirão da Veiga, Annalisa Buffa, Giancarlo Sangalli, and Rafael Vázquez. Mathematical analysis of variational isogeometric methods. *Acta Numer.*, 23:157–287, 2014.

[BG16a] Annalisa Buffa and Eduardo M. Garau. Refinable spaces and local approximation estimates for hierarchical splines. *IMA J. Numer. Anal.*, 37(3):1125–1149, 2016.

[BG16b] Annalisa Buffa and Carlotta Giannelli. Adaptive isogeometric methods with hierarchical splines: error estimator and convergence. *Math. Models Methods Appl. Sci.*, 26(01):1–25, 2016.

[BG17] Annalisa Buffa and Carlotta Giannelli. Adaptive isogeometric methods with hierarchical splines: Optimality and convergence rates. *Math. Models Methods Appl. Sci.*, 27(14):2781–2802, 2017.

[BG18] Annalisa Buffa and Eduardo M. Garau. *A posteriori error estimators for hierarchical B-spline discretizations*. *Math. Models Methods Appl. Sci.*, 28(8):1453–1480, 2018.

[BGG$^{+22}$] Annalisa Buffa, Gregor Gantner, Carlotta Giannelli, Dirk Praetorius, and Rafael Vázquez. Mathematical foundations of adaptive isogeometric analysis. *Arch. Comput. Methods Eng.*, 29:4479–4555, 2022.

[BGV18] Cesare Bracco, Carlotta Giannelli, and Rafael Vázquez. Refinement algorithms for adaptive isogeometric methods with hierarchical splines. *Axioms*, 7(3):43, 2018.
[BMV20] Jan Blechta, Josef Málek, and Martin Vohralík. Localization of the $W^{-1,q}$ norm for local a posteriori efficiency. *IMA J. Numer. Anal.*, 40(2):914–950, 2020.

[BPS09] Dietrich Braess, Veronika Pillwein, and Joachim Schöberl. Equilibrated residual error estimates are $p$-robust. *Comput. Methods Appl. Mech. Engrg.*, 198(13-14):1189–1197, 2009.

[BS08] Dietrich Braess and Joachim Schöberl. Equilibrated residual error estimator for edge elements. *Math. Comp.*, 77(262):651–672, 2008.

[CDD08] M. Costabel, M. Dauge, and L. Demkowicz. Polynomial extension operators for $H^1$, $H(\text{curl})$ and $H(\text{div})$-spaces on a cube. *Math. Comp.*, 77(264):1967–1999, 2008.

[CF00] Carsten Carstensen and Stefan A. Funken. Fully reliable localized error control in the FEM. *SIAM J. Sci. Comput.*, 21(4):1465–1484, 1999/00.

[CHB09] J. Austin Cottrell, Thomas J. R. Hughes, and Yuri Bazilevs. *Isogeometric analysis: toward integration of CAD and FEA*. John Wiley & Sons, New York, 2009.

[CM10] Martin Costabel and Alan McIntosh. On Bogovskiǐ and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains. *Math. Z.*, 265(2):297–320, 2010.

[dB86] Carl de Boor. B(asic)-spline basics. Technical report, Mathematics Research Center, University of Wisconsin-Madison, 1986.

[dB01] Carl de Boor. *A practical guide to splines*. Springer, New York, 2001.

[DEV16] Vít Dolejší, Alexandre Ern, and Martin Vohralík. hp-adaptation driven by polynomial-degree-robust a posteriori error estimates for elliptic problems. *SIAM J. Sci. Comput.*, 38(5):A3220–A3246, 2016.

[DGS12] Leszek Demkowicz, Jayadeep Gopalakrishnan, and Joachim Schöberl. Polynomial extension operators. Part III. *Math. Comp.*, 81(279):1289–1326, 2012.

[DLP13] Tor Dokken, Tom Lyche, and Kjell F. Pettersen. Polynomial splines over locally refined box-partitions. *Comput. Aided Geom. Design*, 30(3):331–356, 2013.

[DM99] Philippe Destuynder and Brigitte Métivet. Explicit error bounds in a conforming finite element method. *Math. Comp.*, 68(228):1379–1396, 1999.

[EG21] Alexandre Ern and Jean-Luc Guermond. *Finite Elements I. Approximation and Interpolation*. Springer, Cham, 2021.

[ESV17a] Alexandre Ern, Iain Smears, and Martin Vohralík. Discrete $p$-robust $H(\text{div})$-liftings and a posteriori estimates for elliptic problems with $H^{-1}$ source terms. *Calcolo*, 54(3):1009–1025, 2017.

[ESV17b] Alexandre Ern, Iain Smears, and Martin Vohralík. Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems. *SIAM J. Numer. Anal.*, 55(6):2811–2834, 2017.

[ET99] Ivar Ekeland and Roger Témam. *Convex analysis and variational problems*. SIAM, Philadelphia, 1999.

[EV15] Alexandre Ern and Martin 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(2):1058–1081, 2015.

[EV20] Alexandre Ern and Martin Vohralík. Stable broken $H^1$ and $H(\text{div})$ polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions. *Math. Comp.*, 89(322):551–594, 2020.

[GHP17] Gregor Gantner, Daniel Haberlik, and Dirk Praetorius. Adaptive IGAFEM with optimal convergence rates: Hierarchical B-splines. *Math. Models Methods Appl. Sci.*, 27(14):2631–2674, 2017.

[GJS12] Carlotta Giannelli, Bert Jüttler, and Hendrik Speleers. THB-splines: The truncated basis for hierarchical splines. *Comput. Aided Geom. Design*, 29(7):485–498, 2012.

[GP20] Gregor Gantner and Dirk Praetorius. Adaptive IGAFEM with optimal convergence rates: T-splines. *Comput. Aided Geom. Design*, 81:101906, 2020.

[HCB05] Thomas J. R. Hughes, J. Austin Cottrell, and Yuri Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. *Comput. Methods Appl. Mech. Engrg.*, 194(39):4135–4195, 2005.

[HKMP17] Paul Hennig, Markus Kästner, Philipp Morgenstern, and Daniel Peterseim. Adaptive mesh refinement strategies in isogeometric analysis—A computational comparison. *Comput. Methods Appl. Mech. Engrg.*, 316:424–448, 2017.

[JKD14] Kjetil A. Johannessen, Trond Kvamsdal, and Tor Dokken. Isogeometric analysis using LR B-splines. *Comput. Methods Appl. Mech. Engrg.*, 269:471–514, 2014.
[JRK15] Kjetil A. Johannessen, Filippo Remonato, and Trond Kvamsdal. On the similarities and differences between classical hierarchical, truncated hierarchical and LR B-splines. *Comput. Methods Appl. Mech. Engrg.*, 291:64–101, 2015.

[KT15a] Stefan K. Kleiss and Satyendra K. Tomar. Guaranteed and sharp a posteriori error estimates in isogeometric analysis. *Comput. Math. Appl.*, 70(3):167–190, 2015.

[KT15b] Stefan K. Kleiss and Satyendra K. Tomar. Two-sided robust and sharp a posteriori error estimates in isogeometric discretization of elliptic problems. In *Isogeometric analysis and applications 2014*, volume 107 of *Lect. Notes Comput. Sci. Eng.*, pages 231–246. Springer, 2015.

[LL83] P. Ladevèze and D. Leguillon. Error estimate procedure in the finite element method and applications. *SIAM J. Numer. Anal.*, 20(3):485–509, 1983.

[Mat18] S. Matculevich. Functional approach to the error control in adaptive IgA schemes for elliptic boundary value problems. *J. Comput. Appl. Math.*, 344:394–423, 2018.

[PS47] William Prager and John L. Synge. Approximations in elasticity based on the concept of function space. *Quart. Appl. Math.*, 5:241–269, 1947.

[PV22] Jan Papež and Martin Vohralík. Inexpensive guaranteed and efficient upper bounds on the algebraic error in finite element discretizations. *Numer. Algorithms*, 89:371–407, 2022.

[Rep99] Sergey I. Repin. A posteriori error estimates for approximate solutions to variational problems with strongly convex functionals. *J. Math. Sci.*, 97(4):4311–4328, 1999.

[Rep00a] Sergey I. Repin. A posteriori error estimation for nonlinear variational problems by duality theory. *J. Math. Sci.*, 99(1):927–935, 2000.

[Rep00b] Sergey I. Repin. A posteriori error estimation for variational problems with uniformly convex functionals. *Math. Comp.*, 69(230):481–500, 2000.

[Rep08] Sergey Repin. *A posteriori estimates for partial differential equations*. De Gruyter, Berlin, 2008.

[Sch07] Larry Schumaker. *Spline functions: Basic theory*. Cambridge University Press, Cambridge, 2007.

[SLSH12] Michael A. Scott, Xin Li, Thomas W. Sederberg, and Thomas J. R. Hughes. Local refinement of analysis-suitable T-splines. *Comput. Methods Appl. Mech. Engrg.*, 213:206–222, 2012.

[SM16] Hendrik Speleers and Carla Manni. Effortless quasi-interpolation in hierarchical spaces. *Numer. Math.*, 132(1):155–184, 2016.

[SMPZ08] Joachim Schöberl, Jens M. Melenk, Clemens Pechstein, and Sabine Zaglmayr. Additive Schwarz preconditioning for *p*-version triangular and tetrahedral finite elements. *IMA J. Numer. Anal.*, 28(1):1–24, 2008.

[SZBN03] Thomas W. Sederberg, Jianmin Zheng, Almaz Bakenov, and Ahmad Nasri. T-splines and T-NURCCs. *ACM Trans. Graph.*, 22(3):477–484, 2003.

[TCHM19] Hoang Phuong Thai, Ludovic Chamoin, and Cuong Ha-Minh. *A posteriori error estimation for isogeometric analysis using the concept of constitutive relation error*. *Comput. Methods Appl. Mech. Engrg.*, 355:1062–1096, 2019.

[Ver13] Rüdiger Verfürth. *A posteriori error estimation techniques for finite element methods*. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, 2013.

[VGJS11] Anh-Vu Vuong, Carlotta Giannelli, Bert Jüttler, and Bernd Simeon. A hierarchical approach to adaptive local refinement in isogeometric analysis. *Comput. Methods Appl. Mech. Engrg.*, 200(49):3554–3567, 2011.

[Voh05] Martin Vohralík. On the discrete Poincaré–Friedrichs inequalities for nonconforming approximations of the Sobolev space $H^1$. *Numer. Funct. Anal. Optim.*, 26(7-8):925–952, 2005.

[VV12] Andreas Veeser and Rüdiger Verfürth. Poincaré constants for finite element stars. *IMA J. Numer. Anal.*, 32(1):30–47, 2012.

INSTITUTE OF ANALYSIS AND SCIENTIFIC COMPUTING, TU WIEN, WIEDNER HAUPТSTRASSE 8-10,
1040 VIENNA, AUSTRIA (GREGOR.GANTNER@ASC.TUWIEN.AC.AT). CURRENT ADDRESS: INSTITUTE FOR
NUMERICAL SIMULATION, UNIVERSITY OF BONN, FRIEDRICH-HIRZEBRUCH-ALLEE 7, 53115 BONN,
GERMANY
INRIA, 2 RUE SIMONE IFF, 75589 PARIS, FRANCE & CERMICS, ECOLE DES PONTS, 77455 MARNE-LA-VALLÉE, FRANCE (MARTIN.VOHRALIK@INRIA.FR).