# An equilibrated estimator for mixed finite element discretizations of the curl-curl problem

Théophile Chaumont-Frelet

► To cite this version:

Théophile Chaumont-Frelet. An equilibrated estimator for mixed finite element discretizations of the curl-curl problem. IMA Journal of Numerical Analysis, 2024, 10.1093/imanum/drae007. hal-04177080

HAL Id: hal-04177080

https://inria.hal.science/hal-04177080

Submitted on 3 Aug 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
# AN EQUILIBRATED ESTIMATOR FOR MIXED FINITE ELEMENT DISCRETIZATIONS OF THE CURL-CURL PROBLEM

T. CHAUMONT-FRELET†

**ABSTRACT.** We propose a new a posteriori error estimator for mixed finite element discretizations of the curl-curl problem. This estimator relies on a Prager-Synge inequality, and therefore leads to fully guaranteed constant-free upper bounds on the error. The estimator is also locally efficient and polynomial-degree-robust. The construction is based on patch-wise divergence-constrained minimization problems, leading to a cheap embarrassingly parallel algorithm. Crucially, the estimator operates without any assumption on the topology of the domain, and unconventional arguments are required to establish the reliability estimate. Numerical examples illustrate the key theoretical results, and suggest that the estimator is suited for mesh adaptivity purposes.

**KEY WORDS.** a posteriori error estimate; electromagnetics finite element method; high order method; potential reconstruction; Prager-Synge

## 1. INTRODUCTION

This work develops an equilibrated a posteriori error estimator for mixed finite element discretizations of the curl-curl problem. The curl-curl equation is the prototypical elliptic PDE in $H(\text{curl})$, and constitutes the basic model problem for magnetostatics. In contrast to the $H(\text{curl})$ setting, the development of equilibrated estimators in $H^1$ with application to, e.g., electrostatics, is much more advanced. This introduction reviews the key concepts of equilibrated estimators in $H^1$, highlights the challenges that arise in $H(\text{curl})$, presents the construction of the new estimator and summarizes the key results of the present work.

**1.1. Equilibrated estimators in $H^1$.** A density of static charges $\rho : \Omega \to \mathbb{R}$ generates an electric field $E : \Omega \to \mathbb{R}^3$. The electric field is linked to the electric displacement $D = \varepsilon E$, where $\varepsilon$ is the electric permittivity tensor, and Gauss' law ensures that $\nabla \cdot D = \rho$. Assuming a steady state, Faraday's law stipulates the existence an electric potential $\phi : \Omega \to \mathbb{R}$ such that $E = \nabla \phi$. This leads to the electrostatic problem of finding $\phi$ such that

