1 Introduction

We consider an a posteriori error estimator for finite element methods for solving equations of the form \(\nabla \times \mu ^{-1}(\nabla \times {\mathbf{u}})={\mathbf{j}}\). These equations are related to magnetostatics but also appear in eddy current models for non-conductive media.

The first a posteriori error estimator in this context was introduced and analysed in [4]. It is a residual-type estimator and provides bounds of the form

$$\begin{aligned} c_0\, \text {estimator} \le \text {error}&\le c_1\,\text {estimator}, \end{aligned}$$

up to some higher-order data oscillation terms, where \(c_0,c_1\) are positive constants that do not depend on the mesh resolution. Similar bounds can be obtained by hierarchical error estimators; see, e.g., [5], under the assumption of a saturation condition, and by Zienkiewicz–Zhu-type error estimators; see, e.g., [19]. A drawback of these estimators is that the constants \(c_0\) and \(c_1\) are usually unknown, resulting in significant overestimation or underestimation of the real error.

Equilibration-based error estimators can circumvent this problem. Often attributed to Prager and Synge [20], these estimators have become a major research topic; for a recent overview, see, for example, [13] and the references therein. An equilibration-based error estimator was introduced for magnetostatics in [6] and provides bounds of the form

$$\begin{aligned} c_0\, \text {estimator} \le \text {error}&\le \text {estimator}, \end{aligned}$$

up to some higher-order data oscillation terms. In other words, it provides a constant-free upper bound on the error. A different equilibration-based error estimator for magnetostatics was introduced in [21] and, for an eddy current problem, in [10, 11]. Constant-free upper bounds are also obtained by the functional estimate in [18], when selecting a proper function y in their estimator, and by the recovery-type error estimator in [7], in case the equations contain an additional term \(\beta {\mathbf{u}}\), with \(\beta >0\).

A drawback of the estimators in [10, 11, 18, 21] is that they require solving a global problem. The estimator in [6], on the other hand, only involves solving local problems related to vertex patches. However, the latter estimator is defined for Nédélec elements of the lowest degree only. In this paper, we present a new equilibration-based constant-free error estimator that can be applied to Nédélec elements of arbitrary degree. Furthermore, our estimator involves solving problems on only single elements, single faces, and very small sets of nodes.

The paper is constructed as follows: We firstly introduce a finite element method for solving magnetostatic problems in Sect. 2. We then derive our error estimator step by step in Sect. 3, with a summary given in Sect. 3.3, and prove its reliability and efficiency in Sect. 3.4. Numerical examples confirming the reliability and efficiency of our estimator are presented in Sect. 4, and an overall summary is given in Sect. 5.

2 A Finite Element Method for Magnetostatic Problems

Let \(\varOmega \subset {\mathbb {R}}^3\) be an open, bounded, simply connected, polyhedral domain with a connected Lipschitz boundary \(\partial \varOmega \). In case of a linear, isotropic medium and a perfectly conducting boundary, the static magnetic field \({\mathbf{H}}:\varOmega \rightarrow {\mathbb {R}}^3\) satisfies the equations

$$\begin{aligned} \nabla \times {\mathbf{H}}&= {\mathbf{j}}&\text {in }\varOmega , \\ \nabla \cdot \mu {\mathbf{H}}&= 0&\text {in }\varOmega , \\ \hat{{\mathbf{n}}}\cdot \mu {\mathbf{H}}&= 0&\text {on }\partial \varOmega , \end{aligned}$$

where \(\nabla \) is the vector of differential operators \((\partial _1,\partial _2,\partial _3)\), \(\times \) and \(\cdot \) denote the outer- and inner product, respectively, (therefore, \(\nabla \times \) and \(\nabla \cdot \) are the curl- and divergence operator, respectively), \(\hat{{\mathbf{n}}}\) denotes the outward pointing unit normal vector, \(\mu :\varOmega \rightarrow {\mathbb {R}}^+\), with \(\mu _0\le \mu \le \mu _1\) for some positive constants \(\mu _0\) and \(\mu _1\), is a scalar magnetic permeability, and \({\mathbf{j}}:\varOmega \rightarrow {\mathbb {R}}^3\) is a given divergence-free current density. The first equality is known as Ampère’s law and the second as Gauss’s law for magnetism.

These equations can be solved by writing \({\mathbf{H}}=\mu ^{-1}\nabla \times {\mathbf{u}}\), where \({\mathbf{u}}:\varOmega \rightarrow {\mathbb {R}}^3\) is a vector potential, and by solving the following problem for \({\mathbf{u}}\):

$$\begin{aligned} \nabla \times (\mu ^{-1}\nabla \times {\mathbf{u}})&= {\mathbf{j}}&\text {in }\varOmega , \end{aligned}$$
(1a)
$$\begin{aligned} \nabla \cdot {\mathbf{u}}&= 0&\text {in }\varOmega , \end{aligned}$$
(1b)
$$\begin{aligned} \hat{{\mathbf{n}}}\times {\mathbf{u}}&= {\mathbf{0}}&\text {on }\partial \varOmega . \end{aligned}$$
(1c)

The second condition is only added to ensure uniqueness of \({\mathbf{u}}\) and is known as Coulomb’s gauge.

Now, for any domain \(D\in {\mathbb {R}}^n\), let \(L^2(D)^m\) denote the standard Lebesque space of square-integrable vector-valued functions \({\mathbf{u}}:D\rightarrow {\mathbb {R}}^m\) equipped with norm \(\Vert {\mathbf{u}}\Vert _D^2:=\int _D {\mathbf{u}}\cdot {\mathbf{u}}\;\mathrm {d}{{\mathbf{x}}}\) and inner product \(({\mathbf{u}},{\mathbf{w}})_D=\int _D {\mathbf{u}}\cdot {\mathbf{w}}\;\mathrm {d}{{\mathbf{x}}}\), and define the following Sobolev spaces:

$$\begin{aligned} H^1(\varOmega )&:= \{\phi \in L^2(\varOmega ) \;|\; \nabla \phi \in L^2(\varOmega )^3 \}, \\ H_0^1(\varOmega )&:= \{\phi \in H^1(\varOmega ) \;|\; \phi =0 \text { on }\partial \varOmega \}, \\ H(\mathrm {curl}; \varOmega )&:= \{{\mathbf{u}}\in L^2(\varOmega )^3 \;|\; \nabla \times {\mathbf{u}}\in L^2(\varOmega )^3 \}, \\ H_0(\mathrm {curl}; \varOmega )&:= \{{\mathbf{u}}\in H(\mathrm {curl};\varOmega ) \;|\; \hat{{\mathbf{n}}}\times {\mathbf{u}}={\mathbf{0}}\text { on }\partial \varOmega \}, \\ H(\mathrm {div}; \varOmega )&:= \{{\mathbf{u}}\in L^2(\varOmega )^3 \;|\; \nabla \cdot {\mathbf{u}}\in L^2(\varOmega ) \}, \\ H(\mathrm {div}^0; \varOmega )&:= \{{\mathbf{u}}\in H(\mathrm {div};\varOmega ) \;|\; \nabla \cdot {\mathbf{u}}=0 \}. \end{aligned}$$

The weak formulation of problem (1) is finding \({\mathbf{u}}\in H_0(\mathrm {curl};\varOmega )\cap H(\mathrm {div}^0;\varOmega )\) such that

$$\begin{aligned} (\mu ^{-1}\nabla \times {\mathbf{u}},\nabla \times {\mathbf{w}})_{\varOmega }&= ({\mathbf{j}},{\mathbf{w}})_{\varOmega }&\forall {\mathbf{w}}\in H_0(\mathrm {curl};\varOmega ), \end{aligned}$$
(2)

which is a well-posed problem [15, Theorem 5.9].

The solution of the weak formulation can be approximated using a finite element method. Let T be a tetrahedron and define \(P_k(T)\) to be the space of polynomials on T of degree k or less. Also, define the Nédélec space of the first kind \(R_k(T)\) and the Raviart-Thomas space \(D_k(T)\) by

$$\begin{aligned} R_k(T)&:= \{{\mathbf{u}}\in P_k(T)^3 \;|\; {\mathbf{u}}({\mathbf{x}}) = {\mathbf{v}}({\mathbf{x}}) + {\mathbf{x}}\times {\mathbf{w}}({\mathbf{x}}) \text { for some }{\mathbf{v}},{\mathbf{w}}\in P_{k-1}(T)^3\}, \\ D_k(T)&:= \{{\mathbf{u}}\in P_k(T)^3 \;|\; {\mathbf{u}}({\mathbf{x}}) = {\mathbf{v}}({\mathbf{x}}) + {\mathbf{x}}w({\mathbf{x}}) \text { for some }{\mathbf{v}}\in P_{k-1}(T)^3, w\in P_{k-1}(T)\}. \end{aligned}$$

Finally, let \(\mathcal {T}_h\) denote a tessellation of \(\varOmega \) into tetrahedra with a diameter smaller than or equal to h, let \(P_k^{-1}(\mathcal {T}_h)\), \(R_k^{-1}(\mathcal {T}_h)\), and \(D_k^{-1}(\mathcal {T}_h)\) denote the discontinuous spaces given by

$$\begin{aligned} P_k^{-1}(\mathcal {T}_h):= & {} \{\phi \in L^2(\varOmega ) \;|\; \phi |_T \in P_k(T) \text { for all }T\in \mathcal {T}_h\}, \\ R_k^{-1}(\mathcal {T}_h):= & {} \{{\mathbf{u}}\in L^2(\varOmega )^3 \;|\; {\mathbf{u}}|_T \in R_k(T) \text { for all }T\in \mathcal {T}_h\}, \\ D_k^{-1}(\mathcal {T}_h):= & {} \{{\mathbf{u}}\in L^2(\varOmega )^3 \;|\; {\mathbf{u}}|_T \in D_k(T) \text { for all }T\in \mathcal {T}_h\}, \end{aligned}$$

