We consider an a posteriori error estimator for finite element methods for solving equations of the form . 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 . It is a residual-type estimator and provides bounds of the form
up to some higher-order data oscillation terms, where are positive constants that do not depend on the mesh resolution. Similar bounds can be obtained by hierarchical error estimators; see, e.g., , under the assumption of a saturation condition, and by Zienkiewicz–Zhu-type error estimators; see, e.g., . A drawback of these estimators is that the constants and 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 , these estimators have become a major research topic; for a recent overview, see, for example,  and the references therein. An equilibration-based error estimator was introduced for magnetostatics in  and provides bounds of the form
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  and, for an eddy current problem, in [11, 10]. Constant-free upper bounds are also obtained by the functional estimate in , when selecting a proper function in their estimator, and by the recovery-type error estimator in , in case the equations contain an additional term , with .
A drawback of the estimators in [19, 11, 10, 16] is that they require solving a global problem. The estimator in , 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 Section 2. We then derive our error estimator step by step in Section 3, with a summary given in Section 3.3, and prove its reliability and efficiency in Section 3.4. Numerical examples confirming the reliability and efficiency of our estimator are presented in Section 4, and an overall summary is given in Section 5.
2. A finite element method for magnetostatic problems
Let be an open, bounded, simply connected, polyhedral domain with a connected Lipschitz boundary . In case of a linear, isotropic medium and a perfectly conducting boundary, the static magnetic field satisfies the equations
is the vector of differential operators, and denote the outer- and inner product, respectively, (therefore, and are the curl- and divergence operator, respectively), denotes the outward pointing unit normal vector, , with for some positive constants and , is a scalar magnetic permeability, and is a given divergence-free current density. The first equality is known as Ampères law and the second as Gauss’s law for magnetism.
These equations can be solved by writing , where is a vector potential, and by solving the following problem for :
The second condition is only added to ensure uniqueness of and is known as Coulomb’s gauge.
Now, for any domain , let denote the standard Lebesque space of square-integrable vector-valued functions equipped with norm and inner product , and define the following Sobolev spaces:
The weak formulation of problem (1) is finding 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 be a tetrahedron and define to be the space of polynomials on of degree or less. Also, define the Nédélec space of the first kind and the Raviart-Thomas space by
Finally, let denote a tessellation of into tetrahedra with a diameter smaller than or equal to , let , , and denote the discontinuous spaces given by
We define the finite element approximation for the magnetic vector potential as the vector field that solves
The approximation of the magnetic field is then given by
which converges quasi-optimally as the mesh width tends to zero [15, Theorem 5.10].
In the next section, we show how we can obtain a reliable and efficient estimator for .
3. An equilibration-based a posteriori error estimator
We follow  and present an a posteriori error estimator that is based on the following result.
Theorem 3.1 ([7, Thm. 10]).
The result follows from the orthogonality of and :
and Pythagoras’s theorem
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 to obtain a solution to this equation is called equilibration of .
Since (7) is an identity of distributions, we can equivalently write
where denotes the application of a distribution to a function in . Now, set . Using the definition , we obtain
An error estimator of this type was first introduced in , where it is referred to as an equilibrated residual error estimator. There, is decomposed into a sum of local divergence-free current distributions that have support on only a single vertex patch. The error estimator is then obtained by solving local problems of the form for each vertex patch and by then taking the sum of all local fields . It is, however, not straightforward to decompose into local divergence-free current distributions. An explicit expression for is given in  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 Sections 3.1 and 3.2 below, a brief summary is given in Section 3.3, and reliability and efficiency are proven in Section 3.4.
In order to compute the error estimator, we use polynomial function spaces of degree , where denotes the degree of the finite element approximation , and assume that:
The magnetic permeability is piecewise constant. In particular, the domain can be partitioned into a finite set of polyhedral subdomains such that is constant on each subdomain . Furthermore, the mesh is assumed to be aligned with this partition so that is constant within each element.
The current density is in .
Although assumption A2 does not hold in general, we can always replace by a suitable projection by taking, for example,
as the standard quasi-interpolation operator of thespace. The error is in that case bounded by
where is some positive constant that does not depend on the mesh width . If is sufficiently smooth, then the data oscillation term is of order , so if , then this term converges with a higher rate than and we may assume that it is negligible.
3.2. Derivation of the error estimator
Before we derive the error estimator, we first write in terms of element and face distributions. For every , we can write
where denotes the application of a distribution to a function, denotes the set of all internal faces, denotes the normal unit vector to pointing outward of , and denotes the tangential jump operator, with , , and and the two adjacent elements of .
Since and is piecewise constant, we have that . Therefore, if and if , and , where is a normal unit vector of and is given by
In other words, can be represented by functions on the elements and face distributions on the internal faces.
We define and can write
We look for a solution of (7) of the form , with and and where denotes the element-wise gradient operator. The term will take care of the element distributions of and the term will take care of the remaining face distributions.
In the following, we firstly describe how to compute in Section 3.2.1 and characterize the remainder in Section 3.2.2. We then describe how to compute the jumps of on internal faces in Section 3.2.3 and explain how to reconstruct from its jumps in Section 3.2.4.
3.2.1. Computation of
We compute by solving the local problems
for each element . This problem is well-defined and has a unique solution due to the discrete exact sequence property
and since and . This last property follows from the fact that due to assumption A2 and due to assumption A1.
3.2.2. Representation of the remainder
For every , we can write
for all . This means that can be represented by only face distributions, and since and , we have that and therefore .
3.2.3. Computation of the jumps of on internal faces
It now remains to find a such that
For every , we can write
Therefore, we need to find a such that
To do this, we define, for each internal face , the scalar jump with , two orthogonal unit tangent vectors and such that , differential operators , and the gradient operator restricted to the face: . We can then write
for all . We therefore introduce an auxiliary variable and solve
for each , 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 such that for all . Now, we will prove that problem (13) uniquely defines . We start by showing that (13) corresponds to a 2D curl problem on a face. To see this, note that . If we take the inner product of (13a) with and , we obtain
where , which is equivalent to a 2D curl problem on . To show that (13) is well-posed, we use the discrete exact sequence in 2D:
where . Since , it suffices to show that . To prove this, we use that, for every ,
Then, for every , we can write