An Equilibrated a Posteriori Error Estimator for Arbitrary-Order Nédélec Elements for Magnetostatic Problems

We present a novel a posteriori error estimator for Nédélec elements for magnetostatic problems that is constant-free, i.e. it provides an upper bound on the error that does not involve a generic constant. The estimator is based on equilibration of the magnetic field and only involves small local problems that can be solved in parallel. Such an error estimator is already available for the lowest-degree Nédélec element (Braess and Schöberl in Math Comput 77(262):651-672, 2008) and requires solving local problems on vertex patches. The novelty of our estimator is that it can be applied to Nédélec elements of arbitrary degree. Furthermore, our estimator does not require solving problems on vertex patches, but instead requires solving problems on only single elements, single faces, and very small sets of nodes. We prove reliability and efficiency of the estimator and present several numerical examples that confirm this.


Introduction
We consider an a posteriori error estimator for finite element methods for solving equations of the form ∇ × μ −1 (∇ × u) = 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 c 0 estimator ≤ error ≤ c 1 estimator, 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 c 0 estimator ≤ error ≤ estimator, 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 βu, with β > 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.
in Ω, ∇ · μH = 0 i n Ω, where ∇ is the vector of differential operators (∂ 1 , ∂ 2 , ∂ 3 ), × and · denote the outer-and inner product, respectively, (therefore, ∇× and ∇· are the curl-and divergence operator, respectively),n denotes the outward pointing unit normal vector, μ : Ω → R + , with μ 0 ≤ μ ≤ μ 1 for some positive constants μ 0 and μ 1 , is a scalar magnetic permeability, and j : Ω → 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 H = μ −1 ∇ × u, where u : Ω → R 3 is a vector potential, and by solving the following problem for u: The second condition is only added to ensure uniqueness of u and is known as Coulomb's gauge. Now, for any domain D ∈ R n , let L 2 (D) m denote the standard Lebesque space of squareintegrable vector-valued functions u : D → R m equipped with norm u 2 D := D u · u dx and inner product (u, w) D = D u · w dx, and define the following Sobolev spaces: The weak formulation of problem (1) is finding u ∈ H 0 (curl; Ω) ∩ H (div 0 ; Ω) such that 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 Finally, let T h denote a tessellation of Ω into tetrahedra with a diameter smaller than or equal to h, let P −1 k (T h ), R −1 k (T h ), and D −1 k (T h ) denote the discontinuous spaces given by We define the finite element approximation for the magnetic vector potential as the vector field u h ∈ R k,0 (T h ) that solves The approximation of the magnetic field is then given by 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

An Equilibration-Based a Posteriori Error Estimator
We follow [6] and present an a posteriori error estimator that is based on the following result. then Proof The result follows from the orthogonality of μ 1/2 (H − H) and μ 1/2 (H − H h ): and Pythagoras's theorem (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 H h to obtain a solution to this equation is called equilibration of H h .

Corollary 3.3
Let u be the solution to (2), let u h be the solution of (3), set H := μ −1 ∇ ×u and H h := μ −1 ∇ ×u h , and let j h := ∇ ×H h be the discrete current distribution. IfH Δ ∈ L 2 (Ω) 3 satisfies the (residual) equilibrium condition which is an identity of distributions, then Proof Since (7) is an identity of distributions, we can equivalently write where ·, · denotes the application of a distribution to a function in C ∞ 0 (Ω) 3 . Now, set From this, it follows that ∇×H = j ∈ L 2 (Ω) 3 , soH is in H (curl; Ω) and satisfies equilibrium condition (4). Inequality (8) then follows from Theorem 3.1.
From Corollary 3.3, it follows that a constant-free upper bound on the error can be obtained from any fieldH Δ 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, j − j h is decomposed into a sum of local divergence-free current distributions j Δ i that have support on only a single vertex patch. The error estimator is then obtained by solving local problems of the form ∇ ×H Δ i = j Δ i for each vertex patch and by then taking the sum of all local fieldsH Δ i . It is, however, not straightforward to decompose j − j h into local divergence-free current distributions. An explicit expression for j Δ 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.