and define

$$\begin{aligned} P_k(\mathcal {T}_h)&:=P^{-1}_k(\mathcal {T}_h)\cap H^1(\varOmega ),&P_{k,0}(\mathcal {T}_h)&:=P^{-1}_k(\mathcal {T}_h)\cap H^1_0(\varOmega ),\\ R_k(\mathcal {T}_h)&:=R^{-1}_k(\mathcal {T}_h)\cap H(\mathrm {curl};\varOmega ),&R_{k,0}(\mathcal {T}_h)&:=R^{-1}_k(\mathcal {T}_h)\cap H_0(\mathrm {curl};\varOmega ),\\ D_k(\mathcal {T}_h)&:=D^{-1}_k(\mathcal {T}_h)\cap H(\mathrm {div};\varOmega ).&\end{aligned}$$

We define the finite element approximation for the magnetic vector potential as the vector field \({\mathbf{u}}_h\in R_{k,0}(\mathcal {T}_h)\) that solves

$$\begin{aligned} (\mu ^{-1}\nabla \times {\mathbf{u}}_h,\nabla \times {\mathbf{w}})_{\varOmega }&= ({\mathbf{j}},{\mathbf{w}})_{\varOmega }&\forall {\mathbf{w}}\in R_{k,0}(\mathcal {T}_h), \end{aligned}$$
(3a)
$$\begin{aligned} ({\mathbf{u}}_h,\nabla \psi )_{\varOmega }&=0,&\forall \psi \in P_{k,0}(\mathcal {T}_h). \end{aligned}$$
(3b)

The approximation of the magnetic field is then given by

$$\begin{aligned} {\mathbf{H}}_h:=\mu ^{-1}\nabla \times {\mathbf{u}}_h, \end{aligned}$$

which converges quasi-optimally as the mesh width h tends to zero [15, Theorem 5.10].

In the next section, we show how we can obtain a reliable and efficient estimator for \(\Vert {\mathbf{H}}-{\mathbf{H}}_h\Vert _{\varOmega }\).

3 An Equilibration-Based a Posteriori Error Estimator

We follow [6] and present an a posteriori error estimator that is based on the following result.

Theorem 3.1

( [6, Theorem 10]) Let \({\mathbf{u}}\) be the solution to (2), let \({\mathbf{u}}_h\) be the solution of (3), and set \({\mathbf{H}}:=\mu ^{-1}\nabla \times {\mathbf{u}}\) and \({\mathbf{H}}_h:=\mu ^{-1}\nabla \times {\mathbf{u}}_h\). If \({\tilde{{\mathbf{H}}}}\in H(\mathrm {curl};\varOmega )\) satisfies the equilibrium condition

$$\begin{aligned} \nabla \times {\tilde{{\mathbf{H}}}}&= {\mathbf{j}}, \end{aligned}$$
(4)

then

$$\begin{aligned} \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{\varOmega }&\le \Vert \mu ^{1/2}({\tilde{{\mathbf{H}}}}-{\mathbf{H}}_h)\Vert _{\varOmega }. \end{aligned}$$
(5)

Proof

The result follows from the orthogonality of \(\mu ^{1/2}({\tilde{{\mathbf{H}}}}-{\mathbf{H}})\) and \(\mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\):

$$\begin{aligned} \big (\mu ^{1/2}({\tilde{{\mathbf{H}}}}-{\mathbf{H}}),\mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\big )_{\varOmega }&=\big (\mu ^{1/2}({\tilde{{\mathbf{H}}}}-{\mathbf{H}}), \mu ^{-1/2}\nabla \times ({\mathbf{u}}-{\mathbf{u}}_h)\big )_{\varOmega } \\&= \big ({\tilde{{\mathbf{H}}}}-{\mathbf{H}},\nabla \times ({\mathbf{u}}-{\mathbf{u}}_h) \big )_{\varOmega } \\&= \big (\nabla \times ({\tilde{{\mathbf{H}}}}-{\mathbf{H}}),{\mathbf{u}}-{\mathbf{u}}_h \big )_{\varOmega } \\&= ({\mathbf{j}}-{\mathbf{j}}, {\mathbf{u}}-{\mathbf{u}}_h)_{\varOmega } \\&= 0 \end{aligned}$$

and Pythagoras’s theorem

$$\begin{aligned} \Vert \mu ^{1/2}({\tilde{{\mathbf{H}}}}-{\mathbf{H}}_h)\Vert ^2_{\varOmega }&= \Vert \mu ^{1/2}({\tilde{{\mathbf{H}}}}-{\mathbf{H}})\Vert ^2_{\varOmega } + \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert ^2_{\varOmega }. \end{aligned}$$
(6)

\(\square \)

Remark 3.2

Equation (6) is also known as a Prager–Synge type equation and obtaining an error estimator from such an equation is also known as the hypercircle method. Furthermore, equation (4) is known as the equilibrium condition and using the numerical approximation \({\mathbf{H}}_h\) to obtain a solution to this equation is called equilibration of \({\mathbf{H}}_h\).

Corollary 3.3

Let \({\mathbf{u}}\) be the solution to (2), let \({\mathbf{u}}_h\) be the solution of (3), set \({\mathbf{H}}:=\mu ^{-1}\nabla \times {\mathbf{u}}\) and \({\mathbf{H}}_h:=\mu ^{-1}\nabla \times {\mathbf{u}}_h\), and let \({\mathbf{j}}_h:=\nabla \times {\mathbf{H}}_h\) be the discrete current distribution. If \({\tilde{{\mathbf{H}}}}^{\varDelta }\in L^2(\varOmega )^3\) satisfies the (residual) equilibrium condition

$$\begin{aligned} \nabla \times {\tilde{{\mathbf{H}}}}^{\varDelta }&= {\mathbf{j}}-{\mathbf{j}}_h, \end{aligned}$$
(7)

which is an identity of distributions, then

$$\begin{aligned} \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{\varOmega }&\le \Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert _{\varOmega }. \end{aligned}$$
(8)

Proof

Since (7) is an identity of distributions, we can equivalently write