$$ (1.1) \qquad \left\{ \begin{array}{ll} -\nabla \cdot (\varepsilon \nabla \phi) & = \rho & \text{in } \Omega, \\ \phi & = 0 & \text{on } \Gamma_T, \\ \nabla \phi \cdot \mathbf{n} & = 0 & \text{on } \Gamma_N, \end{array} \right. $$

where the partition $\{\Gamma_T, \Gamma_N\}$ of the boundary $\partial\Omega$ depends on the properties of the materials surrounding $\Omega$.

<footnote>†Inria Université Côte d'Azur, LJAD, CNRS</footnote>
If $E_h$ is any square-integrable vector field, then the following generalized Prager-Synge inequality holds true

$$ (1.2) \qquad \|\boldsymbol{E} - \boldsymbol{E}_h\|_{\boldsymbol{\varepsilon}, \Omega}^2 \le \min_{\tilde{\phi} \in H_{\Gamma_T}^1(\Omega)} \|\boldsymbol{E}_h - \nabla \tilde{\phi}\|_{\boldsymbol{\varepsilon}, \Omega}^2 + \min_{\tilde{D} \in H_{\Gamma_N}^{(\mathrm{div}, \Omega)}} \|\boldsymbol{E}_h - \boldsymbol{\varepsilon}^{-1} \tilde{D}\|_{\boldsymbol{\varepsilon}, \Omega}^2, $$

as shown for instance in [25, Theorem 3.3]. The name Prager-Synge is after the seminal work [41], where Prager and Synge proved (1.2) in the context of elastostatics for fields $E_h = \nabla \phi_h$, $\phi_h \in H_{\Gamma_T}^1(\Omega)$, hence with the first term vanishing. It is noteworthy that the two terms of the right-hand of (1.2) respectively quantify the inability of $E_h$ to satisfy Gauss' and Faraday's laws.

When numerically solving (1.1) with a conforming finite element method, one of the two terms in the right-hand side of (1.2) automatically vanishes by construction. The remaining term can then be bounded using an error estimator, leading to guaranteed error bounds [25].

Specifically, when Lagrange elements are used to approximate (1.1) in primal form, a conforming electric potential $\phi_h \in H_{\Gamma_T}^1(\Omega)$ is computed, and the resulting electric field approximation is $E_h := \nabla \phi_h$. The first term in the right-hand side of (1.2) then vanishes, and the second term is controlled using an equilibrated flux reconstruction whereby a field $D_h$ such that $\nabla \cdot D_h = \rho$ is computed [3, 19, 25, 35, 36, 39].

On the other hand, when employing Raviart-Thomas elements for the mixed form of (1.1), an electric displacement $D_h$ satisfying $\nabla \cdot D_h = \rho$ is immediately computed. In this case, the second term in the right-hand side of (1.2) vanishes. The first term is then estimated by building a field $\phi_h$ such that $\boldsymbol{\varepsilon}D_h - \nabla \phi_h$ is small: this process is called a potential reconstruction [2, 25, 26, 45].

In practice, for both primal and mixed forms, the identity $\nabla \cdot D_h = \rho$ only holds point-wise if $\rho$ is piecewise polynomial on the finite element mesh. In the general case, an additional term corresponding to the approximation of $\rho$ is added to (1.2). This term is fully computable, and when $\rho$ is piecewise smooth, it is actually of higher-order and is usually called a “data oscillation” term.

A chief question is then: How such equilibrated flux and potential should be constructed? For primal discretizations, several approaches have been proposed in the past for the flux reconstruction [3, 35, 36, 39], and in this work, we will follow the technique introduced in [19]. Specifically, instead of solving the global minimization problem

$$ \min_{\substack{D_h \in H_{\Gamma_N}^{(\mathrm{div}, \Omega) \cap \mathrm{RT}_p}(\mathcal{T}_h) \\ \nabla \cdot D_h = \rho}} \|\boldsymbol{E}_h - \boldsymbol{\varepsilon}^{-1} \boldsymbol{D}_h\|_{\boldsymbol{\varepsilon}, \Omega}^2, $$

with Raviart-Thomas elements, where $p$ is the polynomial degree of the Lagrange finite elements, the localized version

$$ (1.3) \qquad D_h^a := \min_{\substack{v_h \in H_0^{(\mathrm{div}, \boldsymbol{\omega}^a) \cap \mathrm{RT}_{p+1}}(\mathcal{T}_h^a) \\ \nabla \cdot v_h = \psi^a \rho - \nabla \psi^a \cdot \nabla \phi_h}} \|\boldsymbol{E}_h - \boldsymbol{\varepsilon}^{-1} v_h\|_{\boldsymbol{\varepsilon}, \omega^a}^2, $$

is considered for each vertex $a$ of the mesh, where $\psi^a$ is the associated hat function (these notations are rigorously introduced in Section 3.4 below). Crucially, the compatibility condition implied by Stokes' formula in (1.3) is satisfied due to Galerkin orthogonality. Indeed $\phi_h$ is the Galerkin finite element solution and $\psi^a$ is a lowest-order Lagrange finite element
function. One easily checks that the corresponding field

$$D_h := \sum_{a \in V_h} D_h^a$$

is an equilibrated flux, so that (1.2) leads to a constant-free error estimate.

Similar localization strategies involving the hat functions of the finite element mesh have also been introduced for the potential reconstruction, and we refer the reader to [25] for more details.

Another key aspect of equilibrated estimators built from local minimization problems is that they are “polynomial-degree-robust” [11, 25, 26]. It means that the estimator is also a lower bound to the error up to a generic constant, and that this constant does not depend on the polynomial degree $p$ of the finite element space. Therefore, equilibrated estimators are (at least theoretically) more attractive than traditional residual-based estimators for high-order finite element methods and/or hp-adaptive algorithms.

## 1.2. Prager-Synge estimates in $H(\text{curl})$.
In this work, we are rather interested in magnetostatics [32]. In this case, a magnetic field $H: \Omega \to \mathbb{R}^3$ is generated by a static current $J: \Omega \to \mathbb{R}^3$ with $\nabla \cdot J = 0$. The magnetic field is linked to the induction field $B$ through the constitutive relation $B = \mu H$, where $\mu$ is the magnetic permeability tensor. Ampère's law and Gauss' law for magnetism then respectively demand that $\nabla \times H = J$ and $B = \nabla \times A$ for some vector-potential $A: \Omega \to \mathbb{R}^3$ with $\nabla \cdot A = 0$. The curl-curl equation is then obtained by reformulating the problem in terms of $A$ as follows:

$$ (1.4) \qquad \left\{ \begin{array}{lll} \nabla \times (\mu^{-1} \nabla \times A) & = & J \quad \text{in } \Omega, \\ \nabla \cdot A & = & 0 \quad \text{in } \Omega, \\ (\nabla \times A) \times n & = & o \quad \text{on } \Gamma_T, \\ A \times n & = & o \quad \text{on } \Gamma_N. \end{array} \right. $$

Depending on the topology of $\Omega$, $\Gamma_T$ and $\Gamma_N$ a finite number of additional gauge conditions may be imposed to uniquely define $A$. This is detailed in Section 2.3 below.

The magnetostatic problem in (1.4) has a structure similar to the electrostatic problem. As a result, it also enjoys a generalized Prager-Synge inequality. Specifically, if $H_h$ is any square-integrable vector field, we have

$$ (1.5) \quad \|H - H_h\|_{\mu, \Omega}^2 \le \min_{\tilde{A} \in H_{\Gamma_N}(\text{curl}, \Omega)} \|H_h - \mu^{-1} \nabla \times \tilde{A}\|_{\mu, \Omega}^2 + \min_{\tilde{H} \in H_{\Gamma_T}(\text{curl}, \Omega)} \|H_h - \tilde{H}\|_{\mu, \Omega}^2, $$
$$ \nabla \times \tilde{H} = J $$

where, as a counterpart to (1.2), the two terms in the right-hand side respectively measure the failure of Ampère's and Gauss' laws.

The curl-curl equation (1.4) can be immediately discretized in primal form. In this case, an approximation $A_h$ of $A$ is computed with Nédélec elements, and the divergence constraint is imposed with (the gradient of) Lagrange elements as Lagrange multipliers. As for primal discretizations of (1.1), in this case, we set $H_h := \mu^{-1}A_h$, and the first term in the right-hand side of (1.5) vanishes. The error is then estimated by reconstructing a field $\tilde{H}$ such that $\nabla \times \tilde{H} = J$, up to data oscillations. These curl-constrained equilibrated reconstructions turn out to be substantially more complex than in the $H^1$ setting, and general polynomial degrees $p$ have only been handled very recently [12, 17, 28, 29].

Here, we focus on the discretization of the mixed form of (1.4) which is not currently covered in the literature. Namely, $H$ is immediately discretized with Nédélec elements, so
that $\nabla \times H_h = J$ and the second term in the right-hand side of (1.5) vanishes up to data oscillation. A potential reconstruction is thus required to control the first term.

1.3. **A novel error estimator.** To analyze the efficiency of estimators based on a potential reconstruction in $H^1$, one of the key tools is a broken Poincaré inequality for scalar functions with jumps of vanishing mean value, see e.g. [25, Lemma 3.13 and Theorem 3.17] and [26, Corollary 4.1] as well as [13, 34, 44]. However, to the best of the author's knowledge, such an inequality is not available in the $H(\text{curl})$ context.

In this work, we therefore follow an alternative strategy. Specifically, instead of constructing a potential $A_h$ with Nédélec elements and then take its curl in (1.5), we immediately construct a Raviart-Thomas function $B_h$ in the range of the curl operator. $B_h$ is then used in place of $\nabla \times A_h$ in (1.5).

If $\Omega$ is homotopy equivalent to a ball, and either $\Gamma_T$ or $\Gamma_N = \emptyset$, then it is in fact sufficient to build a divergence-free field $B_h \in H_{\Gamma_N}(\text{div}, \Omega)$. This observation leads us to consider a divergence-constrained minimization problem of the form

$$ \min_{ \tilde{B} \in H_{\Gamma_N}(\text{div}, \Omega) \atop \nabla \cdot \tilde{B} = 0 } \| H_h - \tilde{\mu}^{-1} \tilde{B} \|_{\mu, \Omega}^2, $$

which is very similar to the second term in the right-hand side of (1.2). In fact, the localization technique introduced in [19] and further analyzed in [11, 25] can be accommodated, leading to the localized divergence-constrained problems

$$ (1.6a) \qquad B_h^a := \min_{ \substack{ v_h \in H_0(\text{div}, \omega^a) \cap \mathcal{R} T_{p+2}(\mathcal{T}_h^a) \\ \nabla \cdot v_h = \nabla \psi^a \cdot H_h } } \| \psi^a H_h - \mu^{-1} v_h \|_{\mu, \omega^a}^2. $$

These local contributions are then assembled

$$ (1.6b) \qquad B_h := \sum_{a \in \mathcal{V}_h} B_h^a $$

leading to the estimator

$$ (1.6c) \qquad \eta_K := \| H_h - \mu^{-1} B_h \|_{\mu, K}, $$

for all $K \in \mathcal{T}_h$. It is clear that $\nabla \cdot B_h = 0$, so that $B_h$ is an admissible field to plug in (4.1) when the cohomology of $\Omega$ is trivial, thus leading to guaranteed error bounds.

The interesting (and perhaps surprising) part of this work, is that this simple procedure is still valid without any assumption on the topology of $\Omega$, $\Gamma_T$ or $\Gamma_N$. Specifically, not only does the localization procedure produce a divergence-free field $B_h$, but it also guarantees that this field is in the range of the curl operator. This result is obtained by checking that the flux of $B_h$ vanishes through all closed surfaces. This is a striking result, since it is a global property that is not actively enforced in the construction based on local problems.

1.4. **Main results.** The key theoretical results of this work may be summarized as follows. The local problems (1.6a) are well-posed, and the vector field $B_h$ constructed following (1.6b) always sits in the range of the curl operator, without any assumption on the topology. The corresponding error estimator in (1.6c) is reliable up to data oscillations and locally efficient:

$$ (1.7) \qquad \| H - H_h \|_{\mu, \Omega}^2 \le \sum_{K \in \mathcal{T}_h} \eta_K^2 + \text{osc}^2, \qquad \eta_K \lesssim \| H - H_h \|_{\mu, \tilde{K}}, $$
for all $K \in \mathcal{T}_h$. In addition, the hidden constant in the lower bound does not depend on the polynomial degree $p$, so that the estimator is polynomial-degree-robust. These estimates are rigorously stated and established in Theorem 4.1 and Corollary 5.2 below.

We also illustrate the theoretical findings with numerical examples. Crucially, we observe that (1.7) indeed sharply holds numerically. We also investigate the accuracy of the reconstructed field $B_h$ as opposed to $\mu H_h$ to approximate the magnetic induction $B$. Our observation is that although $B_h$ does not exhibit super convergence properties, it is 2 to 5 times more accurate than $\mu H_h$ in our examples. We finally employ the estimator for adaptive mesh refinements, and observe optimal convergence rates.

**1.5. Outline.** The remainder of this work is organized as follows. Section 2 properly states our model problem, whereas Section 3 introduces the discretization setting. Sections 4 and 5 respectively contain the reliability and efficiency proofs. We present numerical examples in Section 6, before drawing concluding remarks in Section 7.

## 2. CONTINUOUS SETTING

**2.1. Domain and coefficient.** We consider a polyhedral domain $\Omega \subset \mathbb{R}^3$ with a Lipschitz boundary $\partial\Omega$. $\partial\Omega$ is split into two disjoint, relatively open, and polytopal components $\Gamma_N$ and $\Gamma_T$ with Lipschitz boundaries (the case where these manifolds do not have boundaries is allowed).

We also fix a symmetric matrix-valued coefficient $\mu : \Omega \to \mathbb{R}^{3 \times 3}$. Classically, we require that $\mu$ is uniformly elliptic and bounded, and for the sake of simplicity, $\chi := \mu^{-1}$ will denote the inverse of $\mu$. We further assume that there exists a partition $\mathcal{P}$ of $\Omega$ into non-overlapping polyhedral subsets $P$ such that for each $P \in \mathcal{P}$, $\mu|_P$ and $\chi|_P$ are constant matrices.

The “coefficient contrast” will appear in the efficiency bounds of the proposed estimator. For any open subset $\omega \subset \Omega$, it is defined as

$$ \mathcal{C}_{\mu, \omega} := \operatorname*{ess\,sup}_{x \in \omega} \max_{\substack{\xi \in \mathbb{R}^3 \\ |\xi|=1}} \frac{\mu\xi \cdot \xi}{\operatorname*{ess\,inf}_{x \in \omega} \min_{\substack{\xi \in \mathbb{R}^3 \\ |\xi|=1}} \mu\xi \cdot \xi}. $$

**2.2. Functional spaces.** For an open set $\omega \subset \mathbb{R}^3$, $L^2(\omega)$ is the usual Lebesgue space of real-valued square integrable functions, and $L^2(\omega) := [L^2(\omega)]^3$ contains vector-valued functions [1]. The standard inner-products of both $L^2(\omega)$ and $L^2(\omega)$ are denoted by $(\cdot, \cdot)_\omega$. We employ the notation $\|\cdot\|_\omega$ for the usual norm of $L^2(\omega)$ associated with $(\cdot, \cdot)_\omega$. We will also frequently used the (equivalent) weighted norms $\|\cdot\|_{\mu, \omega}^2 := (\mu \cdot, \cdot)_\omega$ and $\|\cdot\|_{\chi, \omega}^2 := (\chi \cdot, \cdot)_\omega$.

For standard Sobolev spaces, we will use the notations

$$ H^1(\omega) := \{v \in L^2(\omega); \nabla v \in L^2(\omega)\}, $$

$$ H(\text{curl}, \omega) := \{v \in L^2(\omega); \nabla \times v \in L^2(\omega)\}, $$

$$ H(\text{div}, \omega) := \{v \in L^2(\omega); \nabla \cdot v \in L^2(\omega)\}, $$

where $\nabla$, $\nabla \times$ and $\nabla \cdot$ are the weak gradient, curl and divergence operators defined in the sense of distributions, see [1, 31].

If $\gamma \subset \partial\omega$, we respectively denote by $H_\gamma^1(\omega)$, $H_\gamma(\text{curl}, \omega)$ and $H_\gamma(\text{div}, \omega)$ the closure of smooth functions that vanish on $\gamma$ in $H^1(\omega)$, $H(\text{curl}, \omega)$ and $H(\text{div}, \omega)$. Assuming that $\omega$ has a Lipschitz boundary and that $\gamma$ is relatively open, these spaces may be interpreted as containing functions with vanishing trace, tangential trace and normal trace on $\gamma$, see [27].

Finally, we introduce the short-hand notation $\Lambda_\gamma(\omega) := \nabla \times H_\gamma(\text{curl}, \omega)$.
2.3. **Cohomology.** The space $\Lambda_{\Gamma_N}(\Omega)$ will play a key role in the following. It is clear that for all $v \in H_{\Gamma_N}(\text{curl}, \Omega)$, $w := \nabla \times v \in H_{\Gamma_N}(\text{div}, \Omega)$ with $\nabla \cdot w = 0$, so that

$$ \Lambda_{\Gamma_N}(\Omega) \subset H_{\Gamma_N}(\text{div}^0, \Omega) := \{w \in H_{\Gamma_N}(\text{div}, \Omega) \mid \nabla \cdot w = 0\}. $$

In fact, if $\Omega$ is homotopy equivalent to a ball and either $\Gamma_T = \emptyset$ or $\Gamma_N = \emptyset$, then the identity $\Lambda_{\Gamma_N}(\Omega) = H_{\Gamma_N}(\text{div}^0, \Omega)$ holds true.

In the general case, a finite number of linear constraints must be imposed on $w \in H_{\Gamma_N}(\text{div}^0, \Omega)$ to ensure the existence of $v \in H_{\Gamma_N}(\text{curl}, \Omega)$ such that $w = \nabla \times v$. Specifically, following [9, 23, 33], there exists a finite number of (relatively open) oriented surfaces $\Sigma_j \subset \bar{\Omega}$ with unit normal vector $n_{\Sigma_j}$ and $\partial\Sigma_j \subset \bar{\Gamma}_N$, $1 \le j \le N$, such that $w \in \Lambda_{\Gamma_N}(\text{div}, \Omega)$ if and only if $w \in H_{\Gamma_N}(\text{div}^0, \Omega)$ and

$$ (2.1) \qquad \int_{\Sigma_j} w \cdot n_{\Sigma_j} = 0 $$

for $1 \le j \le N$. The integrals in the left-hand side of (2.1) are called the “periods” of $w$. A function $w \in H_{\Gamma_N}(\text{div}^0, \Omega)$ thus admits a vector potential if and only if its periods vanish.

**Remark 2.1 (Practical implementation).** The surfaces $\{\Sigma_j\}_{j=1}^N$ are only needed to develop the theory. They are not used in the implementation and do not explicitly enter the construction of the proposed error estimator.

2.4. **Model problem.** Given a right-hand side $J \in \Lambda_{\Gamma_T}(\Omega)$, our model problem is to find $H \in H_{\Gamma_T}(\text{curl}, \Omega)$ and $A \in \Lambda_{\Gamma_T}(\Omega)$ such that

$$ (2.2) \qquad \left\{ \begin{array}{l} (\mu H, v)_{\Omega} + (A, \nabla \times v)_{\Omega} = 0, \\ (\nabla \times H, w)_{\Omega} = (J, w)_{\Omega}, \end{array} \right. $$

for all $v \in H_{\Gamma_T}(\text{curl}, \Omega)$ and $w \in \Lambda_{\Gamma_T}(\Omega)$. The system of equations in (2.2) is a saddle-point problem. By definition of $\Lambda_{\Gamma_T}(\Omega)$, the so-called “inf-sup” condition is trivially satisfied, so that there exists a unique solution [14].

### 3. DISCRETE SETTING

3.1. **Computational mesh.** Throughout this work, we consider a fixed mesh $\mathcal{T}_h$ of $\Omega$ consisting of (open) tetrahedral elements $K$. The set of all mesh vertices and faces are respectively denoted by $\mathcal{V}_h$ and $\mathcal{F}_h$, and we further split $\mathcal{F}_h$ into the set $\mathcal{F}_h^e$ of exterior faces $F \in \mathcal{F}_h$ such that $F \subset \partial\Omega$ and the set interior faces $\mathcal{F}_h^i := \mathcal{F}_h \setminus \mathcal{F}_h^e$. For each element $K \in \mathcal{T}_h$ and face $F \in \mathcal{F}_h$, $\mathcal{V}(K)$, $\mathcal{V}(F) \subset \mathcal{V}_h$ are the sets of vertices of $K$ and $F$.

We require that the mesh is conforming in the sense of [10, Section 2.1.2] and [24, Definition 6.11]. Specifically, we assume that $\bar{\Omega} = \cup_{K \in \mathcal{T}_h} \bar{K}$, and that if $K_{\pm} \in \mathcal{T}_h$ are two distinct elements, $\bar{K}_+ \cap \bar{K}_-$ is either empty, or a single vertex, edge or face of both $K_+$ and $K_-$. This assumption is standard, and although it rules out hanging nodes, it does not prevent heavily localized refinements with strong mesh grading [8].

We further demand that the mesh fits the coefficient $\mu$, meaning that for all $K \in \mathcal{T}_h$, their exists $P \in \mathcal{P}$ such that $\bar{K} \subset \bar{P}$. In effect, $\mu|_K$ is a constant value for all $K \in \mathcal{T}_h$.

We finally assume that the surfaces $\{\Sigma_j\}_{j=1}^N$ characterizing the space $\Lambda_{\Gamma_N}(\Omega)$ are aligned with the mesh, meaning that for $1 \le j \le N$,

$$ (3.1) \qquad \bar{\Sigma}_j := \bigcup_{F \in \mathcal{F}_h^{\Sigma_j}} F $$
for some $\mathcal{F}_h^{\Sigma_j} \subset \mathcal{F}_h$. Notice that because the particular choice of surfaces $\{\Sigma_j\}_{j=1}^N$ is not important, but only their topological properties, (3.1) is by no means a restrictive assumption in practice [33]. In fact, software packages are available to automatically generate a set of surfaces $\{\Sigma_j\}_{j=1}^N$ aligned with faces, given any conforming mesh as input, see e.g. [20, 40].

3.2. **Mesh size and shape-regularity parameters.** For all $K \in \mathcal{T}_h$, $h_K$ and $\rho_K$ respectively denote the diameters of the smallest ball containing $K$ and the largest ball contained in $\bar{K}$. We employ the notation $h_{\max} := \max_{K \in \mathcal{T}_h} h_K$ for the mesh size. The quantity $\kappa_K := h_K/\rho_K \ge 1$ is called the shape-regularity parameter of $K$, and for $\mathcal{T} \subset \mathcal{T}_h$ we set $\kappa_{\mathcal{T}} := \max_{K \in \mathcal{T}} \kappa_K$.

3.3. **Jumps.** Each face $F \in \mathcal{F}_h$ is equipped with a unit normal vector $n_F$. If $F \in \mathcal{F}_h^i$ the orientation of $n_F$ is arbitrarily fixed, and we assume that $n_F = n$ if $F \in \mathcal{F}_h^e$. If $v: \Omega \to \mathbb{R}$ is a piecewise smooth function, its jump through $F = \bar{K}_+ \cap \bar{K}_- \in \mathcal{F}_h^i$ is defined by

$$[v]_F := v_+|_F(n_+ \cdot n_F) + v_-|_F(n_- \cdot n_F)$$

where $v_{\pm} := v|_{K_{\pm}}$ and $n_{\pm}$ is the unit normal vector of $K_{\pm}$. We then define the jump $[v]_F$ a piecewise smooth vector-valued function $v: \Omega \to \mathbb{R}^3$ using the same formula for each component.

3.4. **Vertex patches.** For all $a \in \mathcal{V}_h$, $\mathcal{T}_h^a \subset \mathcal{T}_h$ collects the elements $K \in \mathcal{T}_h$ such that $a \in \mathcal{V}(K)$. We also denote by $\omega^a$ the corresponding open domain. The hat function $\psi^a: \Omega \to \mathbb{R}$ is the only function affine in each element $K \in \mathcal{T}_h$ such that $\psi^a(a) = 1$ and $\psi^a(b) = 0$ for all $b \in \mathcal{V}_h \setminus \{a\}$. Notice that $\bar{\omega}^a$ then corresponds to the support of the hat function $\psi^a$.

We let $\Gamma^a := \Gamma_N \cup \{\psi^a = 0\}$ and

$$H_0(\text{div}, \omega^a) := \{v \in H(\text{div}, \omega^a) \mid v \cdot n = 0 \text{ on } \Gamma^a\}.$$

We also introduce $L^2_*(\omega^a) := \nabla \cdot H_0(\text{div}, \omega^a)$, so that $L^2_*(\omega^a)$ consists of functions with zero mean value when $\Gamma^a = \partial\omega^a$, and coincide with the whole $L^2(\omega^a)$ otherwise. Finally, we set

$$H^1_*(\omega^a) := \{v \in H^1(\omega^a) \mid v = 0 \text{ on } \partial\omega \setminus \Gamma^a\} \cap L^2_*(\omega^a).$$

The inequality

$$(3.2) \qquad \|\nabla(\psi^a q)\|_{\mu, \omega^a} \le C_{\text{cont}, a} \mathcal{C}_{\mu, \omega^a}^{1/2} \|\nabla q\|_{\mu, \omega^a} \quad \forall q \in H^1_*(\omega^a)$$

holds true with a constant $C_{\text{cont}, a}$ only depending on the shape-regularity parameter of $\mathcal{T}_h^a$, see, e.g., [25, Equation (3.29)].

3.5. **Surfaces and half patches.** Consider a relatively open Lipschitz surface $\Sigma \subset \bar{\Omega}$ with $\partial\Sigma \subset \partial\Omega$ that corresponds to a collection of faces of the mesh $\mathcal{T}_h$. We assume that $\Sigma$ is oriented by a unit normal vector $n_{\Sigma}$. Consider a mesh vertex $a \in \mathcal{V}_h$ such that $a \in \bar{\Sigma}$ and the associated patch $\omega^a$. Then we may uniquely decompose $\omega^a$ into two connected subsets $\omega_{\Sigma, \pm}^a$ with outward unit normal $n_{\omega_{\Sigma, \pm}^a}$ in such way that $n_{\omega_{\Sigma, \pm}^a} = \pm n_{\Sigma}$ on $\Sigma \cap \omega^a$. Notice that one of the two subsets may be empty if $\Sigma \cap \omega^a \subset \partial\Omega$. These two halves of the patch correspond to unions of mesh cells that we denote by $\mathcal{T}_h^{a, \Sigma, \pm}$. We also introduce the notation $\Gamma_{\Sigma, \pm}^a := \Gamma^a \cap \omega_{\Sigma, \pm}^a$.

Finally, we note that if $a \in \bar{\Gamma}_N$ and $\partial\Sigma \subset \bar{\Gamma}_N$, then at least one of the two halves $\omega_{\Sigma, \pm}^a$ of the patch $\omega^a$ is not empty with the boundary $\Gamma_{\Sigma, \pm}$ lying in $\Gamma_N$.
3.6. **Finite element spaces.** If $K \in \mathcal{T}_h$ is an element of the mesh and $q \ge 0$, then $\mathcal{P}_q(K)$ is the set of polynomials mapping $K$ into $\mathbb{R}$ of degree at most $q$, and $\mathcal{P}_q(K) := [\mathcal{P}_q(K)]^3$ collects vector-valued polynomial functions. $\mathcal{N}_q(K) := \mathcal{P}_q(K) \times \mathbf{x} + \mathcal{P}_q(K)$ and $\mathbf{RT}_q(K) := \mathcal{P}_q(K)\mathbf{x} + \mathcal{P}_q(K)$ are then the space of Nédélec [37] and Raviart-Thomas [42] polynomials of order $q$, see also [10, Chapter 2.3] and [24, Chapters 11 and 12].

For a subset of elements $\mathcal{T} \subset \mathcal{T}_h$, $\mathcal{P}_q(\mathcal{T})$ stands for piecewise polynomial functions on $\mathcal{T}$, i.e. $v \in \mathcal{P}_q(\mathcal{T})$ if and only if $v|_K \in \mathcal{P}_q(K)$ for all $K \in \mathcal{T}$. $\mathcal{P}_q(\mathcal{T})$, $\mathcal{N}_q(\mathcal{T})$ and $\mathbf{RT}_q(\mathcal{T})$ are defined similarly.

3.7. **Stable discrete minimization.** Following [26, Theorem 2.3] (see also [16, Proposition 3.1] for general configurations of boundary patches), for all vertices $\mathbf{a} \in \mathcal{V}_h$, there exists a constant $C_{\text{st},\mathbf{a}}$ only depending on the shape-regularity parameter $\kappa_{\mathcal{T}_h^\mathbf{a}}$ of $\mathcal{T}_h^\mathbf{a}$ such that

$$ (3.3) \qquad \min_{\substack{\mathbf{v}_h \in \mathbf{RT}_q(\mathcal{T}_h^\mathbf{a}) \cap \mathcal{H}_0(\text{div}, \omega^\mathbf{a}) \\ \nabla \cdot \mathbf{v}_h = r_h}} \|\mathbf{v}_h - \boldsymbol{\xi}_h\| \le C_{\text{st},\mathbf{a}} \min_{\substack{\mathbf{v} \in \mathcal{H}_0(\text{div}, \omega^\mathbf{a}) \\ \nabla \cdot \mathbf{v} = r_h}} \|\mathbf{v} - \boldsymbol{\xi}_h\| $$

for all polynomial degree $q \ge 0$, $r_h \in \mathcal{P}_q(\mathcal{T}_h^\mathbf{a}) \cap L_\star^2(\omega^\mathbf{a})$ and $\boldsymbol{\xi}_h \in \mathbf{RT}_q(\mathcal{T}_h^\mathbf{a})$. Crucially, both minimizers in (3.3) are uniquely defined, and $C_{\text{st},\mathbf{a}}$ does not depend on $q$.

3.8. **Discrete solution.** Fix $p \ge 0$. We denote by $V_h := \mathcal{N}_p(\mathcal{T}_h) \cap \mathbf{H}_{\Gamma_T}(\text{curl}, \Omega)$ the usual space of Nédélec finite elements and set $\mathbf{\Lambda}_h := \nabla \times V_h$. There exists a unique couple $(\mathbf{H}_h, \mathbf{A}_h) \in V_h \times \mathbf{\Lambda}_h$ such that

$$ (3.4) \qquad \left\{ \begin{array}{l} (\mu \mathbf{H}_h, \mathbf{v}_h)_\Omega + (\mathbf{A}_h, \nabla \times \mathbf{v}_h)_\Omega = 0, \\ (\nabla \times \mathbf{H}_h, \mathbf{w}_h)_\Omega = (\mathbf{J}, \mathbf{w}_h)_\Omega. \end{array} \right. $$

for all $\mathbf{v}_h \in V_h$ and $\mathbf{w}_h \in \mathbf{\Lambda}_h$. Similar to (2.2), the existence and uniqueness of the solution is due to the definition of $\mathbf{\Lambda}_h$ and the results in [14].

**Remark 3.1** (Practical construction of $\mathbf{\Lambda}_h$). Let us consider the set of Raviart-Thomas elements $\mathbf{W}_h := \mathbf{RT}_p(\mathcal{T}_h) \cap \mathbf{H}_{\Gamma_T}(\text{div}, \Omega)$. Then, we have $\mathbf{\Lambda}_h \subset \mathbf{W}_h$, and (at least) two approaches are available to construct $\mathbf{\Lambda}_h$.

First (i), similar to (2.1), it is known that $\mathbf{w}_h \in \mathbf{\Lambda}_h$ if and only if $\mathbf{w}_h \in \mathbf{W}_h$ with $\nabla \cdot \mathbf{w}_h = 0$ and

$$ (3.5) \qquad \int_{\tilde{\Sigma}_j} \mathbf{w}_h \cdot \mathbf{n}_{\tilde{\Sigma}_j} = 0 $$

for $1 \le j \le \tilde{N}$, where each $\tilde{\Sigma}_j \subset \bar{\Omega}$ is a (relatively open) surface with $\partial\tilde{\Sigma}_j \subset \Gamma_T$. We may assume without generality that each $\tilde{\Sigma}_j$ is exactly covered by mesh faces, so that (3.5) may be efficiently checked numerically. As previously mentioned, efficient algorithms are available to automatically construct $\{\tilde{\Sigma}_j\}_{j=1}^{\tilde{N}}$, see e.g. [20, 33, 40]. Then, instead of explicitly constructing a basis of $\mathbf{\Lambda}_h$ we augment the discrete formulation (3.4) with two Lagrange multipliers, namely $q_h \in \nabla \cdot \mathbf{W}_h \subset \mathcal{P}_p(\mathcal{T}_h)$ to impose the divergence constraint, and $\zeta \in \mathbb{R}^{\tilde{N}}$ to enforce (3.5).

As an alternative (ii), it is possible to employ a spanning-tree approach to explicitly construct a basis of $\mathbf{\Lambda}_h$, see e.g. [4, 5, 43].

3.9. **Oscillation term.** For $\phi \in L^2(\Omega)$, consider the following orthogonal projector onto the subset $\mathbf{\Lambda}_h$ of Raviart-Thomas elements

$$ \pi_h^0 \phi := \arg \min_{\mathbf{v}_h \in \mathbf{\Lambda}_h} \|\phi - \mathbf{v}_h\|_{\mathcal{X}, \Omega}. $$
Then, the following quantity will be central for the data oscillation term

$$ \beta_h := \sup_{\substack{\phi \in H_{\Gamma_T}(\text{curl}, \Omega) \cap \Lambda_{\Gamma_N}(\Omega) \\ \|\nabla \times \phi\|_{\chi, \Omega} = 1}} \|\phi - \pi_h^0 \phi\|_{\chi, \Omega}. $$

Notice that $\beta_h$ has dimension (length)$^{-1}$. Actually, when the domain is convex (or smooth) and either $\Gamma_T$ or $\Gamma_N$ is empty, the inclusion $H_{\Gamma_T}(\text{curl}, \Omega) \cap \Lambda_{\Gamma_N} \subset H^1(\Omega)$ holds true [7, Theorem 2.17], and $\beta_h \sim h_{\max}$. A similar property still holds true in more general domains, if the mesh is properly graded in the vicinity of re-entrant corners and edges [8, 38]. In the general case, however, we have $\beta_h \sim \ell_{\Omega}^{1-s} h_{\max}^s$ for some $s > 0$ for which the inclusion $H_{\Gamma_T}(\text{curl}, \Omega) \cap \Lambda_{\Gamma_N} \subset H^s(\Omega)$ holds true [18]. When either $\Gamma_T$ or $\Gamma_N$ is empty, we can even select $s > 1/2$, see [7, Proposition 3.7].

## 4. RELIABILITY

### 4.1. General upper bound.
We start with a general upper-bound given by any field $B \in \Lambda_{\Gamma_N}(\Omega)$. The bound (4.1) is to be compared with the generalized Prager-Synge bound (1.2) for the electrostatic problem.

**Theorem 4.1 (General upper-bound).** The estimate

$$ (4.1) \qquad \|H - H_h\|_{\mu, \Omega}^2 \le \|B - \mu H_h\|_{\chi, \Omega}^2 + \beta_h^2 \|J - \pi_h^0 J\|_{\mu, \Omega}^2 $$

holds true for all $B \in \Lambda_{\Gamma_N}(\Omega)$.

*Proof.* Consider an arbitrary test function $v \in L^2(\Omega)$. The space $\chi\Lambda_{\Gamma_N}(\Omega)$ is closed in $L^2(\Omega)$, so that we may introduce the decomposition $v = v_0 + v^{\perp}$, with $v_0 \in \chi\Lambda_{\Gamma_N}(\Omega)$ and $v^{\perp} \in (\chi\Lambda_{\Gamma_N}(\Omega))^{\perp}$, orthogonal for the $(\mu, \cdot)_{\Omega}$ inner-product. Then, we have

$$ (\mu(H - H_h), v)_{\Omega} = (H - H_h, \mu v_0)_{\Omega} + (\mu(H - H_h), v^{\perp})_{\Omega}. $$

The first equation of the magnetostatic problem (2.2) ensures that $\mu H \in \Lambda_{\Gamma_N}(\Omega)$. Hence, $H \in \chi\Lambda_{\Gamma_N}(\Omega)$, and $(\mu H, v^{\perp})_{\Omega} = 0$. As a result, we have

$$ (\mu(H - H_h), v^{\perp})_{\Omega} = (\mu(\chi B - H_h), v^{\perp})_{\Omega} $$

for all $B \in \Lambda_{\Gamma_N}(\Omega)$.

On the other hand, since $v_0 \in \chi\Lambda_{\Gamma_N}(\Omega)$, there exists $\phi \in H_{\Gamma_N}(\text{curl}, \Omega)$ such that $\mu v_0 = \nabla \times \phi$. In addition, since only $\nabla \times \phi$ (as opposed to $\phi$ itself) intervenes, we may select $\phi$ such that $\phi \in \Lambda_{\Gamma_T}(\Omega)$. Specifically, this may be done by selecting $\psi \in H_{\Gamma_T}(\text{curl}, \Omega)$ such that $\nabla \times \nabla \times \psi = v_0$ in $\Omega$, e.g. by solving (2.2) with $v_0$ instead of $J$ as a right-hand side, and then letting $\phi = \nabla \times \psi$. Then, we have

$$ (H - H_h, \mu v_0)_{\Omega} = (H - H_h, \nabla \times \phi)_{\Omega} = (J - \pi_h^0 J, \phi)_{\Omega}. $$

Since $\pi_h^0$ is an orthogonal projector, we can then write that

$$ (H - H_h, \mu v_0)_{\Omega} = (J - \pi_h^0 J, \phi - \pi_h^0 \phi)_{\Omega} \le \|J - \pi_h^0 J\|_{\mu, \Omega} \|\phi - \pi_h^0 \phi\|_{\chi, \Omega}, $$

and it follows from the definition of $\beta_h$ and the fact that $\nabla \times \phi = \mu v_0$ that

$$ (H - H_h, \mu v_0)_{\Omega} \le \beta_h \|J - \pi_h^0 J\|_{\mu, \Omega} \|\nabla \times \phi\|_{\chi, \Omega} = \beta_h \|J - \pi_h^0 J\|_{\mu, \Omega} \|v_0\|_{\mu, \Omega}. $$
We can then conclude the proof with

$$ 
\begin{aligned}
(\mu(H - H_h), v) &\le \beta_h \|J - \pi_h^0 J\|_{\mu, \Omega} \|v_0\|_{\mu, \Omega} + \|\chi B - H_h\|_{\mu, \Omega} \|v^{\perp}\|_{\mu, \Omega}, \\
&\le (\beta_h^2 \|J - \pi_h^0 J\|_{\mu, \Omega}^2 + \|\chi B - H_h\|_{\mu, \Omega}^2)^{1/2} (\|v_0\|_{\mu, \Omega}^2 + \|v^{\perp}\|_{\mu, \Omega}^2)^{1/2} \\
&= (\beta_h^2 \|J - \pi_h^0 J\|_{\mu, \Omega}^2 + \|\chi B - H_h\|_{\mu, \Omega}^2)^{1/2} \|v\|_{\mu, \Omega},
\end{aligned}
$$

since $v$ was arbitrary in $L^2(\Omega)$.

□

**4.2. Magnetic field reconstruction.** We now present a practical construction of a discrete field $B_h \in \Lambda_{\Gamma_N}(\Omega)$ that is obtained as local post-processing of $H_h$. Our construction is based on divergence-constrained patch-wise minimization problems, and in fact closely follows the equilibrated-flux reconstruction of $D_h$ in electrostatic problems [11, 19, 25].

For each vertex $a \in \mathcal{V}_h$, the local contribution to our induction field reconstruction is given by

$$ (4.2) \qquad B_h^a := \arg \min_{\substack{v_h \in \mathcal{R}T_{p+2}(\mathcal{T}_h^a) \cap \mathcal{H}_0(\text{div}, \omega^a) \\ \nabla \cdot v_h = \nabla \psi^a \cdot (\mu H_h)}} \|v_h - \psi^a(\mu H_h)\|_{\chi, \omega^a} $$

This definition indeed makes sense: the compatibility condition is satisfied due to the first equation in (3.4), since $\nabla \psi^a \in \mathcal{N}_0(\mathcal{T}_h) \cap \mathcal{H}_0(\text{curl}, \Omega)$. Introducing

$$ (4.3) \qquad B_h := \sum_{a \in \mathcal{V}_h} B_h^a, $$

it is then clear that $B_h \in \mathcal{H}_{\Gamma_N}(\text{div}, \Omega)$ with $\nabla \cdot B_h = 0$.

**Remark 4.2 (Polynomial degree of the reconstruction).** The polynomial degree of the reconstruction is increased by two as compared to the finite element solution. This may seem surprising, as in the case of electrostatic problems, it is sufficient to increase the polynomial degree by one [11, 19, 25]. The reason for increasing to polynomial degree by one in electrostatic problems is because of the multiplication by the hat function $\psi^a$. Here, we have the additional issue that we are trying to represent $a$ (broken) Nédélec polynomial with Raviart-Thomas polynomials, which comes at the expense of another polynomial degree increase. It is possible to only increase the polynomial degree by one if Brezzi-Douglas-Marini are used instead of Raviart-Thomas elements in the reconstruction, but we do not pursue this approach here for the sake of simplicity.

For configurations where cohomology spaces are trivial, the above considerations suffice to infer that $B_h \in \Lambda_h$. In the general case however, we must additionally check that the periods of $B_h$ vanish according to (2.1). Surprisingly, this is automatically satisfied, as we next establish.

**Lemma 4.3 (Periods of the reconstruction).** Let $\Sigma \subset \Omega$ be a relatively open Lipschitz surface consisting of a collection of mesh faces such that $\partial\Sigma \subset \Gamma_N$. If $\Sigma$ if oriented with a unit normal vector $n_\Sigma$, then, we have

$$ \int_{\Sigma} B_h \cdot n_{\Sigma} = 0. $$

*Proof.* Let us denote by $\mathcal{V}_h^\Sigma$ the set of vertices $a \in \mathcal{V}_h$ belonging to $\bar{\Sigma}$. For each $a \in \mathcal{V}_h^\Sigma$, we let $\iota_a \in \{-, +\}$ be such that $\omega_{\Sigma, \iota_a}^a$ is non-empty and, if $a \in \bar{\Gamma}_N$, such that $\Gamma_{\Sigma, \iota_a}^a \subset \bar{\Gamma}_N$.
It is then clear that for each $a \in \mathcal{V}_h^\Sigma$

$$ \int_{\Sigma} \mathbf{B}_h^a \cdot \mathbf{n}_{\Sigma} = \iota_a \int_{\partial\omega_{\Sigma,\iota a}^a} \mathbf{B}_h^a \cdot \mathbf{n}_{\omega_{\Sigma,\iota a}^a} = \iota_a \int_{\omega_{\Sigma,\iota a}^a} \nabla \cdot \mathbf{B}_h^a \\ = \int_{\omega_{\Sigma,\iota a}^a} (\iota_a \nabla \psi^a) \cdot (\mu \mathbf{H}_h) = (\mu \mathbf{H}_h, \iota_a \chi_a \nabla \psi^a)_\Omega $$

where we let $\chi_a \in \mathcal{P}_0(\mathcal{T}_h)$ denote the set function of $\omega_{\Sigma,\iota a}^a$. As a result, introducing the test function

$$ (4.4) \qquad \mathbf{v}_h := \sum_{a \in \mathcal{V}_h^\Sigma} \iota_a \chi_a \nabla \psi^a \in \mathbf{L}^2(\Omega), $$

we have

$$ \int_{\Sigma} \mathbf{B}_h^a \cdot \mathbf{n}_{\Sigma} = (\mu \mathbf{H}_h, \mathbf{v}_h)_{\Omega}, $$

and using the first equation of (3.4), the conclusion follows if we can show that $\mathbf{v}_h \in \mathbf{V}_h$ with $\nabla \times \mathbf{v}_h = \mathbf{0}$.

We first observe that in each cell $K \in \mathcal{T}_h$, $(\iota_a \chi_a)|_K$ is constant, so that

$$ \mathbf{v}_h|_K = \sum_{a \in \mathcal{V}_h^\Sigma} (\iota_a \chi_a)|_K \nabla(\psi^a)|_K = \nabla \tilde{q}_K $$

for some $\tilde{q}_K \in \mathcal{P}_1(K)$. As a result $\mathbf{v}_h \in \mathcal{N}_0(\mathcal{T}_h)$ with $\nabla \times \mathbf{v}_h = \mathbf{0}$ element-wise, and therefore, the conclusion follows if we can show that $[\mathbf{v}_h]_F \times \mathbf{n}_F = \mathbf{0}$ for all $F \in \mathcal{F}_h^i$ and $\mathbf{v}_h \times \mathbf{n} = \mathbf{0}$ on $\Gamma_T$.

Let us first consider an interior face $F \in \mathcal{F}_h^i$. Notice that then $[\nabla \psi^a]_F \times \mathbf{n}_F = \mathbf{0}$ for all $a \in \mathcal{V}_h$, and we have

$$ [\mathbf{v}_h]_F \times \mathbf{n}_F = \sum_{a \in \mathcal{V}_h^\Sigma} \iota_a [\chi_a]_F \nabla \psi^a \times \mathbf{n}_F = \sum_{a \in \mathcal{V}(F)} \iota_a [\chi_a]_F \nabla \psi^a \times \mathbf{n}_F. $$

We then distinguish two cases. (i) If $F \not\subset \bar{\Sigma}$, for each $a \in \mathcal{V}(F)$, then either (ia) $F \subset \partial\omega^a$, in which case, $\nabla \psi^a \times \mathbf{n}_F = \mathbf{0}$ and $[\mathbf{v}_h]_F \times \mathbf{n}_F = \mathbf{0}$ or (ib) $F \not\subset \partial\omega^a$, in which case $[\chi_a] = 0$, again leading to $[\mathbf{v}_h]_F \times \mathbf{n}_F = \mathbf{0}$. On the other hand, (ii) if $F \subset \Sigma$, we have

$$ [\mathbf{v}_h]_F \times \mathbf{n}_F = \sum_{a \in \mathcal{V}(F)} \nabla \psi^a \times \mathbf{n}_F $$

since $[\chi_a]_F = \iota_a$ for all the relevant vertices in the sum. But then, we easily conclude since

$$ [\mathbf{v}_h]_F \times \mathbf{n}_F = \nabla \left( \sum_{a \in \mathcal{V}(F)} \psi^a \right) \times \mathbf{n}_F = \nabla (1) \times \mathbf{n}_F = \mathbf{0}. $$

We conclude the proof by considering a face $F \subset \bar{\Gamma}_T$, whereby we need to show that $\mathbf{v}_h \times \mathbf{n}_F = \mathbf{0}$ on $F$. We first observe that due to the locality of the hat functions $\psi^a$, the value of $\mathbf{v}_h \times \mathbf{n}_F$ on $F$ can only be influenced by boundary vertices $a \subset \partial\Omega$ in (4.4). Since on the other hand $\partial\Sigma \subset \bar{\Gamma}_N$, such vertices must lie in $\bar{\Gamma}_N$. However, for such vertices, we have always chosen the half-patch $\omega_{\Sigma,\iota a}^a$ such that $\Gamma_{\Sigma,\iota a}^a \subset \bar{\Gamma}_N$, meaning that corresponding contribution vanishes on $\Gamma_T$, as $\chi_a = 0$ on $\Gamma_T$. $\square$
A direct consequence of the characterization (2.1) of the image of the curl, the assumption (3.1) that the surfaces $\{\Sigma_j\}_{j=1}^N$ are aligned with the mesh $T_h$ and Lemma 4.3 above is that $B_h$ is indeed the curl of some potential.

**Theorem 4.4** (Localized reconstruction). For the induction field reconstruction provided by (4.2) and (4.3), we have $B_h \in \Lambda_{\Gamma_N}(\Omega)$.

## 5. EFFICIENCY

In this section, we present efficiency results for the estimator constructed with the procedure of Section 4.2. We start with a patch-wise efficiency result that applies to each local contribution $B_h^a$. The proof is similar to seminal results for the electrostatic problem, see [12, 25].

**Theorem 5.1** (Lower-bound). The estimate

$$ (5.1) \qquad \|B_h^a - \psi^a(\mu H_h)\|_{\chi, \omega^a} \le C_{\text{st},a} C_{\text{cont},a} \mathcal{C}_{\mu, \omega^a} \|H - H_h\|_{\mu, \omega^a} $$

holds true for all vertices $a \in \mathcal{V}_h$.

*Proof.* Let $B_h^a$ solve the continuous version of (4.2), i.e.

$$ (5.2) \qquad B_h^a := \arg \min_{\substack{v \in H_0(\text{div}, \omega^a) \\ \nabla \cdot v = \nabla \psi^a(\mu H_h)}} \|v - \psi^a(\mu H_h)\|_{\chi, \omega^a}. $$

Then, (3.3) ensures that

$$ (5.3) \qquad \|B_h^a - \psi^a(\mu H_h)\|_{\chi, \omega^a} \le C_{\text{st},a} \mathcal{C}_{\mu, \omega^a}^{1/2} \|B_h^a - \psi^a(\mu H_h)\|_{\chi, \omega^a}. $$

Indeed, $H_h \in \mathcal{N}_p(T_h) \subset \mathcal{P}_{p+1}(T_h)$, so that $\psi^a(\mu H_h) \in \mathcal{P}_{p+2}(T_h) \subset \mathcal{R}T_{p+2}(T_h)$. As a result, we only need to focus on the continuous minimization problem (5.2). Notice that the additional factor $\mathcal{C}_{\mu, \omega^a}^{1/2}$ as compared to (3.3) is due to the weighted norms.

The Euler-Lagrange equations for (5.2) show that there exists a unique $r^a \in L^2_*(\omega^a)$ such that

$$ \left\{ \begin{array}{l} (\chi B_h^a, v)_{\omega^a} + (r^a, \nabla \cdot v)_{\omega^a} = (\psi^a H_h, v)_{\omega^a} \\ (\nabla \cdot B_h^a, q)_{\omega^a} = (\nabla \psi^a \cdot (\mu H_h), q)_{\omega^a} \end{array} \right. $$

for all $v \in H_0(\text{div}, \omega^a)$ and $q \in L^2_*(\omega^a)$. Integrating by parts in the first equation shows that we actually have $r^a \in H^1_*(\omega^a)$ with

$$ \mu \nabla r^a = B_h^a - \psi^a(\mu H_h), $$

and in particular

$$ \|\nabla r^a\|_{\mu, \omega^a} = \|\chi(B_h^a - \psi^a(\mu H_h))\|_{\mu, \omega^a} = \|B_h^a - \psi^a(\mu H_h)\|_{\chi, \omega^a}. $$

To complete the proof, we estimate the norm of $\nabla r^a$. We start with

$$ \begin{aligned} (\mu \nabla r^a, \nabla v)_{\omega^a} &= (B_h^a, \nabla v)_{\omega^a} - (\psi^a(\mu H_h), \nabla v)_{\omega^a} \\ &= -(\nabla \cdot B_h^a, v)_{\omega^a} - (\psi^a(\mu H_h), \nabla v)_{\omega^a} \\ &= -(\nabla \psi^a \cdot (\mu H_h), v)_{\omega^a} - (\psi^a(\mu H_h), \nabla v)_{\omega^a} \\ &= -(\mu H_h, \nabla (\psi^a v))_{\omega^a}, \end{aligned} $$

and upon recalling that $\nabla \cdot (\mu H) = 0$, we arrive at

$$ (\mu \nabla r^a, \nabla v)_{\omega^a} = (\mu(H - H_h), \nabla(\psi^a v))_{\omega^a}. $$
At that point, we easily finish the proof with (3.2) since

$$ \|\nabla r^a\|_{\mu, \omega^a}^2 \le \|H - H_h\|_{\mu, \omega^a} \|\nabla(\psi^a r^a)\|_{\mu, \omega^a} \le C_{\text{cont},a} \mathscr{C}_{\mu, \omega^a}^{1/2} \|H - H_h\|_{\mu, \omega^a} \|\nabla r^a\|_{\mu, \omega^a}, $$

so that

$$ \|\nabla r^a\|_{\mu, \omega^a} \le C_{\text{cont},a} \mathscr{C}_{\mu, \omega^a}^{1/2} \|H - H_h\|_{\mu, \omega^a}. $$

□

From the efficiency of the local contribution $B_h^a$, we finally establish an efficiency result for the induction field reconstruction $B_h$.

**Corollary 5.2** (Lower-bounds for the element-wise estimator). The efficiency estimate

$$ (5.4) \qquad \|B_h - \mu H_h\|_{\chi,K} \le 2C_{\text{st},K} C_{\text{cont},K} \mathscr{C}_{\mu, \tilde{K}} \|H - H_h\|_{\mu, \tilde{K}} $$

holds true for all $K \in \mathcal{T}_h$ with

$$ \tilde{K} := \bigcup_{a \in \mathcal{V}(K)} \omega^a, \quad C_{\text{st},K} := \max_{a \in \mathcal{V}(K)} C_{\text{st},a}, \quad C_{\text{cont},K} := \max_{a \in \mathcal{V}(K)} C_{\text{cont},K}. $$

*Proof.* Let $K \in \mathcal{T}_h$. The four hat functions associated with the vertices of K form a partition of unity over K, and $B_h^a = o$ on K unless $a \in \mathcal{V}(K)$. As a result, we have

$$ B_h - \mu H_h = \sum_{a \in \mathcal{V}(K)} B_h^a - \psi^a(\mu H_h) $$

on K, and the triangle inequality immediately gives that

$$ \|B_h - \mu H_h\|_{\chi,K} = \left\| \sum_{a \in \mathcal{V}(K)} B_h^a - \psi^a(\mu H_h) \right\|_{\chi,K} \le \sum_{a \in \mathcal{V}(K)} \|B_h^a - \psi^a(\mu H_h)\|_{\chi,K} \\ \le \sum_{a \in \mathcal{V}(K)} \|B_h^a - \psi^a(\mu H_h)\|_{\chi, \omega^a}. $$

The estimate in (5.4) then follows from (5.1) since

$$ \begin{aligned} \|B_h - \mu H_h\|_{\chi,K} &\le \sum_{a \in \mathcal{V}(K)} C_{\text{st},a} C_{\text{cont},a} \mathscr{C}_{\mu, \omega^a} \|H_h - H_h\|_{\mu, \omega^a} \\ &\le C_{\text{st},K} C_{\text{cont},K} \mathscr{C}_{\mu, \tilde{K}} \sum_{a \in \mathcal{V}(K)} \|H_h - H_h\|_{\mu, \omega^a} \\ &\le 2C_{\text{st},K} C_{\text{cont},K} \mathscr{C}_{\mu, \tilde{K}} \|H_h - H_h\|_{\mu, \tilde{K}}, \end{aligned} $$

as K has 4 vertices.

□

# 6. NUMERICAL EXAMPLES

## 6.1. Settings.
In all numerical examples, we consider the case where $\Gamma_N = \emptyset$, and thus, $\Gamma_T = \partial\Omega$. For simplicity, we also set $\mu = I$. Notice that then $H = B$.

In the graphs below, we employ the notations

$$ \text{err}_H := \|H - H_h\|_{\Omega}, \quad \text{err}_B := \|H - B_h\|_{\Omega}, \quad \eta := \|H_h - B_h\|_{\Omega}, $$

where $H_h$ is the finite element solution, and $B_h$ is the reconstruction introduced in Section 4.2. Observe that for the sake of simplicity, we have not introduced the oscillation term in $\eta$.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Trivial cohomology, $p = 0$</figcaption></figure>

Element-wise versions of err<sub>H</sub> and η are also plotted in the adaptivity example. The notation $N_{\text{dofs}}$ stands for the dimension of the Nédélec finite element space $V_h$.

We employ gmsh to generate meshes [30], and mmg3d for adaptive mesh refinements [21]. The linear systems are solved using the mumps software package [6].

6.2. **Trivial cohomology.** Here, $\Omega = (0,1)^3$ is the unit cube. The right-hand side and solution respectively read

$$ (6.1) \qquad \mathbf{J} := \begin{pmatrix} \sin(\pi \mathbf{x}_1) \cos(\pi \mathbf{x}_2) \\ -\cos(\pi \mathbf{x}_1) \sin(\pi \mathbf{x}_2) \\ 0 \end{pmatrix}, \qquad \mathbf{H} := \begin{pmatrix} 0 \\ 0 \\ \sin(\pi \mathbf{x}_1) \sin(\pi \mathbf{x}_2) \end{pmatrix}. $$

In Figures 1 and 2, we respectively fix the polynomial degree to $p = 0$ and $1$, and let the mesh size vary. The expected convergence rates (namely linear and quadratic) are observed for the discrete solution $H_h$, and the estimator closely follows the error. Indeed, as can be seen on the right panels of Figures 1 and 2, the effectivity indices stay slightly above one. We also see that the reconstruction $B_h$ converges toward $B$ with the same rates as $H_h$ converges to $H$. The post-processing is about 3 to 5 times more accurate than the original discrete solution in this case.

Figure 3 displays a $p$-convergence result whereby a mesh of 176 tetrahedra is fixed ($h_{\max} = 0.76$), and $p$ varies from 0 to 5. On the left panel, we note that the convergence is exponential in the number of dofs, as is to be expected. The post-processed solution $B_h$ is 2 to 4 times more accurate than $H_h$. The behaviour of the effectivity index is reported on the right-panel. It slightly drops below one because the mesh is fairly coarse and the data oscillation term have not been computed. Yet, the agreement between the error and the estimator is excellent, and the effectivity index remains essentially constant for all polynomial degrees, highlighting the $p$-robustness of the estimator.

6.3. **Non-trivial cohomology.** We then consider an example where the cohomology space is non-trivial by setting $\Omega := R \setminus C$ with $R := (-2, 2) \times (-1, 1) \times (-2, 2)$ and $C := (-1, 1)^3$. $\Omega$ thus has the topology of a torus. To easily construct $\Lambda_h$, we require that the surface $\Sigma := (1, 2) \times (-1, 1) \times \{0\}$ is exactly meshed, and impose that $(w_h \cdot n, 1)_\Sigma = 0$ for all
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Trivial cohomology example, $p = 1$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3. Trivial cohomology example, $p$-convergence</figcaption></figure>

$w_h \in \Lambda_h$. Figures 4 sketches the geometry of the domain. The right-hand side and solution are still defined with (6.1).

Figures 5 and 6 show the behaviour of the post-processed induction field and the estimator for $p = 0$ and $1$ respectively. The behaviour is similar than in the trivial cohomology example, with the difference that the error is slightly underestimated on the first mesh. This is to be expected since the data oscillation term is not included and starting mesh is coarser in this example than in the trivial topology case.

Figure 7 displays a $p$-convergence result on a fixed mesh with 1043 tetrahedra ($h_{max} = 1.19$). As before, we observe that the estimator is $p$-robust. The error is slightly underestimated due to the coarseness of the mesh.

6.4. **Adaptivity.** We close this section with an example where the estimator is employed to drive an adaptive mesh refinement process with Dörfler marking [22]. After the solution and the estimator have been computed, the elements are sorted according to the estimator values,
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Geometry of the non-trivial topology example: cut at $x_2 = 0$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Non-trivial cohomology example, $p = 0$</figcaption></figure>

and we select a minimal set $\mathcal{M}_h$ of elements such that

$$ \sum_{K \in \mathcal{M}_h} \eta_K^2 \le \theta^2 \sum_{K \in \mathcal{T}_h} \eta_K^2, \quad \theta := 0.05. $$

The software `mmg3d` [21] is used to iteratively refine the mesh locally around the marked elements $K \in \mathcal{M}_h$, with the procedure detailed in [15, Section 6.1.4].
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Non-trivial cohomology example, $p = 1$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7. Non-trivial cohomology example, $p$-convergence</figcaption></figure>

The computational domain is the “L-brick” $\Omega := L \times (0, 2)$ where $L = (-2, 2)^2 \setminus (0, 2) \times (-2, 0)$. The right-hand side is given by

$$J := \begin{pmatrix} (x_2 - y_2) \\ -(x_1 - y_1) \\ 0 \end{pmatrix} e^{-|x-y|^2/\sigma^2},$$

with $y = (-1, 1, 1)$ and $\sigma = 0.4$. Notice that this right-hand side is divergence-free but does exactly satisfy the boundary condition $J \cdot n = 0$ on $\partial\Omega$. However, the normal trace is so small that the discrepancy does not affect the numerical results for the mesh sizes considered. The analytical solution $H$ is not available here, so that a reference solution $\tilde{H}$ is computed for each mesh by increasing the polynomial degree by 1.

The adaptive refinement process is started with an initial mesh of 56 tetrahedra ($h_{max} = 2.83$). We perform two simulations with $p = 0$ and $p = 1$. In Figure 8, we compare the convergence rates obtained using the estimator $\eta_K$ and the “true” error $\|\tilde{H} - H_h\|_K$ to mark the elements in the adaptivity refinements. The results are extremely similar, which shows
<figure><img src="image_9.png" /><figcaption>FIGURE 8. err_H in error-driven and estimator-driven adaptivity</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9. Adaptivity example, $p = 0$</figcaption></figure>

that the estimator is perfectly suited to drive the refinements. The optimal convergence rate of $N_{\text{dofs}}^{p/3}$ are also achieved.

Figures 9 and 10 display the behaviour of the error, the estimator and the post-processed numerical solution throughout the adaptive process. As before, we observe that the post-processed solution is more accurate. The effectivity index is very close to one, except on the coarsest mesh where the error is underestimated due to the absence of the data oscillation term.

Figures 11, 12, 13 and 14 represent the meshes obtained at different stages of the iterative refinement process. Specifically the top size of the L-brick is represented on the top of the figures, whereas the two faces sharing the reentrant edges are represented at the bottom. The agreement between the estimator $\eta_K$ and the local error $||\tilde{H} - H_h||_K$ is excellent. We also observe that the mesh is locally refined as expected: close to the source center $\mathbf{y}$ and along the reentrant edge.
<figure><img src="image_11.png" /><figcaption>FIGURE 10. Adaptivity example, $p = 1$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11. Adaptivity example with $p = 0$, iteration #5</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIGURE 12. Adaptivity example with $p = 0$, iteration #15</figcaption></figure>

## 7. CONCLUSION

We propose an a posteriori error estimator for mixed finite element discretizations of the curl-curl problem. Following the framework of Prager-Synge estimates, we engineer a potential reconstruction leading to guaranteed error bounds. A key novelty of our approach is that instead of explicitly constructing the potential, we actually build its curl through patchwise divergence-constrained minimization problems. This construction leads to an estimator that is locally efficient and polynomial-degree-robust. It is noteworthy that the reconstruction technique proposed here for the curl-curl problem is similar to the construction of equilibrated fluxes for the Poisson problem.

We highlight the properties of the reconstructed field and associated estimator with selected numerical examples. These examples show that the reconstructed induction field is typically more accurate than the discrete magnetic field output by the finite element scheme, although not being super-convergent. The key theoretical properties are also observed numerically: the proposed estimator is efficient, polynomial-degree-robust, and fully reliable up to data oscillations. We also employ the estimator to drive adaptive mesh refinements in a domain with a reentrant edge. This example indicates that the estimator is suited for adaptivity purposes: the meshes produced are correctly refined, leading to optimal convergence rates.
<figure><img src="image_14.png" /><figcaption>FIGURE 13. Adaptivity example with $p = 1$, iteration #10</figcaption></figure>

REFERENCES

1. R. Adams and J. Fournier, Sobolev spaces, Academic Press, 2003.

2. M. Ainsworth, A posteriori error estimation for lowest order Raviart-Thomas mixed finite elements, SIAM J. Sci. Comput. **30** (2008), no. 1, 189-204.

3. M. Ainsworth and J.T. Oden, A posteriori error estimation in finite element analysis, Wiley, 2000.

4. A. Alonso Rodríguez, J. Camaño, E. De los Santos, and F. Rapetti, A graph approach for the construciton of high order divergre-free raviart-thomas finite elements, Calcolo **55** (2018).

5. P. Alotto and I. Perugia, Mixed finite element methods and tree-cotree implicit condensation, Calcolo **26** (1999), 233-248.

6. P.R. Amestoy, I.S. Duff, and J.Y. L'Excellent, Multifrontal parallel distributed symmetric and unsymmetric solvers, Comput. Methods Appl. Mech. Engrg. **184** (2000), 501-520.

7. C. Amrouche, C. Bernardi, M. Dauge, and V. Girault, Vector potentials in three-dimensional non-smooth domains, Math. Meth. Appl. Sci. **21** (1998), 823-864.

8. T. Apel, Anisotropic finite elements: local estimates and applications, 1999.

9. F. Assous, P. Ciarlet, and S. Labrunie, Mathematical foundations of computational electromagnetism, Springer, 2018.

10. D. Boffi, F. Brezzi, and M. Fortin, Mixed and hybrid finite element methods and applications, Springer, 1991.

11. D. Braess, V. Pillwein, and J. Schöberl, Equilibrated residual error estimates are p-robust, Comput. Meth. Appl. Mech. Engrg. **198** (2009), 1189-1197.
<figure><img src="image_15.png" /><figcaption>FIGURE 14. Adaptivity example with $p = 1$, iteration #20</figcaption></figure>

12. D. Braess and J. Schöberl, *Equilibrated residual error estimators for edge elements*, Math. Comp. **77** (2008), no. 262, 651–672.

13. S.C. Brenner, *Poincaré-Friedrichs inequalities for piecewise H¹ functions*, SIAM J. Numer. Anal. **41** (2003), no. 1, 306–324.

14. F. Brezzi, *On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers*, ESAIM Math. Model. Numer. Anal. **8** (1974), no. 2, 129–151.

15. T. Chaumont-Frelet, *A simple equilibration procedure leading to polynomial-degree-robust a posteriori error estimators for the curl-curl problem*, Math. Comp. (2023).

16. T. Chaumont-Frelet and M. Vohralík, *Constrained and unconstrained stable discrete minimizations for p-robust local reconstructions in vertex patches in the De Rham complex*, HAL preprint hal-03749682, 2022.

17. —, *p-robust equilibrated flux reconstruction in h(curl) based on local minimizations. application to a posteriori analysis of the curl-curl problem*, HAL preprint hal-03227570, 2022.

18. M. Costabel and M. Dauge, *Singularities of electromagnetic fields in polyhedral domains*, Arch. Ration. Mech. **151** (2000), 221–276.

19. P. Destuynder and B. Métivet, *Explicit error bounds in a conforming finite element method*, Math. Comp. **68** (1999), no. 228, 1379–1396.

20. P. Dlotko, B. Kapidani, and R. Specogna, *TOPOPROCESSOR: An efficient computational topology toolbox for h-oriented eddy current formulations*, IEE Trans. Magn. **53** (2017), no. 6, 7204404.

21. C. Dobrzynski, *MMG3D: User guide*, Tech. Report 422, Inria, 2012.
22. W. Dörfler, *A convergent adaptive algorithm for Poisson's equation*, SIAM J. Numer. Anal. **33** (1996), 1106–1124.

23. G.F.D. Duff, *Differnetial forms in manifolds with boundary*, Ann. Math. **56** (1952), no. 1, 115–127.

24. A. Ern and J.-L. Guermond, *Finite elements I: Approximation and interpolation*, Texts in Applied Mathematics, vol. 72, Springer Nature, Cham, Switzerland, 2021.

25. A. Ern and M. Vohralík, *Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations*, SIAM J. Numer. Anal. **53** (2015), no. 2, 1058–1081.

26. —, *Stable broken H¹ and H(div) polynomial extensions for polynomial-degree-robust potential and flux reconstruction in three space dimensions*, Math. Comp. **89** (2021), 551–594.

27. P. Fernandes and G. Gilardi, *Magnetostatic and electrostatic problems in inhomogeneous anisotropic media with irregular boundary and mixed boundary conditions*, Math. Meth. Appl. Sci. **47** (1997), no. 4, 2872–2896.

28. J. Gedicke, S. Geevers, and I. Perugia, *An equilibrated a posteriori error estimator for arbitrary-order Nédélec elements for magnetostatic problems*, J. Sci. Comput. **83** (2020), no. 3, Paper No. 58, 23.

29. J. Gedicke, S. Geevers, I. Perugia, and J. Schöberl, *A polynomial-degree-robust a posteriori error estimator for Nédélec discretizations of magnetostatic problems*, SIAM J. Numer. Anal. **59** (2021), no. 4, 2237–2253.

30. C. Geuzaine and J.F. Remacle, *Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities*, Int. J. Numer. Meth. Engrg. **79** (2009), 1309–1331.

31. V. Girault and P.A. Raviart, *Finite element methods for Navier-Stokes equations: theory and algorithms*, Springer-Verlag, 1986.

32. D.J. Griffiths, *Introduction to Electrodynamics*, Prentice Hall, 1999.

33. P.W. Gross and P.R. Kotiuga, *Electromagnetic theory and computation: a topological approach*, Cambridge University Press, 2004.

34. P. Knobloch, *Uniform validity of discrete Friedrichs' inequality for general nonconforming finite element spaces*, Numer. Funct. Anal. Optim. **22** (2001), no. 1-2, 107–126.

35. P. Ladevève and D. Leguillon, *Error estimate procedure in the finite element method and applications*, SIAM J. Numer. Anal. **20** (1983), 485–509.

36. R. Luce and B.I. Wohlmuth, *A local a posteriori error estimator based on equilibrated fluxes*, SIAM J. Numer. Anal. **42** (2004), 1394–1414.

37. J.C. Nédélec, *Mixed finite elements in $\mathbb{R}^3$*, Numer. Math. **35** (1980), 315–341.

38. S. Nicaise, *Edge elements on anisotropic meshes and approximation of the Maxwell equations*, SIAM J. Numer. Anal. **39** (2001), no. 3, 784–816.

39. S. Nicaise, K. Witowski, and B.I. Wohlmuth, *An a posteriori error estimator for the Lamé equation based on equilibrated fluxes*, IMA J. Numer. Anal. **28** (2008), 331–353.

40. M. Pellikka, S. Suuriniemi, L. Kettunen, and C. Geuzaine, *Homology and cohomology computation in finite element modeling*, SIAM J. Sci. Comput. **35** (2013), no. 5, B1195–B1214.

41. W. Prager and J.L. Synge, *Approximations in elasticity based on the concept of function space*, Quart. Appl. Math. **5** (1947), no. 3, 241–269.

42. P.A. Raviart and J.M. Thomas, *A mixed finite element method for 2nd order elliptic problems*, Mathematical Aspect of Finite Element Methods, Springer-Verlag, 1977.

43. R. Scheichl, *Decoupling three-dimensional mixed problems using divergence-free finite elements*, SIAM J. Sci. Comput. **23** (2002), no. 5, 1752–1776.

44. M. Vohralík, *On discrete Poincaré-Friedrichs inequalities for nonconforming approximation of the Sobolev space H¹*, Numer. Funct. Anal. Optim. **26** (2005), 925–952.

45. —, *A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations*, SIAM J. Numer. Anal. **45** (2007), no. 4, 1570–1599.