Assumptions
In order to compute the error estimator, we use polynomial function spaces of degree k ≥ k, where k denotes the degree of the finite element approximation u h , and assume that: A1. The magnetic permeability μ is piecewise constant. In particular, the domain Ω can be partitioned into a finite set of polyhedral subdomains {Ω i } such that μ is constant on each subdomain Ω i . Furthermore, the mesh is assumed to be aligned with this partition so that μ is constant within each element. A2. The current density j is in Although assumption A2 does not hold in general, we can always replace j by a suitable projection π h j by taking, for example, π h as the standard Raviart-Thomas interpolation operator corresponding to the D k (T h ) space [16]. The error is in that case bounded by where H := μ −1 ∇ × u and where u is the solution to (2) with j replaced by π h j. If j is sufficiently smooth, i.e. j can be extended to a function j * ∈ H (div; This means that, if k ≥ k, then μ 1/2 (H − H ) Ω converges with a higher rate than μ 1/2 (H − H h ) Ω and so we may assume that the term μ 1/2 (H − H ) Ω is negligible.

Derivation of the Error Estimator
Before we derive the error estimator, we first write j h = ∇ × H h in terms of element and face distributions.
where ·, · denotes the application of a distribution to a C ∞ 0 (Ω) 3 function, F I h denotes the set of all internal faces,n T denotes the normal unit vector to ∂ T pointing outward of T , and [[H]] t | f := (n + × H + +n − × H − )| f denotes the tangential jump operator, withn ± :=n T ± , H ± := H| T ± , and T + and T − the two adjacent elements of f . Since In other words, 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 j Δ := j − j h and can write where We look for a solution of (7) of the formH where ∇ h denotes the element-wise gradient operator. The termĤ Δ will take care of the element distributions of j Δ and the term ∇ h φ will take care of the remaining face distributions.
In the following, we firstly describe how to computeĤ Δ in Sect. 3.2.1 and characterize the remainder j Δ − ∇ ×Ĥ Δ in Sect. 3.2.2. We then describe how to compute the jumps of φ on internal faces in Sect. 3.2.3 and explain how to reconstruct φ from its jumps in Sect. 3.2.4. 1 We computeĤ Δ by solving the local problems

Computation ofĤ
for each element T ∈ T h . This problem is well-defined and has a unique solution due to the discrete exact sequence property . This last property follows from the fact that j| T ∈ D k due to assumption A2 and 3 due to assumption A1.

Representation of the Remainder j
. This means thatĵ Δ can be represented by only face distributions, and since

Computation of the Jumps of on Internal Faces
It now remains to find a φ ∈ P −1 To do this, we define, for each internal face f two orthogonal unit tangent vectorst 1 andt 2 such thatt 1 ×t 2 =n f := n + | f , differential operators ∂ t i :=t i · ∇, and the gradient operator restricted to the face: We therefore introduce an auxiliary variable λ f ∈ P k ( f ) and solve for each f ∈ F I h , 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 φ ∈ P −1 h . Now, we will prove that problem (13) uniquely defines λ f . We start by showing that (13) corresponds to a 2D curl problem on a face. To see this, note thatn f × ∇ f = t 2 ∂ t,1 −t 1 ∂ t,2 . If we take the inner product of (13a) witht 1 andt 2 , we obtain (13) is well-posed, we use the discrete exact sequence in 2D: To prove this, we use that, for every ψ ∈ C ∞ 0 (Ω) 3 , Then, for every ψ ∈ C ∞ 0 (Ω), we can write where E I h denotes the set of all internal edges,n ∂ f denotes the normal unit vector of ∂ f that lies in the same plane as f and points outward of f , andn e, f :=n ∂ f | e . This implies that so ∇ f ·ĵ Δ f = 0 for each internal face and, therefore, problem (13) is well-defined and has a unique solution.

Reconstruction of from its Jumps on the Internal Faces
After computing λ f for all internal faces, it remains to compute φ such that To do this, we use standard Lagrangian basis functions. For each element T , let Q T denote the set of nodes on element T for P k (T ). The barycentric coordinates of these nodes are given by Also, let Q h be the union of all element nodes. We then define the degrees of freedom for φ, denoted by We can decouple this global problem into very local problems. In particular, for each node x ∈ Q h , we can compute the small set of degrees of freedom {φ T ,x } T :T x by solving In this local problem, each degree of freedom corresponds to an element adjacent to x and for any two degrees of freedom corresponding to two adjacent elements, the difference should be equal to λ f (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 x in the interior of an element, there is only one overlapping element T and the above results in φ T ,x = 0. The same applies to a node in the interior of a boundary face. For a node x in the interior of an internal face f , there are only two adjacent elements T + and T − and the above results in φ T + ,x = 1 2 λ f (x) and φ T − ,x = − 1 2 λ f (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 x, the corresponding differences should add up to zero. A cycle means a sequence of elements T 1 → T 2 → · · · → T n+1 with n ≥ 3, such that T 1 , T 2 , . . . , 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 x are the cycles around each internal edge connected to 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 ∈ E I h be an internal edge, lett e be a tangent unit vector of e, and let T 1 → T 2 → · · · → T n+1 be a cycle rotating counter-clockwise when looking towardst e . This means that the normal unit vector of the face T i ∩ T i+1 pointing out of T i+1 is given byn T i+1 | f =t e ×n e, f =:n f ,e (recall thatn e, f =n ∂ f | e ). Also, let ψ ∈ P −1 k (T h ). The sum of the differences of ψ for this cycle can be written as Now, let f = T i ∩ T i+1 and note that Therefore, we can write The sum of the differences of the cycle around e can therefore be rewritten as and from this, we obtain the conditions Problem (16) is therefore well-posed provided that (17) is satisfied.
To prove (17), we first prove that, for each e ∈ E I h , r e is constant and then prove that, for each e ∈ E I h , (r e , 1) e = 0. Let e be an edge and let f be an adjacent face, and consider (14) witht 1 =t e and t 2 = (n f ·n f ,e )n e, f . An illustration oft e ,n e, f , andn f ,e is given in Fig. 1. The condition n f =t 1 ×t 2 is still satisfied, since eithern f =n f ,e orn f = −n f ,e and sô t 1 ×t 2 =t e ×n e, f (n f ·n f ,e ) =n f ,e (n f ·n f ,e ) =n f .

It then follows from (14b) that
where ∂ t e :=t e · ∇. Multiplying the above by −(n f ·n f ,e ) and summing over all faces adjacent to e results in where the last line follows from (15a). We thus have 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 θ e to be the lowest-order Nédélec basis function corresponding to an internal edge e and scaled such that (t e · θ e )| e = 1 and derive where the first two terms in the fifth line follow from the definition of u and u h and the last term in the fifth line follows from (11b), assumption A1, and the fact that ∇ × θ e is piecewise constant. We can then derive where the second line follows from the fact that the tangent components of θ e are zero on all faces that are not adjacent to e, and the last line follows from (13b) and the fact that ∇ f × θ e is constant on f . We continue to obtain where the third line follows from the fact that eithern f =n f ,ẽ orn f = −n f ,ẽ , the fifth line follows from the property (a × b) × c = −a(b · c) + b(a · c) and the fact thattẽ andnẽ , f are orthogonal, and the last line follows from the fact that (tẽ · θ e )|ẽ = 0 for all edgesẽ = e and from (t e · θ e )| e = 1. We then continue to obtain Therefore, (r e , 1) e = 0 and hence r e = 0 for all internal edges and so problem (16) is well-posed.

Summary of Computing the Error Estimator
We fix k ≥ k and compute our error estimator in four steps.
Step 1. We computeĤ Δ ∈ R −1 k (T h ) from the datum j and the numerical solution H h by solving for each T ∈ T h .
Step 2. For each internal face f ∈ F I h , let T + and T − denote the two adjacent elements, let n ± denote the normal unit vector pointing outward of T ± , let H ± := H| T ± denote the vector field restricted to T ± , let [[H]] t | f := (n + × H + +n − × H − )| f denote the tangential jump operator, and let ∇ f denote the gradient operator restricted to face f . We setn f :=n + | f and compute λ f ∈ P k ( f ) by solving for each internal face f ∈ F I h .
Step 3. We compute φ ∈ P −1 k (T h ) by solving, for each x ∈ Q h , the small set of degrees of freedom {φ T ,x } T :T x such that where Q h denotes the set of standard Lagrangian nodes corresponding to the finite element space P k (T h ) and φ T ,x denotes the value of φ| T at node x.
Step 4. We compute the fieldH where ∇ h denotes the element-wise gradient operator, and compute the error estimator

Remark 3.4
For the case k = k = 1, this algorithm requires solving local problems that involve 6 unknowns per element (Step 1), 3 unknowns per face (Step 2), and (#T : T ν) ≈ 24 unknowns per vertex ν (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].

Reliability and Efficiency
The results of the previous sections immediately give the following theorem.
We also prove local efficiency of the estimator and state it as the following theorem.

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 , as we used discrete inverse inequalities and the efficiency result stated in [4]. Whether it is also independent of k is an open issue.

Numerical Experiments
In this section, we present several numerical results for the unit cube and the L-brick domain with constant magnetic permeability μ = 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 u h , i.e. k = k. In the numerical experiments, we do not project the right hand side j onto D k (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 φ 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 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 H h := μ −1 ∇ × 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 (T h ) is based on the hierarchical basis functions from [22].
We evaluate the reliability and efficiency of the a posteriori error estimator η 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 η T of the estimator, we mark elements for refinement based on the bulk marking strategy [12] with bulk parameter θ = 0.5, and finally we refine the marked elements based on a bisection strategy [3].

Unit Cube Example
For the first example we chose the unit cube Ω = (0, 1) 3 with homogeneous boundary conditions and right hand side j according to the polynomial solution In Fig. 2, we present the errors H − H h Ω and efficiency indices η h / H − H h Ω for k = 1, 2, 3, on a sequence of uniformly refined meshes. We observe that the error converges with optimal rates O(h k ) = O(N −k/3 h ), N h = dim(R k (T h )), and that the error estimator is reliable and efficient with efficiency constants between 1 and 2.    Note that, for k = 3, j belongs to D k (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 j does not belong to D k (T h ). Thus the error introduced in the computation of the error estimator by not satisfying A2 is negligible.

L-Brick Example
As second example, we solve the homogeneous Maxwell problem on the (nonconvex) domain As solution, we choose the singular function where (r , ϕ) are the two dimensional polar coordinates in the x-y-plane, and we choose the right hand side j accordingly. Due to the edge singularity, uniform mesh refinement leads to suboptimal convergence rates of O(N −2/9 h ), 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 . 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 η h for k-refinement to those of the residual a posteriori error estimator [4] 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 μ h is not robust in the polynomial degree k, whereas our new estimator appears to be more robust in k.

Example with Discontinuous Permeability
In this example, we choose a discontinuous permeability μ(x, y, z) = μ 1 if y < 1/2 and z < 1/2, μ 2 otherwise, on the unit cube Ω = (0, 1) 3 and the right hand side j = (1, 0, 0) t . We choose μ 1 = 1 and vary μ 2 = 10 for = 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 (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 · 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 .

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.