$$\begin{aligned} \langle \nabla \times {\tilde{{\mathbf{H}}}}^{\varDelta },{\mathbf{w}}\rangle&= \langle {\mathbf{j}}-{\mathbf{j}}_h, {\mathbf{w}}\rangle&\forall {\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3, \end{aligned}$$

where \(\langle \cdot ,\cdot \rangle \) denotes the application of a distribution to a function in \(\mathcal {C}_0^{\infty }(\varOmega )^3\). Now, set \({\tilde{{\mathbf{H}}}}:={\tilde{{\mathbf{H}}}}^{\varDelta }+{\mathbf{H}}_h \in L^2(\varOmega )^{3}\). Using the definition \({\mathbf{j}}_h:=\nabla \times {\mathbf{H}}_h\), we obtain

$$\begin{aligned} \langle \nabla \times {\tilde{{\mathbf{H}}}},{\mathbf{w}}\rangle&= \langle \nabla \times {\tilde{{\mathbf{H}}}}^{\varDelta }+ \nabla \times {\mathbf{H}}_h,{\mathbf{w}}\rangle&\\&= \langle {\mathbf{j}}-{\mathbf{j}}_h + \nabla \times {\mathbf{H}}_h, {\mathbf{w}}\rangle&\\&= \langle {\mathbf{j}}, {\mathbf{w}}\rangle&\forall {\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3. \end{aligned}$$

From this, it follows that \(\nabla \times {\tilde{{\mathbf{H}}}}={\mathbf{j}}\in L^2(\varOmega )^{3}\), so \({\tilde{{\mathbf{H}}}}\) is in \(H(\mathrm {curl};\varOmega )\) and satisfies equilibrium condition (4). Inequality (8) then follows from Theorem 3.1. \(\square \)

From Corollary 3.3, it follows that a constant-free upper bound on the error can be obtained from any field \({\tilde{{\mathbf{H}}}}^{\varDelta }\) that satisfies (7).

An error estimator of this type was first introduced in [6], where it is referred to as an equilibrated residual error estimator. There, \({\mathbf{j}}-{\mathbf{j}}_h\) is decomposed into a sum of local divergence-free current distributions \({\mathbf{j}}^{\varDelta }_i\) that have support on only a single vertex patch. The error estimator is then obtained by solving local problems of the form \(\nabla \times {\tilde{{\mathbf{H}}}}^{\varDelta }_i={\mathbf{j}}^{\varDelta }_i\) for each vertex patch and by then taking the sum of all local fields \({\tilde{{\mathbf{H}}}}^{\varDelta }_i\). It is, however, not straightforward to decompose \({\mathbf{j}}-{\mathbf{j}}_h\) into local divergence-free current distributions. An explicit expression for \({\mathbf{j}}^{\varDelta }_i\) is given in [6] for the lowest-degree Nédélec element, but this expression cannot be readily extended to basis functions of arbitrary degree.

Here, we instead present an error estimator based on equilibration condition (7) that can be applied to elements of arbitrary degree. Furthermore, instead of solving local problems on vertex patches, our estimator requires solving problems on only single elements, single faces, and small sets of nodes. The assumptions and a step-by-step derivation of the estimator are given in Sects. 3.1 and 3.2 below, a brief summary is given in Sect. 3.3, and reliability and efficiency are proven in Sect. 3.4.

3.1 Assumptions

In order to compute the error estimator, we use polynomial function spaces of degree \(k'\ge k\), where k denotes the degree of the finite element approximation \({\mathbf{u}}_h\), and assume that:

  1. A1.

    The magnetic permeability \(\mu \) is piecewise constant. In particular, the domain \(\varOmega \) can be partitioned into a finite set of polyhedral subdomains \(\{\varOmega _i\}\) such that \(\mu \) is constant on each subdomain \(\varOmega _i\). Furthermore, the mesh is assumed to be aligned with this partition so that \(\mu \) is constant within each element.

  2. A2.

    The current density \({\mathbf{j}}\) is in \(D_{k'}(\mathcal {T}_h)\cap H(\mathrm {div}^0;\varOmega )\).

Although assumption A2 does not hold in general, we can always replace \({\mathbf{j}}\) by a suitable projection \(\pi _h{\mathbf{j}}\) by taking, for example, \(\pi _h\) as the standard Raviart–Thomas interpolation operator corresponding to the \(D_{k'}(\mathcal {T}_h)\) space [16]. The error is in that case bounded by

$$\begin{aligned} \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert ^2_{\varOmega }&\le \Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert ^2_{\varOmega } + \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}') \Vert _{\varOmega }, \end{aligned}$$

where \({\mathbf{H}}':=\mu ^{-1}\nabla \times {\mathbf{u}}'\) and where \({\mathbf{u}}'\) is the solution to (2) with \({\mathbf{j}}\) replaced by \(\pi _h{\mathbf{j}}\). If \({\mathbf{j}}\) is sufficiently smooth, i.e. \({\mathbf{j}}\) can be extended to a function \({\mathbf{j}}^*\in H(\mathrm {div};{\mathbb {R}}^3)\cap H^{k'}({\mathbb {R}}^3)^3\) with compact support, then the term \(\Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}') \Vert _{\varOmega }\) is of order \(h^{k'+1}\); see Theorem A.1 in the appendix. This means that, if \(k'\ge k\), then \(\Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}') \Vert _{\varOmega }\) converges with a higher rate than \(\Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h) \Vert _{\varOmega }\) and so we may assume that the term \(\Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}') \Vert _{\varOmega }\) is negligible.

3.2 Derivation of the Error Estimator

Before we derive the error estimator, we first write \({\mathbf{j}}_h=\nabla \times {\mathbf{H}}_h\) in terms of element and face distributions. For every \({\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3\cup R_{k,0}(\mathcal {T}_h)\), we can write

$$\begin{aligned} \langle {\mathbf{j}}_h,{\mathbf{w}}\rangle&= \langle \nabla \times {\mathbf{H}}_h, {\mathbf{w}}\rangle \\&= ({\mathbf{H}}_h,\nabla \times {\mathbf{w}})_{\varOmega } \\&= \sum _{T\in \mathcal {T}_h} \left[ (\nabla \times {\mathbf{H}}_h,{\mathbf{w}})_T + ({\mathbf{H}}_h,\hat{{\mathbf{n}}}_T\times {\mathbf{w}})_{\partial T} \right] \\&= \sum _{T\in \mathcal {T}_h} \left[ (\nabla \times {\mathbf{H}}_h,{\mathbf{w}})_T + ({\mathbf{H}}_h,\hat{{\mathbf{n}}}_T\times {\mathbf{w}})_{\partial T\setminus \partial \varOmega } \right] \\&= \sum _{T\in \mathcal {T}_h} \left[ (\nabla \times {\mathbf{H}}_h,{\mathbf{w}})_T + (-\hat{{\mathbf{n}}}_T\times {\mathbf{H}}_h,{\mathbf{w}})_{\partial T\setminus \partial \varOmega } \right] \\&= \sum _{T\in \mathcal {T}_h} (\nabla \times {\mathbf{H}}_h,{\mathbf{w}})_T + \sum _{f\in \mathcal {F}_{h}^{I}} (-[\![{{\mathbf{H}}_h}]\!]_{t},{\mathbf{w}})_f \\&=: \sum _{T\in \mathcal {T}_h} ({\mathbf{j}}_{h,T},{\mathbf{w}})_T + \sum _{f\in \mathcal {F}_{h}^{I}} ({\mathbf{j}}_{h,f},{\mathbf{w}})_f, \end{aligned}$$

where \(\langle \cdot ,\cdot \rangle \) denotes the application of a distribution to a \(\mathcal {C}_0^{\infty }(\varOmega )^3\) function, \(\mathcal {F}_{h}^{I}\) denotes the set of all internal faces, \(\hat{{\mathbf{n}}}_{T}\) denotes the normal unit vector to \(\partial T\) pointing outward of T, and \([\![{{\mathbf{H}}}]\!]_{t}|_f:=(\hat{{\mathbf{n}}}^+\times {\mathbf{H}}^+ + \hat{{\mathbf{n}}}^-\times {\mathbf{H}}^-)|_{f}\) denotes the tangential jump operator, with \(\hat{{\mathbf{n}}}^{\pm }:=\hat{{\mathbf{n}}}_{T^{\pm }}\), \({\mathbf{H}}^{\pm }:={\mathbf{H}}|_{T^{\pm }}\), and \(T^+\) and \(T^-\) the two adjacent elements of f.

Since \({\mathbf{u}}_h\in R_{k,0}(\mathcal {T}_h)\) and \(\mu \) is piecewise constant, we have that \({\mathbf{H}}_h|_T\in P_{k-1}(T)^3\). Therefore, \({\mathbf{j}}_{h,T}\in P_{k-2}(T)^3\subset D_{k-1}(T)\) if \(k\ge 2\) and \({\mathbf{j}}_{h,T}=0\) if \(k=1\), and \({\mathbf{j}}_{h,f}\in \{{\mathbf{u}}\in P_{k-1}(f)^3 \;|\; \hat{{\mathbf{n}}}_f\cdot {\mathbf{u}}=0 \}\subset D_{k}(f)\), where \(\hat{{\mathbf{n}}}_f:=\hat{{\mathbf{n}}}^+|_f\) is a normal unit vector of f and \(D_{k}(f)\) is given by

$$\begin{aligned} D_{k}(f) := \{{\mathbf{u}}\in L^2(f)^{3} \;|\; {\mathbf{u}}({\mathbf{x}}) = \hat{{\mathbf{n}}}_f\times ({\mathbf{v}}({\mathbf{x}}) + {\mathbf{x}}w({\mathbf{x}})) \text { for some }{\mathbf{v}}\in P_{k-1}(f)^3, w\in P_{k-1}(f)\}. \end{aligned}$$

In other words, \({\mathbf{j}}_h\) can be represented by \(D_{k-1}(T)\) functions on the elements and \(D_{k}(f)\) face distributions on the internal faces.

We define \({\mathbf{j}}^{\varDelta }:={\mathbf{j}}-{\mathbf{j}}_h\) and can write

$$\begin{aligned} \langle {\mathbf{j}}^{\varDelta }, {\mathbf{w}}\rangle&= \sum _{T\in \mathcal {T}_h} ({\mathbf{j}}^{\varDelta }_{T},{\mathbf{w}})_T + \sum _{f\in \mathcal {F}_{h}^{I}} ({\mathbf{j}}^{\varDelta }_{f},{\mathbf{w}})_f&\forall {\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3 \cup R_{k,0}(\mathcal {T}_h), \end{aligned}$$
(9)

where

$$\begin{aligned} {\mathbf{j}}^{\varDelta }_{T} := {\mathbf{j}}|_T - {\mathbf{j}}_{h,T} = {\mathbf{j}}|_T -\nabla \times {\mathbf{H}}_h|_T \quad \text {and}\quad {\mathbf{j}}^{\varDelta }_{f} := -{\mathbf{j}}_{h,f} = [\![{{\mathbf{H}}_h}]\!]_{t}|_f. \end{aligned}$$
(10)

We look for a solution of (7) of the form \({\tilde{{\mathbf{H}}}}^{\varDelta }={\hat{{\mathbf{H}}}}^{\varDelta }+ \nabla _h\phi \), with \({\hat{{\mathbf{H}}}}^{\varDelta }\in R^{-1}_{k'}(\mathcal {T}_h)\) and \(\phi \in P^{-1}_{k'}(\mathcal {T}_h)\) and where \(\nabla _h\) denotes the element-wise gradient operator. The term \({\hat{{\mathbf{H}}}}^{\varDelta }\) will take care of the element distributions of \({\mathbf{j}}^{\varDelta }\) and the term \(\nabla _h\phi \) will take care of the remaining face distributions.

In the following, we firstly describe how to compute \({\hat{{\mathbf{H}}}}^{\varDelta }\) in Sect. 3.2.1 and characterize the remainder \({\mathbf{j}}^{\varDelta }- \nabla \times {\hat{{\mathbf{H}}}}^{\varDelta }\) in Sect. 3.2.2. We then describe how to compute the jumps of \(\phi \) on internal faces in Sect. 3.2.3 and explain how to reconstruct \(\phi \) from its jumps in Sect. 3.2.4.

3.2.1 Computation of \({\hat{\mathbf{H}}}^{\varDelta }\)

We compute \({\hat{{\mathbf{H}}}}^{\varDelta }\) by solving the local problems

$$\begin{aligned} \nabla \times {\hat{{\mathbf{H}}}}^{\varDelta }|_T&= {\mathbf{j}}^{\varDelta }_{T} ,&\end{aligned}$$
(11a)
$$\begin{aligned} (\mu {\hat{{\mathbf{H}}}}^{\varDelta },\nabla \psi )_T&= 0&\forall \psi \in P_{k'}(T), \end{aligned}$$
(11b)

for each element \(T\in \mathcal {T}_h\). This problem is well-defined and has a unique solution due to the discrete exact sequence property

and since \(\nabla \cdot {\mathbf{j}}^{\varDelta }_T =\nabla \cdot {\mathbf{j}}|_T - \nabla \cdot (\nabla \times {\mathbf{H}}_h)|_T= 0\) and \({\mathbf{j}}^{\varDelta }_T={\mathbf{j}}|_T -\nabla \times {\mathbf{H}}_h|_T\in D_{k'}(T)\). This last property follows from the fact that \({\mathbf{j}}|_T\in D_{k'}\) due to assumption A2 and \({\mathbf{H}}_h|_T=\mu ^{-1}\nabla \times {\mathbf{u}}_h|_T\in P_{k-1}(T)^3\subset P_{k'-1}(T)^{3}\) due to assumption A1.

3.2.2 Representation of the Remainder \({\mathbf{j}}^{\varDelta }- \nabla \times {\hat{{\mathbf{H}}}}^{\varDelta }\)

Set

For every \({\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3 \cup R_{k,0}(\mathcal {T}_h)\), we can write

so

(12)

for all \( {\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3 \cup R_{k,0}(\mathcal {T}_h)\). This means that can be represented by only face distributions, and since \({\mathbf{H}}_h\in P^{-1}_{k-1}(\mathcal {T}_h)^3\subset P^{-1}_{k'-1}(\mathcal {T}_h)^3\) and \({\hat{{\mathbf{H}}}}^{\varDelta }\in R^{-1}_{k'}(\mathcal {T}_h)\), we have that \([\![{{\mathbf{H}}_h + {\hat{{\mathbf{H}}}}^{\varDelta }}]\!]_{t}|_f \in D_{k'}(f)\) and therefore .

3.2.3 Computation of the Jumps of \(\phi \) on Internal Faces

It now remains to find a \(\phi \in P^{-1}_{k'}(\mathcal {T}_h)\) such that

For every \({\mathbf{w}}\in \mathcal {C}_0^{\infty }(\varOmega )^3 \cup R_{k,0}(\mathcal {T}_h)\), we can write

$$\begin{aligned} \langle \nabla \times \nabla _h\phi , {\mathbf{w}}\rangle&= (\nabla _h\phi , \nabla \times {\mathbf{w}})_{\varOmega } \\&= \sum _{T\in \mathcal {T}_h} \left[ (\nabla \times \nabla \phi ,{\mathbf{w}})_T + (\nabla \phi ,\hat{{\mathbf{n}}}_T\times {\mathbf{w}})_{\partial T} \right] \\&= \sum _{T\in \mathcal {T}_h} \left[ ({\mathbf{0}},{\mathbf{w}})_T + (\nabla \phi ,\hat{{\mathbf{n}}}_T\times {\mathbf{w}})_{\partial T\setminus \partial \varOmega } \right] \\&= \sum _{T\in \mathcal {T}_h} - (\hat{{\mathbf{n}}}_T\times \nabla \phi ,{\mathbf{w}})_{\partial T\setminus \partial \varOmega } \\&= \sum _{f\in \mathcal {F}_{h}^{I}} (-[\![{ \nabla \phi }]\!]_{t},{\mathbf{w}})_f. \end{aligned}$$

Therefore, we need to find a \(\phi \in P^{-1}_{k'}(\mathcal {T}_h)\) such that

To do this, we define, for each internal face \(f\in \mathcal {F}_{h}^{I}\), the scalar jump \([\![{\phi }]\!]_f:= (\phi ^+ - \phi ^-)|_{f}\) with \(\phi ^{\pm }:=\phi |_{T^{\pm }}\), two orthogonal unit tangent vectors \(\hat{{\mathbf{t}}}_1\) and \(\hat{{\mathbf{t}}}_2\) such that \(\hat{{\mathbf{t}}}_1\times \hat{{\mathbf{t}}}_2=\hat{{\mathbf{n}}}_f:=\hat{{\mathbf{n}}}^+|_f\), differential operators \(\partial _{t_i}:=\hat{{\mathbf{t}}}_i\cdot \nabla \), and the gradient operator restricted to the face: \(\nabla _f := \hat{{\mathbf{t}}}_1\partial _{t_1} + \hat{{\mathbf{t}}}_2\partial _{t_2}\). We can then write

$$\begin{aligned} -[\![{ \nabla \phi }]\!]_{t}|_f&= -(\hat{{\mathbf{n}}}^+\times \nabla \phi ^+ + \hat{{\mathbf{n}}}^-\times \nabla \phi ^-)|_{f} \\&= -(\hat{{\mathbf{n}}}_f\times \nabla \phi ^+ - \hat{{\mathbf{n}}}_f\times \nabla \phi ^-)|_{f} \\&= -(\hat{{\mathbf{n}}}_f\times \nabla _f\phi ^+ - \hat{{\mathbf{n}}}_f\times \nabla _f\phi ^-)|_{f} \\&= -\hat{{\mathbf{n}}}_f\times \nabla _f[\![{\phi }]\!]_f \end{aligned}$$

for all \(f\in \mathcal {F}_{h}^{I}\). We therefore introduce an auxiliary variable \(\lambda _f\in P_{k'}(f)\) and solve

(13a)
(13b)

for each \(f\in \mathcal {F}_{h}^{I}\), where (13b) is only added to ensure a unique solution. In the next section, we will show the existence of and how to construct a \(\phi \in P^{-1}_{k'}(\mathcal {T}_h)\) such that \([\![{\phi }]\!]_f=\lambda _f\) for all \(f\in \mathcal {F}_{h}^{I}\). Now, we will prove that problem (13) uniquely defines \(\lambda _f\). We start by showing that (13) corresponds to a 2D curl problem on a face. To see this, note that \(\hat{{\mathbf{n}}}_f\times \nabla _f = \hat{{\mathbf{t}}}_2\partial _{t,1} - \hat{{\mathbf{t}}}_1\partial _{t,2}\). If we take the inner product of (13a) with \(\hat{{\mathbf{t}}}_1\) and \(\hat{{\mathbf{t}}}_2\), we obtain

(14a)
(14b)

where , which is equivalent to a 2D curl problem on f. To show that (13) is well-posed, we use the discrete exact sequence in 2D:

where \(\mathrm {curl}_f:=(\partial _{t_2},-\partial _{t_1})\). Since , it suffices to show that . To prove this, we use that, for every \(\psi \in \mathcal {C}^\infty _0(\varOmega )^3\),

Then, for every \(\psi \in \mathcal {C}_0^{\infty }(\varOmega )\), we can write

where \(\mathcal {E}_{h}^{I}\) denotes the set of all internal edges, \(\hat{{\mathbf{n}}}_{\partial f}\) denotes the normal unit vector of \(\partial f\) that lies in the same plane as f and points outward of f, and \(\hat{{\mathbf{n}}}_{e,f}:=\hat{{\mathbf{n}}}_{\partial f}|_e\). This implies that

(15a)
(15b)

so for each internal face and, therefore, problem (13) is well-defined and has a unique solution.

3.2.4 Reconstruction of \(\phi \) from its Jumps on the Internal Faces

After computing \(\lambda _f\) for all internal faces, it remains to compute \(\phi \) such that \([\![{\phi }]\!]_f=\lambda _f\) for all \(f\in \mathcal {F}_{h}^{I}\). To do this, we use standard Lagrangian basis functions. For each element T, let \(\mathcal {Q}_T\) denote the set of nodes on element \({\overline{T}}\) for \(P_{k'}(T)\). The barycentric coordinates of these nodes are given by {\((\frac{i_1}{k'},\frac{i_2}{k'},\frac{i_3}{k'},\frac{i_4}{k'})\}_{i_1,i_2,i_3,i_4\ge 0, i_1+i_2+i_3+i_4=k'}\). Also, let \(\mathcal {Q}_h\) be the union of all element nodes. We then define the degrees of freedom for \(\phi \), denoted by \(\{\phi _{T,{\mathbf{x}}}\}_{T\in \mathcal {T}_h,{\mathbf{x}}\in \mathcal {Q}_T}\), as the values of \(\phi |_T\) at the nodes \({\mathbf{x}}\in \mathcal {Q}_T\). Since, for each \(f\in \mathcal {F}_{h}^{I}\), the space \(P_{k'}(f)\) is unisolvent on the nodes \(\mathcal {Q}_h\cap {\overline{f}}\) , we have that \([\![{\phi }]\!]_f=\lambda _f\) for all \(f\in \mathcal {F}_{h}^{I}\) if and only if

$$\begin{aligned} \phi _{T^+,{\mathbf{x}}} - \phi _{T^-,{\mathbf{x}}}&= \lambda _{f}({\mathbf{x}})&\forall f\in \mathcal {F}_{h}^{I}, \, {\forall }{\mathbf{x}}\in \mathcal {Q}_h\cap {\overline{f}}. \end{aligned}$$

We can decouple this global problem into very local problems. In particular, for each node \({\mathbf{x}}\in \mathcal {Q}_h\), we can compute the small set of degrees of freedom \(\{\phi _{T,{\mathbf{x}}}\}_{T:{\overline{T}}\ni {\mathbf{x}}}\) by solving

$$\begin{aligned} \phi _{T^+,{\mathbf{x}}} - \phi _{T^-,{\mathbf{x}}}&= \lambda _{f}({\mathbf{x}})&\forall f\in \mathcal {F}_{h}^{I}:{\overline{f}}\ni {\mathbf{x}}, \end{aligned}$$
(16a)
$$\begin{aligned} \sum _{T:{\overline{T}}\ni {\mathbf{x}}} \phi _{T,{\mathbf{x}}}&= 0.&\end{aligned}$$
(16b)

In this local problem, each degree of freedom corresponds to an element adjacent to \({\mathbf{x}}\) and for any two degrees of freedom corresponding to two adjacent elements, the difference should be equal to \(\lambda _f({\mathbf{x}})\), with f the face connecting the two elements. Condition (16b) is only added in order to ensure a unique local solution.

For a node \({\mathbf{x}}\) in the interior of an element, there is only one overlapping element T and the above results in \(\phi _{T,{\mathbf{x}}}=0\). The same applies to a node in the interior of a boundary face. For a node \({\mathbf{x}}\) in the interior of an internal face f, there are only two adjacent elements \(T^+\) and \(T^-\) and the above results in \(\phi _{T^+,{\mathbf{x}}}=\frac{1}{2}\lambda _f({\mathbf{x}})\) and \(\phi _{T^-,{\mathbf{x}}}=-\frac{1}{2}\lambda _f({\mathbf{x}})\). For a node in the interior of an edge, the degrees of freedom correspond to the ring (for internal edges) or the partial ring (for boundary edges) of elements adjacent to that edge. Finally, for a node on a vertex, the degrees of freedom correspond to the cloud of elements adjacent to that vertex.

For every cycle through elements adjacent to a node \({\mathbf{x}}\), the corresponding differences should add up to zero. A cycle means a sequence of elements \(T_1\rightarrow T_2\rightarrow \cdots \rightarrow T_{n+1}\) with \(n\ge 3\), such that \(T_1,T_2,\dots ,T_n\) are all different from each other, \(T_{n+1}=T_{1}\), and two consecutive elements are connected through a face. For a node in the interior of an internal edge, there is only one possible cycle, which is the cycle through the ring of elements adjacent to that edge. For a node on a vertex, the minimal cycles are the cycles around the internal edges adjacent to that vertex. Nodes in the interior of a face or in the interior of an element only have one or two adjacent elements and therefore have no cycles in their element patches.

Therefore, in general, the minimal cycles for any node \({\mathbf{x}}\) are the cycles around each internal edge connected to \({\mathbf{x}}\). To prove that the overdetermined system (16) is well-posed, it is therefore sufficient to check if, for each internal edge, the differences corresponding to the cycle around that edge sum to zero.

We can write this condition more formally. Let \(e\in \mathcal {E}_{h}^{I}\) be an internal edge, let \(\hat{{\mathbf{t}}}_e\) be a tangent unit vector of e, and let \(T_1\rightarrow T_2\rightarrow \cdots \rightarrow T_{n+1}\) be a cycle rotating counter-clockwise when looking towards \(\hat{{\mathbf{t}}}_e\). This means that the normal unit vector of the face \(T_i\cap T_{i+1}\) pointing out of \(T_{i+1}\) is given by \(\hat{{\mathbf{n}}}_{T_{i+1}}|_f=\hat{{\mathbf{t}}}_e\times \hat{{\mathbf{n}}}_{e,f}=:\hat{{\mathbf{n}}}_{f,e}\) (recall that \(\hat{{\mathbf{n}}}_{e,f}=\hat{{\mathbf{n}}}_{\partial f}|_e\)). Also, let \(\psi \in P^{-1}_{k'}(\mathcal {T}_h)\). The sum of the differences of \(\psi \) for this cycle can be written as

$$\begin{aligned} \sum _{i=1}^n (\psi _{T_{i+1}}-\psi _{T_i})|_e = (\psi _{T_{n+1}}-\psi _{T_1})|_e = (\psi _{T_{1}}-\psi _{T_1})|_e= 0. \end{aligned}$$

Now, let \(f=T_i\cap T_{i+1}\) and note that

$$\begin{aligned} (\psi _{T_{i+1}}-\psi _{T_i})|_f&= {\left\{ \begin{array}{ll} [\![{\psi }]\!]_f ,&{} \text {if } \hat{{\mathbf{n}}}_f=\hat{{\mathbf{n}}}_{T_{i+1}}|_f, \\ -[\![{\psi }]\!]_f ,&{} \text {if } \hat{{\mathbf{n}}}_f=-\hat{{\mathbf{n}}}_{T_{i+1}}|_f. \end{array}\right. } \end{aligned}$$

Therefore, we can write

$$\begin{aligned} (\psi _{T_{i+1}}-\psi _{T_i})|_f&= \big (\hat{{\mathbf{n}}}_f\cdot (\hat{{\mathbf{n}}}_{T_{i+1}}|_f)\big )[\![{\psi }]\!]_f = (\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e})[\![{\psi }]\!]_f. \end{aligned}$$

The sum of the differences of the cycle around e can therefore be rewritten as

$$\begin{aligned} 0=\sum _{i=1}^n (\psi _{T_{i+1}}-\psi _{T_i})|_e = \sum _{f:\partial f\supset e} (\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e})[\![{\psi }]\!]_f|_e \end{aligned}$$

and from this, we obtain the conditions

$$\begin{aligned} r_e:=\sum _{f:\partial f\supset e} (\hat{{\mathbf{n}}}_{f}\cdot \hat{{\mathbf{n}}}_{f,e})\lambda _{f}|_e&= 0&\forall e\in \mathcal {E}_{h}^{I}. \end{aligned}$$
(17)

Problem (16) is therefore well-posed provided that (17) is satisfied.

Fig. 1
figure 1

Illustration of \(\hat{{\mathbf{t}}}_e\), \(\hat{{\mathbf{n}}}_{e,f}\), and \(\hat{{\mathbf{n}}}_{f,e}\). The vector \(\hat{{\mathbf{n}}}_{f}\) either equals \(\hat{{\mathbf{n}}}_{f,e}\) or \(-\hat{{\mathbf{n}}}_{f,e}\)

Full size image

To prove (17), we first prove that, for each \(e\in \mathcal {E}_{h}^{I}\), \(r_e\) is constant and then prove that, for each \(e\in \mathcal {E}_{h}^{I}\), \((r_e,1)_e=0\).

Let e be an edge and let f be an adjacent face, and consider (14) with \(\hat{{\mathbf{t}}}_1=\hat{{\mathbf{t}}}_e\) and \(\hat{{\mathbf{t}}}_2=(\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e})\hat{{\mathbf{n}}}_{e,f}\). An illustration of \(\hat{{\mathbf{t}}}_e\), \(\hat{{\mathbf{n}}}_{e,f}\), and \(\hat{{\mathbf{n}}}_{f,e}\) is given in Fig.  1. The condition \(\hat{{\mathbf{n}}}_f=\hat{{\mathbf{t}}}_1\times \hat{{\mathbf{t}}}_2\) is still satisfied, since either \(\hat{{\mathbf{n}}}_f=\hat{{\mathbf{n}}}_{f,e}\) or \(\hat{{\mathbf{n}}}_f=-\hat{{\mathbf{n}}}_{f,e}\) and so

$$\begin{aligned} \hat{{\mathbf{t}}}_1\times \hat{{\mathbf{t}}}_2 = \hat{{\mathbf{t}}}_e\times \hat{{\mathbf{n}}}_{e,f}(\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e}) = \hat{{\mathbf{n}}}_{f,e}(\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e})= \hat{{\mathbf{n}}}_f. \end{aligned}$$

It then follows from (14b) that

where \(\partial _{t_e}:=\hat{{\mathbf{t}}}_e\cdot \nabla \). Multiplying the above by \(-(\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e})\) and summing over all faces adjacent to e results in

where the last line follows from (15a). We thus have

$$\begin{aligned} \partial _{t_e}r_e&= \sum _{f:\partial f\supset e} (\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e}) \partial _{t_e} \lambda _{f}|_e =0&\forall e\in \mathcal {E}_{h}^{I}. \end{aligned}$$
(18)

This implies that \(r_e\) is a constant. To prove (17), it therefore remains to show that \((r_e,1)_e=0\) for all internal edges.

To prove this, we define \(\varvec{\theta }_e\) to be the lowest-order Nédélec basis function corresponding to an internal edge e and scaled such that \((\hat{{\mathbf{t}}}_e\cdot \varvec{\theta }_e)|_e =1\) and derive

where the first two terms in the fifth line follow from the definition of \({\mathbf{u}}\) and \({\mathbf{u}}_h\) and the last term in the fifth line follows from (11b), assumption A1, and the fact that \(\nabla \times \varvec{\theta }_e\) is piecewise constant. We can then derive

where the second line follows from the fact that the tangent components of \(\varvec{\theta }_e\) are zero on all faces that are not adjacent to e, and the last line follows from (13b) and the fact that \(\nabla _f\times \varvec{\theta }_e\) is constant on f. We continue to obtain

$$\begin{aligned} 0&= \sum _{f:\partial f\supset e}(-\hat{{\mathbf{n}}}_f\times \hat{{\mathbf{n}}}_{\partial f}\lambda _f,\varvec{\theta }_e)_{\partial f} \\&= \sum _{f:\partial f\supset e} \sum _{\tilde{e}:\tilde{e}\subset \partial {f}} (-\hat{{\mathbf{n}}}_f\times \hat{{\mathbf{n}}}_{\tilde{e},f} \lambda _f,\varvec{\theta }_e)_{\tilde{e}} \\&= \sum _{f:\partial f\supset e} \sum _{\tilde{e}:\tilde{e}\subset \partial {f}} \big (-(\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,{\tilde{e}}})\hat{{\mathbf{n}}}_{f,{\tilde{e}}}\times \hat{{\mathbf{n}}}_{{\tilde{e}},f} \lambda _f,\varvec{\theta }_e \big )_{\tilde{e}} \\&= \sum _{f:\partial f\supset e} \sum _{\tilde{e}:\tilde{e}\subset \partial {f}} \big (-(\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,{\tilde{e}}})(\hat{{\mathbf{t}}}_{{\tilde{e}}}\times \hat{{\mathbf{n}}}_{{\tilde{e}},f})\times \hat{{\mathbf{n}}}_{{\tilde{e}},f} \lambda _f,\varvec{\theta }_e \big )_{\tilde{e}} \\&= \sum _{f:\partial f\supset e}\sum _{\tilde{e}:\tilde{e}\subset \partial {f}} \big ( (\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,{\tilde{e}}}) \hat{{\mathbf{t}}}_{\tilde{e}}\lambda _f,\varvec{\theta }_e \big )_{\tilde{e}} \\&=\sum _{f:\partial f\supset e}\sum _{\tilde{e}:\tilde{e}\subset \partial {f}} \big ((\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,{\tilde{e}}})\lambda _f,\hat{{\mathbf{t}}}_{\tilde{e}}\cdot \varvec{\theta }_e \big )_{\tilde{e}} \\&=\sum _{f:\partial f\supset e} \big ((\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e})\lambda _f,1 \big )_{e}, \end{aligned}$$

where the third line follows from the fact that either \(\hat{{\mathbf{n}}}_f=\hat{{\mathbf{n}}}_{f,{\tilde{e}}}\) or \(\hat{{\mathbf{n}}}_f=-\hat{{\mathbf{n}}}_{f,{\tilde{e}}}\), the fifth line follows from the property \(({\mathbf{a}}\times {\mathbf{b}})\times {\mathbf{c}}=-{\mathbf{a}}({\mathbf{b}}\cdot {\mathbf{c}})+{\mathbf{b}}({\mathbf{a}}\cdot {\mathbf{c}})\) and the fact that \(\hat{{\mathbf{t}}}_{{\tilde{e}}}\) and \(\hat{{\mathbf{n}}}_{{\tilde{e}},f}\) are orthogonal, and the last line follows from the fact that \((\hat{{\mathbf{t}}}_{{\tilde{e}}}\cdot \varvec{\theta }_e)|_{{\tilde{e}}}=0\) for all edges \(\tilde{e}\ne e\) and from \((\hat{{\mathbf{t}}}_e\cdot \varvec{\theta }_e)|_e=1\). We then continue to obtain

$$\begin{aligned} 0&= \sum _{f:\partial f\supset e} \big ( (\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e}) \lambda _f,1 \big )_{e} \\&= \left( \sum _{f:\partial f\supset e} (\hat{{\mathbf{n}}}_f\cdot \hat{{\mathbf{n}}}_{f,e}) \lambda _f, 1 \right) _{e} \\&= (r_e,1)_{e}. \end{aligned}$$

Therefore, \((r_e,1)_{e}=0\) and hence \(r_e=0\) for all internal edges and so problem (16) is well-posed.

3.3 Summary of Computing the Error Estimator

We fix \(k'\ge k\) and compute our error estimator in four steps.

Step 1. We compute \({\hat{{\mathbf{H}}}}^{\varDelta }\in R^{-1}_{k'}(\mathcal {T}_h)\) from the datum \({\mathbf{j}}\) and the numerical solution \({\mathbf{H}}_h\) by solving

$$\begin{aligned} \nabla \times {\hat{{\mathbf{H}}}}^{\varDelta }|_T&= {\mathbf{j}}^{\varDelta }_T = {\mathbf{j}}|_T-\nabla \times {\mathbf{H}}_h|_T,&\end{aligned}$$
(19a)
$$\begin{aligned} (\mu {\hat{{\mathbf{H}}}}^{\varDelta },\nabla \psi )_T&= 0&\forall \psi \in P_{k'}(T), \end{aligned}$$
(19b)

for each \(T\in \mathcal {T}_h\).

Step 2. For each internal face \(f\in \mathcal {F}_{h}^{I}\), let \(T^+\) and \(T^-\) denote the two adjacent elements, let \(\hat{{\mathbf{n}}}^{\pm }\) denote the normal unit vector pointing outward of \(T^{\pm }\), let \({\mathbf{H}}^{\pm }:={\mathbf{H}}|_{T^{\pm }}\) denote the vector field restricted to \(T^{\pm }\), let \([\![{{\mathbf{H}}}]\!]_{t}|_f := (\hat{{\mathbf{n}}}^+\times {\mathbf{H}}^+ + \hat{{\mathbf{n}}}^-\times {\mathbf{H}}^-)|_{f}\) denote the tangential jump operator, and let \(\nabla _f\) denote the gradient operator restricted to face f. We set \(\hat{{\mathbf{n}}}_f:=\hat{{\mathbf{n}}}^+|_f\) and compute \(\lambda _f\in P_{k'}(f)\) by solving

(20a)
(20b)

for each internal face \(f\in \mathcal {F}_{h}^{I}\).

Step 3. We compute \(\phi \in P^{-1}_{k'}(\mathcal {T}_h)\) by solving, for each \({\mathbf{x}}\in \mathcal {Q}_h\), the small set of degrees of freedom \(\{\phi _{T,{\mathbf{x}}}\}_{T:{\overline{T}}\ni {\mathbf{x}}}\) such that

$$\begin{aligned} \phi _{T^+,{\mathbf{x}}} - \phi _{T^-,{\mathbf{x}}}&= \lambda _{f}({\mathbf{x}})&\forall f\in \mathcal {F}_{h}^{I}:\overline{f}\ni {\mathbf{x}}, \end{aligned}$$
(21a)
$$\begin{aligned} \sum _{T:{\overline{T}}\ni {\mathbf{x}}} \phi _{T,{\mathbf{x}}}&= 0,&\end{aligned}$$
(21b)

where \(\mathcal {Q}_h\) denotes the set of standard Lagrangian nodes corresponding to the finite element space \(P_{k'}(\mathcal {T}_h)\) and \(\phi _{T,{\mathbf{x}}}\) denotes the value of \(\phi |_T\) at node \({\mathbf{x}}\).

Step 4. We compute the field

$$\begin{aligned} {\tilde{{\mathbf{H}}}}^{\varDelta }= {\hat{{\mathbf{H}}}}^{\varDelta }+ \nabla _h\phi , \end{aligned}$$

where \(\nabla _h\) denotes the element-wise gradient operator, and compute the error estimator

$$\begin{aligned} \eta _h :=\Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert _{\varOmega } = \left( \sum _{T\in \mathcal {T}_h}\eta _T^2\right) ^{1/2},&\quad \text {where}\quad \eta _T :=\Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert _T. \end{aligned}$$
(22)

Remark 3.4

For the case \(k^\prime =k=1\), this algorithm requires solving local problems that involve 6 unknowns per element (Step 1), 3 unknowns per face (Step 2), and \((\# T:{\overline{T}}\ni \nu ) \approx 24\) unknowns per vertex \(\nu \) (Step 3). Hence, the third step is expected to be the most computationally expensive step, although it still involves only 3 times as few unknowns as the vertex-patch problems of [6].

3.4 Reliability and Efficiency

The results of the previous sections immediately give the following theorem.

Theorem 3.5

(reliability) Let \({\mathbf{u}}\) be the solution to (2), let \({\mathbf{u}}_h\) be the solution to (3), and set \({\mathbf{H}}:=\mu ^{-1}\nabla \times {\mathbf{u}}\) and \({\mathbf{H}}_h:=\mu ^{-1}\nabla \times {\mathbf{u}}_h\). Also, fix \(k'\ge k\) and assume that assumptions A1 and A2 hold true. Then the problems in Steps 1–3 are all well-defined and have a unique solution. Furthermore, if \({\tilde{{\mathbf{H}}}}^{\varDelta }\) is computed by following Steps 1–4 in Sect. 3.3, then

$$\begin{aligned} \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{\varOmega }&\le \Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert _{\varOmega }=\eta _h. \end{aligned}$$
(23)

We also prove local efficiency of the estimator and state it as the following theorem.

Theorem 3.6

(local efficiency) Let \({\mathbf{u}}\) be the solution to (2), let \({\mathbf{u}}_h\) be the solution to (3), and set \({\mathbf{H}}:=\mu ^{-1}\nabla \times {\mathbf{u}}\) and \({\mathbf{H}}_h:=\mu ^{-1}\nabla \times {\mathbf{u}}_h\). Also, fix \(k'\ge k\) and assume that assumptions A1 and A2 hold true. If \({\tilde{{\mathbf{H}}}}^{\varDelta }\) is computed by following Steps 1–4 in Sect. 3.3, then

$$\begin{aligned} \eta _T&= \Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert _T \le C\sum _{T':\overline{T'}\cap {\overline{T}}\ne \emptyset } \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{T'} \end{aligned}$$
(24)

for all \(T\in \mathcal {T}_h\), where C is some positive constant that depends on the magnetic permeability \(\mu \), the shape-regularity of the mesh and the polynomial degree \(k'\), but not on the mesh width h.

Proof

In this proof, we always let C denote some positive constant that does not depend on the mesh width h, but may depend on the magnetic permeability \(\mu \), the shape-regularity of the mesh and the polynomial degree \(k'\).

1. Fix \(T\in \mathcal {T}_h\) and let \({\hat{{\mathbf{H}}}}^{\varDelta }_T\in R_{k'}(T)\) be the solution to (19). To obtain an upper bound for \(\Vert \mu ^{1/2}{\hat{{\mathbf{H}}}}^{\varDelta }_T\Vert _T\), we consider (19) for a reference element. Let \({\hat{T}}\) denote the reference tetrahedron, let \(\varvec{\varphi }_T:{\hat{T}}\rightarrow T\) denote the affine element mapping, and let \(J_T := [\frac{\partial \varvec{\varphi }_T}{\partial {\hat{x}}_1}\, \frac{\partial \varvec{\varphi }_T}{\partial {\hat{x}}_2}\, \frac{\partial \varvec{\varphi }_T}{\partial {\hat{x}}_3}]\) be the Jacobian of \(\varvec{\varphi }_T\). Also, let \({\hat{{\mathbf{G}}}}^{\varDelta }_{{{\hat{T}}}}\in R_{k'}({{\hat{T}}})\) be the solution of

$$\begin{aligned} \nabla \times {\hat{{\mathbf{G}}}}^{\varDelta }_{{{\hat{T}}}}&= {\mathbf{j}}^{\varDelta }_{{{\hat{T}}}}&\text {in }{\hat{T}}, \end{aligned}$$
(25a)
$$\begin{aligned} ({\hat{{\mathbf{G}}}}^{\varDelta }_{{{\hat{T}}}},\nabla \psi )_{{{\hat{T}}}}&= 0&\forall \psi \in P_{k'}({{\hat{T}}}), \end{aligned}$$
(25b)

where \({\mathbf{j}}^{\varDelta }_{{{\hat{T}}}}\) is defined as the pull-back of \({\mathbf{j}}^{\varDelta }_{T}\) through the Piola transformation, namely such that \({\mathbf{j}}^{\varDelta }_T\circ \varvec{\varphi }_T= \frac{1}{\mathrm {det}(J_T)}J_T {\mathbf{j}}^{\varDelta }_{{{\hat{T}}}}\), with \(\mathrm {det}(J_T)\) the determinant of \(J_T\). From the discrete Friedrichs inequality, it follows that

$$\begin{aligned} \Vert {\hat{{\mathbf{G}}}}^{\varDelta }_{{{\hat{T}}}}\Vert _{{{\hat{T}}}} \le C \Vert {\mathbf{j}}^{\varDelta }_{{{\hat{T}}}} \Vert _{{{\hat{T}}}}. \end{aligned}$$
(26)

Furthermore, if we set \({\hat{{\mathbf{G}}}}^{\varDelta }_T\) as the push-forward of \({\hat{{\mathbf{G}}}}^{\varDelta }_{{{\hat{T}}}}\) through the covariant transformation, namely \({\hat{{\mathbf{G}}}}^{\varDelta }_T\circ \varvec{\varphi }_T:=J_T^{-t} {\hat{{\mathbf{G}}}}^{\varDelta }_{{{\hat{T}}}}\), where \(J_T^{-t}\) is the inverse of the transpose of the Jacobian \(J_T\), then it can be checked that

$$\begin{aligned} \nabla \times {\hat{{\mathbf{G}}}}^{\varDelta }_T&= {\mathbf{j}}^{\varDelta }_T&\text {in }T. \end{aligned}$$

From (26), it also follows that

$$\begin{aligned} \Vert \mu ^{1/2}{\hat{{\mathbf{G}}}}^{\varDelta }_T\Vert _T&\le Ch_T \Vert {\mathbf{j}}^{\varDelta }_T\Vert _T, \end{aligned}$$
(27)

where \(h_T\) denotes the diameter of T. Now, since \({\hat{{\mathbf{G}}}}^{\varDelta }_T,{\hat{{\mathbf{H}}}}^{\varDelta }_T\in R_{k'}(T)\) and since \(\nabla \times ({\hat{{\mathbf{H}}}}^{\varDelta }_T-{\hat{{\mathbf{G}}}}^{\varDelta }_T)={\mathbf{j}}^{\varDelta }_T-{\mathbf{j}}^{\varDelta }_T={\mathbf{0}}\), we have that \({\hat{{\mathbf{H}}}}^{\varDelta }_T-{\hat{{\mathbf{G}}}}^{\varDelta }_T = \nabla \psi \) for some \(\psi \in P_{k'}(T)\). From (19b) and Pythagoras’s theorem, it then follows that

$$\begin{aligned} \Vert \mu ^{1/2} {\hat{{\mathbf{G}}}}^{\varDelta }_T\Vert _T^2&= \Vert \mu ^{1/2} {\hat{{\mathbf{H}}}}^{\varDelta }_T\Vert _T^2 + \Vert \mu ^{1/2} ({\hat{{\mathbf{G}}}}^{\varDelta }_T-{\hat{{\mathbf{H}}}}^{\varDelta }_T)\Vert _T^2\ge \Vert \mu ^{1/2} {\hat{{\mathbf{H}}}}^{\varDelta }_T\Vert _T, \end{aligned}$$

and from (27), it then follows that

$$\begin{aligned} \Vert \mu ^{1/2} {\hat{{\mathbf{H}}}}^{\varDelta }_T\Vert _T&\le Ch_T\Vert {\mathbf{j}}^{\varDelta }_T\Vert _T. \end{aligned}$$
(28)

2. Fix \(f\in \mathcal {F}_{h}^{I}\) and let \(\lambda _f\) be the solution to (20). In a way, similarly as for \({\hat{{\mathbf{H}}}}^{\varDelta }_T\), we can obtain the bound

(29)

where \(h_f\) denotes the diameter of the face f. Since , we can use (29), the triangle inequality, the trace inequality, discrete inverse inequality, and (28) to obtain

$$\begin{aligned} \Vert \lambda _f\Vert _f&\le Ch_f\left( \Vert {\mathbf{j}}^{\varDelta }_f\Vert _f + \Vert [\![{{\hat{{\mathbf{H}}}}^{\varDelta }}]\!]_{t} \Vert _f \right) \nonumber \\&\le Ch_f\left( \Vert {\mathbf{j}}^{\varDelta }_f\Vert _f + h_{T^+}^{-1/2}\Vert {\hat{{\mathbf{H}}}}^{\varDelta }_{T^+} \Vert _{T^+} + h_{T^-}^{-1/2}\Vert {\hat{{\mathbf{H}}}}^{\varDelta }_{T^-}\Vert _{T^-} \right) \nonumber \\&\le C\left( h_f\Vert {\mathbf{j}}^{\varDelta }_f\Vert _f + h_{T^+}^{3/2}\Vert {\mathbf{j}}^{\varDelta }_{T^+}\Vert _{T^+} + h_{T^-}^{3/2}\Vert {\mathbf{j}}^{\varDelta }_{T^-}\Vert _{T^-} \right) , \end{aligned}$$
(30)

where \(T^+\) and \(T^-\) are the two adjacent elements of f.

3. Fix \({\mathbf{x}}\in \mathcal {Q}_h\) and let \(\{\phi _{T,{\mathbf{x}}}\}_{T:{\overline{T}}\ni {\mathbf{x}}}\) be the solution to (21). Since \(\{\phi _{T,{\mathbf{x}}}\}_{T:{\overline{T}}\ni {\mathbf{x}}}\) depends linearly on \(\{\lambda _f({\mathbf{x}})\}_{f:{\overline{f}}\ni {\mathbf{x}}}\), and since the number of possible nodal patches, which depends on the mesh regularity, is finite, we can obtain the bound

$$\begin{aligned} \left( \sum _{T:{\overline{T}}\ni {\mathbf{x}}} \phi _{T,{\mathbf{x}}}^2\right) ^{1/2} \le C\left( \sum _{f:{\overline{f}}\ni {\mathbf{x}}} \lambda _{f}({\mathbf{x}})^2 \right) ^{1/2}. \end{aligned}$$

Now, fix \(T\in \mathcal {T}_h\). From the above and (30), we can obtain

$$\begin{aligned} \Vert \phi \Vert _T&\le C \sum _{f:{\overline{f}}\cap {\overline{T}}\ne \emptyset } h_f^{1/2} \Vert \lambda _f\Vert _f \\&\le C \sum _{f:{\overline{f}}\cap {\overline{T}}\ne \emptyset } \left[ h_f^{3/2}\Vert {\mathbf{j}}^{\varDelta }_f\Vert _f + h_{T^+}^{2}\Vert {\mathbf{j}}^{\varDelta }_{T^+}\Vert _{T^+} + h_{T^-}^{2}\Vert {\mathbf{j}}^{\varDelta }_{T^-}\Vert _{T^-} \right] . \end{aligned}$$

Note that \(h_f,h_{T'}\le Ch_T\) for all \(f:{\overline{f}}\cap T\ne \emptyset \) and \(T':\overline{T'}\cap {\overline{T}}\ne \emptyset \), due to the regularity of the mesh, which is incorporated in the constant C. Using the above and the discrete inverse inequality, we then obtain

$$\begin{aligned} \Vert \mu ^{1/2}\nabla \phi \Vert _T&\le Ch_T^{-1}\Vert \phi \Vert _T \nonumber \\&\le C\sum _{f:{\overline{f}}\cap {\overline{T}}\ne \emptyset } \left[ h_f^{1/2}\Vert {\mathbf{j}}^{\varDelta }_f\Vert _f + h_{T^+}\Vert {\mathbf{j}}^{\varDelta }_{T^+}\Vert _{T^+} + h_{T^-}\Vert {\mathbf{j}}^{\varDelta }_{T^-}\Vert _{T^-} \right] . \end{aligned}$$
(31)

4. We now use the efficiency estimate of the residual error estimator established in [4]. This estimate can be written as

$$\begin{aligned} h_T\Vert {\mathbf{j}}^{\varDelta }_T\Vert _T&= h_T\Vert {\mathbf{j}}- \nabla \times {\mathbf{H}}_h \Vert _T \le C\Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _T \\ h_f^{1/2}\Vert {\mathbf{j}}^{\varDelta }_f\Vert _f&= h_f^{1/2}\Vert [\![{{\mathbf{H}}_h}]\!]_{t}\Vert _f \le C\left( \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{T^+} + \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{T^-}\right) \end{aligned}$$

for all \(T\in \mathcal {T}_h\), \(f\in \mathcal {F}_{h}^{I}\). Using that \({\tilde{{\mathbf{H}}}}^{\varDelta }|_T = {\hat{{\mathbf{H}}}}^{\varDelta }_T + \nabla \phi |_T\), the triangle inequality, (28), (31), and the above, we obtain

$$\begin{aligned} \Vert \mu ^{1/2}{\tilde{{\mathbf{H}}}}^{\varDelta }\Vert _T&\le \Vert \mu ^{1/2}{\hat{{\mathbf{H}}}}^{\varDelta }\Vert _{T} + \Vert \mu ^{1/2}\nabla \phi \Vert _T \\&\le Ch_T\Vert {\mathbf{j}}^{\varDelta }\Vert _T + C\sum _{f:{\overline{f}}\cap {\overline{T}}\ne \emptyset } \left[ h_f^{1/2}\Vert {\mathbf{j}}^{\varDelta }_f\Vert _f + h_{T^+}\Vert {\mathbf{j}}^{\varDelta }_{T^+}\Vert _{T^+} + h_{T^-}\Vert {\mathbf{j}}^{\varDelta }_{T^-}\Vert _{T^-} \right] \\&\le C\!\!\!\sum _{T':\overline{T'}\cap {\overline{T}}\ne \emptyset } \Vert \mu ^{1/2}({\mathbf{H}}-{\mathbf{H}}_h)\Vert _{T'}, \end{aligned}$$

which completes the proof. \(\square \)

Remark 3.7

In Theorem 3.6, we proved that the constant in the efficiency estimate is independent of the mesh resolution, but may depend on \(k^\prime \), as we used discrete inverse inequalities and the efficiency result stated in [4]. Whether it is also independent of \(k^\prime \) is an open issue.

4 Numerical Experiments

In this section, we present several numerical results for the unit cube and the L-brick domain with constant magnetic permeability \(\mu =1\) for the a posteriori estimator constructed according to Steps 1–4 in Sect. 3.3.

For efficiency of the computations, we choose the same polynomial degree for the computation of the a posteriori error estimator as for the approximation \({\mathbf{u}}_h\), i.e. \(k'=k\). In the numerical experiments, we do not project the right hand side \({\mathbf{j}}\) onto \(D_k(\mathcal {T}_h)\), but solve the local problems of Step 1 in variational form. This introduces small compatibility errors into Step 2 and Step 3, which we observe to be negligible for the computation of the error estimator, cf. also the previous discussion on assumption A2 in Sect. 3.1. The orthogonality conditions of the local problems of Step 1 and Step 2 are incorporated via Lagrange multipliers and the resulting local saddle point problems are solved with a direct solver. For the computation of \(\phi \) in Step 3, we solve the local overdetermined systems with the least-squares method. Since we have shown that the solutions to those discrete problems are unique, the least-squares method computes those discrete unique solutions.

Instead of solving the full discrete problem (3), the numerical approximation \({\mathbf{u}}_h\) is computed by solving the singular system corresponding to (3a) only, since the gauge condition (3b) does not affect the variable of interest \({\mathbf{H}}_h:=\mu ^{-1}\nabla \times {\mathbf{u}}\). In order to do this, we use the preconditioned conjugate gradient algorithm in combination with a multigrid preconditioner [2, 14]. To ensure that, in the presence of quadrature errors, the discretised right-hand side remains discretely divergence free, a small gradient correction is added following [9, Sect. 4.1]. Our implementation of \(R_k(\mathcal {T}_h)\) is based on the hierarchical basis functions from [22].

We evaluate the reliability and efficiency of the a posteriori error estimator \(\eta _h\) defined in (22) for uniform and for adaptive meshes. For adaptive mesh refinement, we employ the standard adaptive finite element algorithm. Firstly we solve the discrete problem (3), then we compute the a posteriori error estimator to estimate the error; based on the local values \(\eta _T\) of the estimator, we mark elements for refinement based on the bulk marking strategy [12] with bulk parameter \(\theta =0.5\), and finally we refine the marked elements based on a bisection strategy [3].

Fig. 2
figure 2

Error and efficiency indices for the unit cube example with polynomial solution and uniformly refined meshes.

Full size image

4.1 Unit Cube Example

For the first example we chose the unit cube \(\varOmega =(0,1)^3\) with homogeneous boundary conditions and right hand side \({\mathbf{j}}\) according to the polynomial solution

$$\begin{aligned} {\mathbf{u}}(x,y,z) = \left( \begin{array}{c} y(1-y)z(1-z)\\ x(1-x)z(1-z) \\ x(1-x)y(1-y) \end{array}\right) . \end{aligned}$$

In Fig. 2, we present the errors \(\Vert {\mathbf{H}}-{\mathbf{H}}_h \Vert _{\varOmega }\) and efficiency indices \(\eta _h/\Vert {\mathbf{H}}-{\mathbf{H}}_h \Vert _{\varOmega }\) for \(k=1,2,3\), on a sequence of uniformly refined meshes. We observe that the error converges with optimal rates \(\mathcal {O}(h^k) = \mathcal {O}(N_h^{-k/3})\), \(N_h = {\text {dim}}(R_k(\mathcal {T}_h))\), and that the error estimator is reliable and efficient with efficiency constants between 1 and 2.

Note that, for \(k=3\), \({\mathbf{j}}\) belongs to \(D_k(\mathcal {T}_h)\), hence in that case assumption A2 is valid. We observe that, in the other cases \(k=1,2\), the estimator is reliable and efficient as well, even though \({\mathbf{j}}\) does not belong to \(D_k(\mathcal {T}_h)\). Thus the error introduced in the computation of the error estimator by not satisfying A2 is negligible.

Fig. 3
figure 3

Error and efficiency indices for adaptive mesh refinement for the L-brick example

Full size image
Fig. 4
figure 4

Polynomial robustness of the equilibrated a posteriori error estimator \(\eta _h\) in comparison to the residual a posteriori error estimator \(\mu _h\) for the L-brick example

Full size image

4.2 L-Brick Example

As second example, we solve the homogeneous Maxwell problem on the (nonconvex) domain

$$\begin{aligned} \varOmega = (-1,1)\times (-1,1)\times (0,1)\setminus \left( [0,1]\times [-1,0]\times [0,1]\right) . \end{aligned}$$

As solution, we choose the singular function

$$\begin{aligned} {\mathbf{u}}(x,y,z) = \nabla \times \left( \begin{array}{c} 0\\ 0 \\ (1-x^2)^2(1-y^2)^2((1-z)z)^2r^{2/3}\cos (\frac{2}{3}\varphi ) \end{array}\right) , \end{aligned}$$

where \((r,\varphi )\) are the two dimensional polar coordinates in the x-y-plane, and we choose the right hand side \({\mathbf{j}}\) accordingly.

Due to the edge singularity, uniform mesh refinement leads to suboptimal convergence rates of \(\mathcal {O}(N_h^{-2/9})\), as shown in Fig. 3, left plot, for \(k=2\). In contrast, adaptive mesh refinement leads to faster convergence rates. Note that, for this example, due to the edge singularity, anisotropic adaptive mesh refinement is needed to observe optimal convergence rates. Hence, the convergence rates in Fig. 3 are limited by the employed isotropic adaptive mesh refinement. In Fig. 3, we observe the best possible rates for adaptive isotropic mesh refinement that is \(\mathcal {O}(N_h^{-1/3})\) for \(k=1\), \(\mathcal {O}((N_h/\ln (N_h))^{-2/3})\) for \(k=2\), and \(\mathcal {O}(N_h^{-2/3})\) for \(k\ge 3\), cf. [1, Sect. 4.2.3]. This shows experimentally that the adaptive algorithm generates meshes which are quasi-optimal for isotropic refinements. Again, the (global) efficiency indices are approximately between 1 and 2, as shown in Fig. 3, right plot.

Next, we compare the efficiency indices of \(\eta _h\) for k-refinement to those of the residual a posteriori error estimator [4]

$$\begin{aligned} \mu _h^2 := \sum _{T\in \mathcal {T}_h} \frac{h_T^2}{k^2} \Vert {\mathbf{j}}- \nabla \times {\mathbf{H}}_h \Vert _{T}^2 + \sum _{f\in \mathcal {F}_{h}^{I}}\frac{h_f}{k} \Vert [\![{ {\mathbf{H}}_h }]\!]_{t} \Vert _{f}^2, \end{aligned}$$

with \(h_T\) the diameter of element T and \(h_f\) the diameter of face f. In Fig. 4, we observe the well known fact that the residual a posteriori error estimator \(\mu _h\) is not robust in the polynomial degree k, whereas our new estimator appears to be more robust in k.

4.3 Example with Discontinuous Permeability

In this example, we choose a discontinuous permeability

$$\begin{aligned} \mu (x,y,z) = \left\{ \begin{array}{ll} \mu _1 &{} \text {if } y<1/2 \text { and } z < 1/2,\\ \mu _2 &{} \text {otherwise}, \end{array} \right. \end{aligned}$$

on the unit cube \(\varOmega = (0,1)^3\) and the right hand side \({\mathbf {j}}=(1,0,0)^t\). We choose \(\mu _1=1\) and vary \(\mu _2= 10^\ell \) for \(\ell =1,2,3\). Since the exact solution is unknown, we approximate the error by comparing the numerical approximations to a reference solution, which is obtained from the last numerical approximation by 8 more adaptive mesh refinements. In these numerical experiments, we restrict ourselves to \(k=2\). As shown in Fig. 5, adaptive mesh refinement leads to optimal convergence rates for isotropic mesh refinement (\(\mathcal {O}((N_h/\ln (N_h))^{-2/3})\)), and the efficiency indices are between 1 and 2, independently of the contrast of the discontinuous permeability. In Fig. 6, we display an adaptive mesh after 14 refinement steps with about \(5\cdot 10^4\) degrees of freedom. We observe strong adaptive mesh refinement towards the edge between the points \((0,1/2,1/2)^t\) and \((1,1/2,1/2)^t\).

Fig. 5
figure 5

Error and efficiency indices for adaptive mesh refinement for the example with discontinuous permeability

Full size image
Fig. 6
figure 6

Adaptive surface meshes of \(\varOmega \) (left) and of the subdomain of the coefficient \(\mu _2=1000\) (right)

Full size image

5 Conclusion

We have presented a novel a posteriori error estimator for arbitrary-degree Nédélec elements for solving magnetostatic problems. This estimator is based on an equilibration principle and is obtained by solving only very local problems (on single elements, on single faces, and on very small sets of nodes). We have derived a constant-free reliability estimate and a local efficiency estimate, and presented numerical tests, involving a smooth solution and a singular solution, that confirm these results. Moreover, the numerical results show an efficiency index between 1 and 2 in all considered cases, also for large polynomial degrees k, and the dependence on k appears to be small. Some remaining questions are how to extend the proposed error estimator for domains with curved boundaries or for domains with a smoothly varying permeability.