# Variational methods for solving numerically magnetostatic systems

Patrick Ciarlet, E. Jamelot

► To cite this version:

Patrick Ciarlet, E. Jamelot. Variational methods for solving numerically magnetostatic systems. Advances in Computational Mathematics, inPress. hal-04152740v2

HAL Id: hal-04152740

https://inria.hal.science/hal-04152740v2

Submitted on 8 Nov 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Variational methods for solving numerically
magnetostatic systems

Patrick Ciarlet Jr.<sup>1*†</sup> and Erell Jamelot<sup>2†</sup>

<sup>1</sup>*POEMS, CNRS, INRIA, ENSTA Paris, Institut Polytechnique de Paris, 91120 Palaiseau, France.*

<sup>2</sup>Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191 Gif-sur-Yvette, France.

*Corresponding author(s). E-mail(s):

patrick.ciarlet@ensta-paris.fr;

Contributing authors: erell.jamelot@cea.fr;

†These authors contributed equally to this work.

**Abstract**

In this paper, we study some techniques for solving numerically magnetostatic systems. We consider fairly general assumptions on the magnetic permeability tensor. It is elliptic, but can be nonhermitian. In particular, we revisit existing classical variational methods and propose new numerical methods. The numerical approximation is either based on the classical edge finite elements, or on continuous Lagrange finite elements. For the first type of discretization, we rely on the design of a new, mixed variational formulation that is obtained with the help of **T**-coercivity. The numerical method can be related to a perturbed approach for solving mixed problems in electromagnetism. For the second type of discretization, we rely on an augmented variational formulation obtained with the help of the Weighted Regularization Method.

**Keywords:** Magnetostatic systems, variational formulations, *T*-coercivity, edge finite elements, Lagrange finite elements
## 2 Solving magnetostatic systems

# 1 Introduction

In this paper, we revisit existing classical variational methods for solving magnetostatic models, and we also propose a new method. The magnetic permeability tensor $\mu$ is assumed to be tensor-valued, and elliptic. We study the methods from the mathematical and discretization points of view. Unless otherwise specified, the model is set in a region of $\mathbb{R}^3$, that is in a 3D domain, and we give all definitions assuming this geometrical setting. At some point, we consider models set in a region of $\mathbb{R}^2$, and leave to the reader the straightforward adaptation of the mathematical tools.

The outline is as follows. In the next section, we define the mathematical setting. In section 3, we introduce the model, whose solution is the magnetic induction, and the mathematical setting. Then, in section 4, we recall how the magnetic induction can be split into two parts. First, a div- and curl-free part, that belongs to a finite dimensional vector space. As claimed for instance in [1], "it is an easy finite dimensional problem", so we do not dwell on this issue. Second, a part that can be recast as the curl of a vector field, the so-called vector potential. For the sake of comparison, we propose two strategies to study and discretize the problem of finding the vector potential and the corresponding magnetic induction. In section 5, we start with a formulation in $H_0(\text{curl}; \Omega)$ that is obtained through the use of the so-called T-coercivity [2, 3], and that can be discretized with the help of the classical edge finite elements [4, 5]. In section 6, we consider formulations that can be discretized with continuous vector finite elements. For that, one needs to take into account a measure of the divergence of the fields: to that aim, we choose the Weighted Regularization Method (WRM) of Costabel-Dauge [6]. Finally in section 7, we give some concluding remarks on the relative merits of both approaches.

## 2 Mathematical setting

We consider function spaces of complex-valued functions. Vector-valued function spaces are written in boldface character. The index $zmv$ indicates zero-mean-value fields. Duality brackets between a Banach space $X$ and its topological dual $X'$ are denoted by $\langle \cdot, \cdot \rangle_X$. Given a non-empty open set $O$ of $\mathbb{R}^3$ with a Lipschitz boundary $\partial O$, $n$ denotes the unit outward normal vector field to $\partial O$. It is assumed that the reader is familiar with function spaces related to Maxwell's equations, such as $H(\text{curl}; O)$, $H_0(\text{curl}; O)$, $H(\text{div}; O)$, $H_0(\text{div}; O)$ etc. We consider a priori that $H(\text{curl}; O)$ is endowed with the “natural” norm $v \mapsto (\|v\|_{L^2(O)}^2 + \|\text{curl } v\|_{L^2(O)}^2)^{1/2}$, etc. We refer to the monographs [5, 7, 8] for details.

The model is set in a domain $\Omega$ in $\mathbb{R}^3$, i.e. an open, connected and bounded subset of $\mathbb{R}^3$ with a Lipschitz-continuous boundary $\partial\Omega$. The domain $\Omega$ can be simply connected or not [9]. This property is mathematically characterized by the value of the first Betti number, $\beta_1(\Omega)$ : if $\beta_1(\Omega) = 0$, the domain $\Omega$ is
simply connected and, if $\beta_1(\Omega) > 0$, it is non-simply connected. This means that we assume that one of the two conditions below holds:

* $\beta_1(\Omega) = 0$ 'for all vector fields $v \in C^0(\Omega)$ such that $\text{curl } v = 0$ in $\Omega$, there exists $p \in C^1(\Omega)$ such that $v = \nabla p$ in $\Omega'$;

* $\beta_1(\Omega) > 0$ 'there exist $I = \beta_1(\Omega)$ non-intersecting, piecewise plane manifolds $(\Sigma_j)_{j=1,\dots,I}$, called cuts, with boundaries $\partial\Sigma_i \subset \partial\Omega$, such that, if one introduces the connected set $\dot{\Omega} = \Omega \setminus \bigcup_{i=1}^I \Sigma_i$, for all vector fields $v \in C^0(\Omega)$ such that $\text{curl } v = 0$ in $\Omega$, there exists $\dot{p} \in C^1(\dot{\Omega})$ such that $v = \nabla \dot{p}$ in $\dot{\Omega}$.'

Given $v \in L^2(\dot{\Omega})$ (resp. $v \in L^2(\Omega)$), we denote by $\tilde{v}$ (resp. $\tilde{v}$) its continuation to $L^2(\Omega)$ (resp. $L^2(\Omega)$).

The boundary $\partial\Omega$ is split into $K + 1$ maximal connected components $(\Gamma_k)_{k=0,K}$, where $K = \beta_2(\Omega)$ is the second Betti number. We denote by $\Gamma_0$ the boundary of the unbounded component of $\mathbb{R}^3 \setminus \bar{\Omega}$.

Given $\xi \in (L^\infty(\Omega))^{3 \times 3}$, we use the notation $\xi_+ = \|\xi\|_{(L^\infty(\Omega))^{3 \times 3}}$, and define

$$ H_0(\text{div } \xi; \Omega) = \{v \in L^2(\Omega) \text{ such that } \xi v \in H_0(\text{div}, \Omega)\}. $$

Furthermore, we say that $\xi \in (L^2(\Omega))^{3 \times 3}$ is piecewise smooth if there exists a partition<fnref n="1" /> $(\Omega_p^\xi)_{p=1,P}$ of $\Omega$ so that $\xi|_{\Omega_p^\xi} \in (W^{1,\infty}(\Omega_p^\xi))^{3 \times 3}$ for $p = 1, P$. In this case, the partition is called compatible (with respect to $\xi$). Finally, for $s > 0$, we introduce

$$ PH^s(\Omega) = \{f \in L^2(\Omega) \text{ such that } f|_{\Omega_p^\xi} \in H^s(\Omega_p^\xi), \forall p = 1, P\}. $$

The symbol $C$ is used to denote a generic positive constant which is independent of the meshsize, the mesh and the fields of interest; $C$ may depend on the geometry, or on the coefficients defining the model. We use the notation $A \lesssim B$ for the inequality $A \le CB$, where $A$ and $B$ are two scalar fields, and $C$ is a generic constant. If the inequality depends on some parameter $s$, we write $A \lesssim_s B$.

## 3 Model and assumptions

We consider a material, characterized by the magnetic permeability tensor $\mu$. Neglecting the time-variation of the electric displacement, one finds that the magnetic induction $B \in L^2(\Omega)$ is governed by the so-called magnetostatic system of equations. We refer to [10, Chapter 6] or [11, Chapter 3] for its

<footnote n="1">¹We recall that a partition of $\Omega$ is $(\Omega_p)_{p=1,P}$ such that</footnote>

<footnote>$\Omega_p$ is a domain, for $p = 1, P$; $\Omega_p \cap \Omega_q = \emptyset$ for $p \ne q$; $\bar{\Omega} = \cup_{p=1,P} \bar{\Omega}_p$.</footnote>
4 Solving magnetostatic systems

definition. See also [12–17] and References therein for related models.

$$
\left\{
\begin{array}{l}
\text{Find } \mathbf{B} \in \mathbf{L}^2(\Omega) \text{ such that} \\
\text{curl } (\mu^{-1}\mathbf{B}) = \mathbf{J} \text{ in } \Omega \\
\text{div } \mathbf{B} = 0 \text{ in } \Omega \\
\mathbf{B} \cdot \mathbf{n} = 0 \text{ on } \partial\Omega \\
\langle \mathbf{B} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Sigma_i)} = c_i, \forall i = 1, I
\end{array}
\right.
\tag{1}
$$

Above, the pair of data $(\mathbf{J}, (c_i)_{i=1,I})$ is: $\mathbf{J} \in \mathbf{L}^2(\Omega)$ the electric current density in $\Omega$; the numbers $(c_i)_{i=1,I}$ are related to current intensities flowing through some ad hoc exterior surfaces, see [11, Section 3.6] or section 4.2 below. It is also possible to use a formulation with the magnetic field $\mathbf{H} = \mu^{-1}\mathbf{B}$ as the main unknown. Because of the first equation, written $\mathbf{J} = \text{curl } \mathbf{H}$ with $\mathbf{H} \in \mathbf{H}(\text{curl}; \Omega)$, we note that $\mathbf{J} \in \mathbf{H}(\text{div}; \Omega)$ with $\text{div} \mathbf{J} = 0$ in $\Omega$. Also, one may check that $\langle \mathbf{J} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_k)} = 0$, for all $k=1, K$, as in Remark 3.4.2 in [8].

In what follows, we therefore consider any pair of data $(\mathbf{J}, (c_i)_{i=1,I})$ with

$$
\mathbf{J} \in \mathbf{H}(\mathrm{div}; \Omega) \text{ s.t. div } \mathbf{J} = 0 \text{ in } \Omega, \langle \mathbf{J} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_k)} = 0, k = 1, K. \quad (2)
$$

Mathematically, we assume that the permeability $\mu$ is a complex-valued, measurable, bounded, elliptic (see next definition) and smooth or piecewise smooth tensor field. We observe that, in this case, the magnetic induction and the data are a priori complex-valued. We recall definitions and properties that can be found e.g. in [18, 19]. From now on, piecewise smoothness is understood with respect to a given compatible partition $(\Omega_p^\mu)_{p=1,P}$.

**Definition 1** The tensor-valued field $\xi \in (L^\infty(\Omega))^{3 \times 3}$ is elliptic if

$$
\exists (\theta_{\xi}, \xi_{-}) \in \mathbb{R} \times \mathbb{R}_{>0}, \text{ a.e. in } \Omega, \forall \mathbf{z} \in \mathbb{C}^{3}, \xi_{-} | | \mathbf{z} |^{2} \leq \Re[e^{i \theta_{\xi}} \xi \mathbf{z} \cdot \bar{\mathbf{z}}]. \quad (3)
$$

In (3), $\theta_{\xi}$ is called an ellipticity direction. We denote by $\Theta_{\xi}$ the set of admissible ellipticity directions

$$
\Theta_{\xi} = \{ \theta_{\xi} \in \mathbb{R}, (\theta_{\xi}, \xi_{-}) \text{ fulfills (3) for some } \xi_{-} \in \mathbb{R}_{>0} \}.
$$

**Proposition 1** Let $\zeta \in (L^\infty(\Omega))^{3 \times 3}$. If $\zeta$ is elliptic, one has $\zeta^{-1} \in (L^\infty(\Omega))^{3 \times 3}$ with

$$
(\zeta^{-1})_{+} \le \inf_{(\theta_{\xi}, \xi_{-}) \text{ s.t. (3) with } \xi = \zeta \text{ holds}} ((\zeta^{-1})^{-1}).
$$

Moreover, $\Theta_{\zeta^{-1}} = \{-\theta_{\zeta}, \theta_{\zeta} \in \Theta_{\zeta}\}$ and, given any $\theta_{\zeta^{-1}} \in \Theta_{\zeta^{-1}}$, one can choose $(\zeta^{-1})_{-} = \zeta_{-}(\zeta_{+})^{-2}$ in (3) with $\xi = \zeta^{-1}$, where $(-\theta_{\zeta^{-1}}, \zeta_{-})$ is such that (3) with $\xi = \zeta$ holds.

For conciseness, we assume that $0 \in \Theta_{\mu}$, so that $0 \in \Theta_{\mu-1}$ according to the above (see footnote² page 10 for the case where $0 \notin \Theta_{\mu-1}$). So, for the possibly
nonhermitian elliptic tensor field $\mu$, denoting by $\mu_-$ a value in (3) with $\xi = \mu$ for $\theta_\mu = 0$, one knows that $(\mu^{-1})_+ \le (\mu_-)^{-1}$, and that $(\mu^{-1})_- = \mu_-(\mu_+)^{-2}$ is an admissible value in (3) with $\xi = \mu^{-1}$ for $\theta_{\mu^{-1}} = 0$. We use those values throughout the manuscript.

For the ease of exposition and as far as discretization is concerned, we assume that $\Omega$ is a Lipschitz polyhedron. As a matter of fact, if the boundary is of class $C^2$, using the theory developed e.g. in [20, 21], one ends up with convergence results that are identical to those we obtain. The polyhedron $\bar{\Omega}$ is triangulated by a shape regular family of meshes $(\mathcal{T}_h)_h$, made up of (closed) simplices $K$. Further, we assume that the meshes are conforming with the compatible partition $(\Omega_p^\mu)_{p=1,P}$ and/or with respect to the cuts. In other words, for all $h$, for all $K \in \mathcal{T}_h$, there exists $p \in \{1, \dots, P\}$ such that $K \subset \bar{\Omega}_p$ and $\text{int}(K) \cap (\cup_{i=1,I} \Sigma_i) = \emptyset$. This is a realistic assumption that allows to simplify the numerical analysis presented here.

# 4 A reformulation of the problem

By linearity, we see that solving the magnetostatic problem (1) can be done in two independent steps, by splitting the data: a first solve with the pair $(\mathbf{J}, 0)$ with a solution called $\mathbf{B}_J$, and a second one with the pair $(0, (c_i)_{i=1,I})$, with a solution called $\mathbf{B}_c$. The total solution is then $\mathbf{B} = \mathbf{B}_J + \mathbf{B}_c$. This is this approach that we follow below.

## 4.1 Introducing the vector potential

The pair of data is equal to $(\mathbf{J}, 0)$, with $\mathbf{J} \in L^2(\Omega)$. In this case, we note that the magnetic induction fulfills the conditions

$$ \mathbf{B}_J \in H_0(\text{div}; \Omega), \text{ div } \mathbf{B}_J = 0 \text{ in } \Omega \text{ and } \langle \mathbf{B}_J \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Sigma_i)} = 0, \forall i = 1, I. $$

According to e.g. Theorem 3.5.1. in [8], there exists a vector potential $\mathbf{A} \in H_0(\text{curl}; \Omega)$ such that $\mathbf{B}_J = \text{curl} \mathbf{A}$ in $\Omega$. Moreover, one may characterize the vector potential by choosing a divergence-free, flux-free potential $\mathbf{A}$, i.e. $\text{div} \mathbf{A} = 0$ in $\Omega$ and $\langle \mathbf{A} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_k)} = 0$, for all $k = 0, K$. Obviously, the converse assertion is true: such a field $\mathbf{A}$ is such that $\mathbf{B}_J = \text{curl} \mathbf{A}$ solves the above. So we conclude that, with the pair of data $(\mathbf{J}, 0)$, the magnetostatic system is equivalently recast as

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{A} \in H_0(\text{curl}; \Omega) \text{ such that} \\ \text{curl } (\mu^{-1}\text{curl } \mathbf{A}) = \mathbf{J} \text{ in } \Omega \\ \text{div } \mathbf{A} = 0 \text{ in } \Omega \\ \langle \mathbf{A} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_k)} = 0, \forall k = 1, K \end{array} \right. \qquad (4) $$
6 Solving magnetostatic systems

and then one defines $B_J = \text{curl } A$. We explain how to solve the problem with unknown $A$, using variational formulations set either in $H_0(\text{curl}; \Omega)$, or in a proper subspace of $H_0(\text{curl}; \Omega)$. More precisely, in section 5, we start with a formulation in $H_0(\text{curl}; \Omega)$, discretized with the help of the classical edge finite elements [4, 5]. Then, in section 6, we consider formulations in

$$ \{v \in H_0(\text{curl}; \Omega) \text{ s.t. } \operatorname{div} v \in Y(\Omega)\}, \qquad (5) $$

where $Y(\Omega)$ is a function space such that $L^2(\Omega) \subset Y(\Omega) \subset H^{-1}(\Omega)$. We focus on the so-called Weighted Regularization Method (WRM) of Costabel-Dauge [6], that can be discretized with continuous vector finite elements. To our knowledge, this is the first method that was designed to approximate successfully singular fields in 3D domains by such finite elements. For other methods using the same mathematical framework, we refer to [22–26].

## 4.2 Finding curl- and div-free fields

The pair of data is now equal to $(0, (c_i)_{i=1,I})$, with $(c_i)_{i=1,I} \in \mathbb{C}^I$. In this second case, we note that $\mu^{-1}B_c \in H(\text{curl}; \Omega)$ fulfills in particular the condition $\text{curl}(\mu^{-1}B_c) = 0$ in $\Omega$. We let

$$ P(\dot{\Omega}) = \{\dot{q} \in H^1(\dot{\Omega}) \text{ such that } [\dot{q}]_{\Sigma_i} \in \mathbb{C}, \forall i = 1, I\}. $$

One uses the norm $q \mapsto \|\nabla \dot{q}\|_{L^2(\dot{\Omega})}$ in $P_{zmv}(\dot{\Omega})$. According to theorem 3.2.2 in [8], there exists one, and only one, $\dot{p} \in P_{zmv}(\dot{\Omega})$ such that $\mu^{-1}B_c = \nabla \dot{p}$ in $\Omega$. Hence, with the pair of data $(0, (c_i)_{i=1,I})$, the magnetostatic system is equivalently recast as

$$ \left\{ \begin{array}{l} \text{Find } \dot{p} \in P_{zmv}(\dot{\Omega}) \text{ such that} \\ \operatorname{div} \mu \nabla \dot{p} = 0 \text{ in } \Omega \\ \mu \nabla \dot{p} \cdot \mathbf{n} = 0 \text{ on } \partial \Omega \\ \langle \mu \nabla \dot{p} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Sigma_i)} = c_i, \forall i = 1, I \end{array} \right. , \qquad (6) $$

and then one defines $B_c = \mu \nabla \dot{p}$. This is a finite dimensional problem. One can solve it with unknown $\dot{p}$, using either a variational formulation in $H^1(\dot{\Omega})$ discretized with conforming Lagrange finite elements, or a mixed approach, discretized with Raviart-Thomas finite elements. The latter yields a conforming approximation $H_0(\operatorname{div}; \Omega)$. We observe that, once $\dot{p}$ is known, so are the values $([\dot{p}]_{\Sigma_i})_{i=1,P}$ of the jumps across the cuts, which are equal to the current intensities flowing through the exterior surfaces, cf. [11, Section 3.6].
## 5 Conforming variational formulations in $H_0(\text{curl}; \Omega)$

Given data $J$ that fulfills assumption (2), we now build a first variational formulation, that is equivalent to (4). Let us start by simple observations. We notice that, given $u \in H_0(\text{curl}; \Omega)$, $u$ is such that $\text{curl}(\mu^{-1}\text{curl}u) = J$ in $\Omega$ if, and only if,

$$ \forall v \in H_0(\text{curl}; \Omega), \quad (\mu^{-1}\text{curl}u, \text{curl}v)_{L^2(\Omega)} = (J, v)_{L^2(\Omega)}. \quad (7) $$

Then, let

$$ H_{\partial\Omega}^{1}(\Omega) = \{q \in H^{1}(\Omega) \text{ such that } q|_{\Gamma_{0}} = 0, q|_{\Gamma_{k}} \in \mathbb{C}, \forall k = 1, K\}. $$

One uses the norm $q \mapsto \|\nabla q\|_{L^2(\Omega)}$ in $H_{\partial\Omega}^1(\Omega)$.

One checks easily that, given $u \in L^2(\Omega)$, $u$ is such that $\text{div} u = 0$ in $\Omega$ and $(u \cdot n, 1)_{H^{1/2}(\Gamma_k)} = 0$, for $k = 1, K$ if, and only if,

$$ \forall q \in H_{\partial\Omega}^{1}(\Omega), \quad (u, \nabla q)_{L^2(\Omega)} = 0. \quad (8) $$

Hence, it is equivalent to solve (4), or the variational formulation

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{A} \in H_0(\text{curl}; \Omega) \text{ such that} \\ \forall \mathbf{v} \in H_0(\text{curl}; \Omega), \ (\mu^{-1}\text{curl}\mathbf{A}, \text{curl}\mathbf{v})_{L^2(\Omega)} = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}. \quad (9) \\ \forall q \in H_{\partial\Omega}^{1}(\Omega), \quad (\mathbf{A}, \nabla q)_{L^2(\Omega)} = 0 \end{array} \right. $$

### 5.1 A symmetrized formulation with artificial pressure

Classically, we build a second variational formulation, that is symmetrized via the introduction of an artificial pressure. As a matter of fact, given $\gamma > 0$, we aim at solving the symmetrized, complex-valued saddle-point formulation

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{A}, p) \in H_0(\text{curl}; \Omega) \times H_{\partial\Omega}^1(\Omega) \text{ such that} \\ \forall \mathbf{v} \in H_0(\text{curl}; \Omega), \ (\mu^{-1}\text{curl}\mathbf{A}, \text{curl}\mathbf{v})_{L^2(\Omega)} \\ \qquad \qquad \qquad \qquad + \gamma (\mathbf{v}, \nabla p)_{L^2(\Omega)} = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}. \quad (10) \\ \forall q \in H_{\partial\Omega}^1(\Omega), \quad \gamma (\mathbf{A}, \nabla q)_{L^2(\Omega)} = 0 \end{array} \right. $$

Clearly, if $\mathbf{A}$ solves (9), then $(\mathbf{A}, 0)$ solves (10). And, thanks to the assumption (2) on $J$, we have that any solution $(\mathbf{A}, p)$ to (10) is such that $p = 0$. As a matter of fact, we note that, since $p \in H_{\partial\Omega}^1(\Omega)$ it holds that $\mathbf{v} = \nabla p \in H_0(\text{curl}; \Omega)$. Inserting in (10), we find by integrating by parts that

$$ \gamma \|\nabla p\|_{L^2(\Omega)}^2 = (\mathbf{J}, \nabla p)_{L^2(\Omega)} = -(\text{div} \mathbf{J}, p)_{L^2(\Omega)} + \sum_{k=0,K} \langle \mathbf{A} \cdot \mathbf{n}, p \rangle_{H^{1/2}(\Gamma_k)} = 0. $$

So, one concludes that
8 Solving magnetostatic systems

**Theorem 2** Let $\gamma > 0$ be given, and let the data $J$ fulfill assumption (2). Then, the vector potential $A$ solves (4) if, and only if, $(A, 0)$ solves (10).

The next step is to prove that the variational formulation (10) is well-posed. This is again a classical result. We use below the theory of $T$-coercivity, that will help design a new approach to solve the problem, both at the theoretical level, and numerically.

## 5.2 Reminders about $T$-coercivity

We recall here the so-called $T$-coercivity theory, cf. [2, 3]. Let $V$ and $W$ be two Hilbert spaces, and consider the variational formulation

$$ \text{Find } u \in V \text{ such that } \quad \forall w \in W, a(u, w) = \ell(w), \qquad (11) $$

where $a$ is a continuous, sesquilinear form on $V \times W$, and $\ell \in V'$. Classically, the variational formulation (11) is well-posed if the form $a$ satisfies the stability (or inf-sup) condition and the solvability condition of the Banach–Nečas–Babuška (BNB) Theorem (see for instance [27, Thm. 25.9]). One may introduce an equivalent condition.

**Definition 2** Let $V$ and $W$ be two Hilbert spaces and $a$ be a continuous and sesquilinear form on $V \times W$. It is $T$-coercive if

$$ \exists T \in \mathcal{L}(V, W) \text{ bijective, } \exists \alpha > 0, \forall v \in V, |a(v, Tv)| \ge \alpha \|v\|_V^2. \qquad (12) $$

In other words, the form $a(\cdot, T\cdot)$ is coercive on $V \times V$. Since $T$ is bijective, one recovers well-posedness through the use of Lax-Milgram theorem for the variational formulation

$$ \text{Find } u \in V \text{ such that } \quad \forall v \in V, a(u, Tv) = \ell(Tv), \qquad (13) $$

which is equivalent to (11). Whereas the BNB theorem relies on an abstract inf-sup condition, finding a mapping $T$ allows one to realise explicitly both the stability (or inf-sup) condition and the solvability condition. Interestingly, in general one can derive discrete counterparts of the mapping $T$ to obtain a uniform discrete inf-sup condition, which guarantees convergence. Among others, this method has been used by the two co-authors and co-workers for classical models, such as neutron diffusion problems [28–30] and Stokes problems [31, 32].

## 5.3 Using $T$-coercivity

Let us apply the theory of $T$-coercivity to solve the variational formulation (10). We refer to the classroom notes [33], or to [31], for details. Let $V =$
$H_0(\text{curl}; \Omega) \times H_{\partial\Omega}^1(\Omega)$, endowed with the norm

$$ \|(v, q)\|_{V} = (\|v\|_{H(\text{curl}; \Omega)}^{2} + \|\nabla q\|_{L^{2}(\Omega)}^{2})^{1/2}. $$

On can rewrite (10) as

$$ \left\{ \begin{array}{l} \text{Find } (\mathbf{A}, p) \in V \text{ such that} \\ \forall (\mathbf{v}, q) \in V, a((\mathbf{A}, p), (\mathbf{v}, q)) = \ell((\mathbf{v}, q)), \end{array} \right. \qquad (14) $$

with

$$ \begin{aligned} a((\mathbf{u}, p), (\mathbf{v}, q)) &= (\mu^{-1}\mathbf{curl}\mathbf{u}, \mathbf{curl}\mathbf{v})_{L^2(\Omega)} + \gamma(\mathbf{v}, \nabla p)_{L^2(\Omega)} + \gamma(\mathbf{u}, \nabla q)_{L^2(\Omega)} \\ \ell((\mathbf{v}, q)) &= (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}. \end{aligned} $$

Above, $a$ is a continuous, sesquilinear form on $V \times V$, respectively $\ell$ an antilinear continuous form on $V$. We shall prove that the form $a(\cdot, \cdot)$ is T-coercive. In $H_0(\text{curl}; \Omega)$, let us introduce the orthogonal subspace to the range of the gradient operator from $H_{\partial\Omega}^1(\Omega)$:

$$ H_0(\text{curl}; \Omega) = \nabla[H_{\partial\Omega}^1(\Omega)]^{\perp} \oplus K_N^{-}(\Omega), \qquad (15) $$

where the orthogonality is understood with respect to the "natural" inner product $(\cdot, \cdot)_{H(\text{curl}; \Omega)}$. The decomposition (15) is usually called a Helmholtz decomposition. It is easily checked that

$$ \begin{aligned} K_N^{-}(\Omega) &= \{\mathbf{k} \in H_0(\text{curl}; \Omega) \text{ s.t. } \forall q \in H_{\partial\Omega}^1(\Omega), (\mathbf{k}, \nabla q)_{L^2(\Omega)} = 0\}, \\ &= \{\mathbf{k} \in H_0(\text{curl}; \Omega) \text{ s.t. } \text{div } \mathbf{k} = 0 \text{ in } \Omega, \langle \mathbf{k} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_k)} = 0, \forall \mathbf{k}\}. \end{aligned} $$

Above, the first line is an instance of the famous double orthogonality property in electromagnetism. In our case, it holds with respect to the "natural" $H(\text{curl}; \Omega)$ and $L^2(\Omega)$ inner products. This double orthogonality property is crucial to establish T-coercivity. Before we proceed, we recall a first Weber inequality [34]. There exists $C_W > 0$ such that, for all $\mathbf{v} \in H_0(\text{curl}; \Omega) \cap H(\text{div}; \Omega)$, one has

$$ \|\mathbf{v}\|_{L^2(\Omega)} \le C_W \left\{ \|\mathbf{curl}\,\mathbf{v}\|_{L^2(\Omega)} + \|\text{div}\,\mathbf{v}\|_{L^2(\Omega)} + \sum_{k=1,K} |\langle\mathbf{v} \cdot \mathbf{n}, 1\rangle_{H^{1/2}(\Gamma_k)}| \right\}. $$

It follows that, for all $\mathbf{k} \in K_N^{-}(\Omega)$, one has the bound

$$ \|\mathbf{k}\|_{L^2(\Omega)} \le C_W \|\mathbf{curl}\,\mathbf{k}\|_{L^2(\Omega)}. \qquad (16) $$

**Theorem 3** The form $a$ is T-coercive.
10 Solving magnetostatic systems

Proof Let $(\mathbf{u}, p) \in \mathbb{V}$. Our goal is to find $(\mathbf{v}^*, q^*) \in \mathbb{V}$ that depends linearly (and continuously) on $(\mathbf{u}, p)$, and such that

$$|a((\mathbf{u}, p), (\mathbf{v}^*, q^*))| \geq \alpha^* \|(\mathbf{u}, p)\|_{\mathbb{V}}^2,$$

with $\alpha^* > 0$ that is independent of $(\mathbf{u}, p)$. In this case, defining $T((\mathbf{u}, p)) = (\mathbf{v}^*, q^*)$ yields T-coercivity. To that aim, we decompose $\mathbf{u}$ into $\mathbf{u} = \mathbf{k}_u + \nabla\phi_u$, with $\mathbf{k}_u \in \mathbf{K}_N^-(\Omega)$ and $\phi_u \in H_{\partial\Omega}^1(\Omega)$, and we choose $\mathbf{v}^* = \mathbf{k}_u + \nabla p$, respectively $q^* = \phi_u$. In particular, it holds that $\mathbf{curl}\mathbf{v}^* = \mathbf{curl}\mathbf{k}_u = \mathbf{curl}\mathbf{u}$. We get, using the double orthogonality properties

$$a((\mathbf{u}, p), (\mathbf{v}^*, q^*)) = (\mu^{-1}\mathbf{curl}\mathbf{u}, \mathbf{curl}\mathbf{v}^*)_{L^2(\Omega)} + \gamma(\mathbf{v}^*, \nabla p)_{L^2(\Omega)} + \gamma(\mathbf{u}, \nabla q^*)_{L^2(\Omega)} \\ = (\mu^{-1}\mathbf{curl}\mathbf{k}_u, \mathbf{curl}\mathbf{k}_u)_{L^2(\Omega)} + \gamma\|\nabla p\|_{L^2(\Omega)}^2 + \gamma\|\nabla\phi_u\|_{L^2(\Omega)}^2.$$

By assumption, $0 \in \Theta_{\mu-1}$, with $(\mu^{-1})_- = \mu_-(\mu_+)^{-2}$.<fnref n="2" /> So one finds that

$$\Re(a((\mathbf{u}, p), (\mathbf{v}^*, q^*)) \geq (\mu^{-1})_- \|\mathbf{curl}\mathbf{k}_u\|_{L^2(\Omega)}^2 + \gamma \|\nabla p\|_{L^2(\Omega)}^2 + \gamma \|\nabla\phi_u\|_{L^2(\Omega)}^2 \\ \geq \frac{(\mu^{-1})_-}{1 + C_W^2} \|\mathbf{k}_u\|_{H(\mathbf{curl}; \Omega)}^2 + \gamma \|\nabla p\|_{L^2(\Omega)}^2 + \gamma \|\nabla\phi_u\|_{L^2(\Omega)}^2 \\ \geq \alpha^* \left\{ \|\mathbf{k}_u\|_{H(\mathbf{curl}; \Omega)}^2 + \|\nabla p\|_{L^2(\Omega)}^2 + \|\nabla\phi_u\|_{L^2(\Omega)}^2 \right\} \\ = \alpha^* \left\{ \|\mathbf{k}_u + \nabla\phi_u\|_{H(\mathbf{curl}; \Omega)}^2 + \|\nabla p\|_{L^2(\Omega)}^2 \right\} \\ = \alpha^* \|(\mathbf{u}, p)\|_{\mathbb{V}}^2$$

$$\text{where } \alpha^* = \min \left( (\mu^{-1})_-(1 + C_W^2)^{-1}, \gamma \right) = \min \left( \mu_-(\mu_+)^{-2}(1 + C_W^2)^{-1}, \gamma \right) > 0.$$

Let $(\mathbf{u}, p) \in \mathbb{V}$. Using the above notation, we note that $\mathbf{v}^* = \mathbf{k}_u + \nabla p$ is the orthogonal decomposition (15) of $\mathbf{v}^* \in H_0(\mathbf{curl}; \Omega)$ since $\mathbf{k}_u \in \mathbf{K}_N^-(\Omega)$ and $q^* \in H_{\partial\Omega}^1(\Omega)$, so that $T(\mathbf{v}^*, q^*) = (\mathbf{k}_u + \nabla\phi_u, p) = (\mathbf{u}, p)$. Hence, $T \circ T = I_{\mathbb{V}}$, and T is a bijective operator. This concludes the proof. □

According to the abstract T-coercivity theory, one has the

**Corollary 4** The variational formulation (14) is well-posed.

## 5.4 A perturbed approach

Let T be defined as in the proof of theorem 3: $T(\mathbf{v}, q) = (\mathbf{k}_v + \nabla q, \phi_v)$. Since T is bijective, the variational formulation (14) is equivalent to

$$\left\{ \begin{array}{l} \text{Find } (\mathbf{A}, p) \in \mathbb{V} \text{ such that} \\ \forall (\mathbf{v}, q) \in \mathbb{V}, a((\mathbf{A}, p), T(\mathbf{v}, q)) = \ell(T(\mathbf{v}, q)). \end{array} \right. \qquad (17)$$

$$a((\mathbf{u}, p), (\mathbf{v}, q)) = e^{i\theta}(\mu^{-1}\mathbf{curl}\mathbf{u}, \mathbf{curl}\mathbf{v})_{L^2(\Omega)} + \gamma(\mathbf{v}, \nabla p)_{L^2(\Omega)} + \gamma(\mathbf{u}, \nabla q)_{L^2(\Omega)} \\ \ell((\mathbf{v}, q)) = e^{i\theta}(\mathbf{J}, \mathbf{v})_{L^2(\Omega)},$$

<footnote n="2">² If $0 \notin \Theta_{\mu-1}$, one chooses some $\theta \in \Theta_{\mu-1}$ to solve (14) with the tilted forms</footnote>

<footnote>and the proof can be carried out similarly. Furthermore, one can check easily that, if $0 \notin \Theta_{\mu-1}$, all proofs given hereafter still hold by considering appropriately tilted forms.</footnote>
In the spirit of section 2.3.2 in [3], we note that if (17) can be explicitly discretized, then one has to solve numerically a "simpler" discrete problem than the one that approximates (14).

By design, the form $a' : ((u, p), (v, q)) \mapsto a((u, p), T(v, q))$ is coercive. Using $T(v, q) = (k_v + \nabla q, \phi_v)$ one can compute its expression:

$$a'((u, p), (v, q)) = (\mu^{-1}\mathrm{curl}\,u, \mathrm{curl}\,v)_{L^2(\Omega)} + \gamma(\nabla q, \nabla p)_{L^2(\Omega)} + \gamma(\nabla \phi_u, \nabla \phi_v)_{L^2(\Omega)}.$$

Computationally speaking, the first two terms are explicit, while one has to know the gradient parts of $u$ and $v$ to compute the last term.

Then, we note that, according to (2), the data $J$ is orthogonal to $\nabla[H_{\partial\Omega}^1(\Omega)]$ in $L^2(\Omega)$. So, the expression of $\ell' : (v, q) \mapsto \ell(T(v, q))$ is

$$\begin{aligned}
\ell'((v, q)) &= (J, k_v + \nabla q)_{L^2(\Omega)} \\
&= (J, k_v + \nabla q + \nabla(\phi_v - q))_{L^2(\Omega)} \\
&= (J, v)_{L^2(\Omega)} \\
&= \ell((v, q)).
\end{aligned}$$

Hence one may rewrite (17) as the variational formulation

$$\begin{cases} \text{Find } (A, p) \in \mathbb{V} \text{ such that} \\ \forall (v, q) \in \mathbb{V}, a'((A, p), (v, q)) = \ell((v, q)). \end{cases}$$

Splitting the unknowns, it writes

$$\left\{
\begin{array}{l}
\text{Find } (A, p) \in H_0(\text{curl}; \Omega) \times H_{\partial\Omega}^1(\Omega) \text{ such that} \\
\forall v \in H_0(\text{curl}; \Omega), (\mu^{-1}\text{curl } A, \text{curl } v)_{L^2(\Omega)} \\
\qquad \qquad \qquad \qquad + \gamma(\nabla\phi_u, \nabla\phi_v)_{L^2(\Omega)} = (J, v)_{L^2(\Omega)} \\
\forall q \in H_{\partial\Omega}^1(\Omega), \qquad \gamma(\nabla p, \nabla q)_{L^2(\Omega)} = 0
\end{array}
\right.
$$

For the part in $p$, one has obviously $p = 0$, and we end up with the new variational formulation

$$\left\{
\begin{array}{l}
\text{Find } A \in H_0(\text{curl}; \Omega) \text{ such that} \\
\forall v \in H_0(\text{curl}; \Omega), c(A, v) = (J, v)_{L^2(\Omega)}
\end{array}
\right.
\qquad (18)$$

with

$$c(u, v) = (\mu^{-1}\mathrm{curl}\,u, \mathrm{curl}\,v)_{L^2(\Omega)} + \gamma(\nabla \phi_u, \nabla \phi_v)_{L^2(\Omega)}.$$

From a computational point of view, we note that, in the original variational formulation (14), all terms are explicitly computable, but the form $a$ is not coercive. Moreover, the solution is made up of two parts, a physical part $A$ and an artificial part $p$. On the other hand, in the new variational formulation (18), one term is not explicitly known (second term in the expression of $c$), while the form $c$ is coercive, and only the physical part $A$ of the solution remains. At the discrete level, forgetting for a moment the non-explicit term
12 Solving magnetostatic systems

in c, the cost of linear solvers should be higher in the discrete version of (14) than in the discrete version of (18). Below, we propose a perturbed approach, that allows one to recover explicit terms in the definition of the coercive sesquilinear form, at the expense of solving an inexact problem with a "small" perturbation.

In the new variational formulation (18), there are two options: either one can evaluate simply the second term in the expression of c, that is evaluate simply the gradient part in the Helmholtz decomposition (15) of any field in $H_0(\text{curl}; \Omega)$. Or, one has to modify this second term. We study next the second option.

We observe that the solution $\mathbf{A}$ is independent of the value of the parameter $\gamma$. So, a natural idea is to choose a "small" value of $\gamma$, and to add a perturbation in the order of $\gamma$. How so? Let us introduce

$$c_{\gamma}(\mathbf{u}, \mathbf{v}) = c(\mathbf{u}, \mathbf{v}) + \gamma(\mathbf{k}_{u}, \mathbf{k}_{v})_{L^2(\Omega)}.$$

Using again the orthogonality properties, we find the following expression to the perturbed sesquilinear form

$$c_{\gamma}(\mathbf{u}, \mathbf{v}) = (\mu^{-1}\text{curl } \mathbf{u}, \text{curl } \mathbf{v})_{L^2(\Omega)} + \gamma(\mathbf{u}, \mathbf{v})_{L^2(\Omega)}.$$

*Remark 1* Since by assumption $0 \in \Theta_{\mu^{-1}}$, the form $c_{\gamma}$ is coercive for all $\gamma > 0$. We note that, for $\gamma$ smaller than $(\mu^{-1})_-$, the largest coercivity constant is equal to $\gamma$, while the norm remains bounded by $(\mu^{-1})_+$.

Then, we are solving the perturbed variational formulation

$$\left\{ \begin{array}{l} \text{Find } \mathbf{A}_{\gamma} \in \mathbf{H}_0(\text{curl}; \Omega) \text{ such that} \\ \forall \mathbf{v} \in \mathbf{H}_0(\text{curl}; \Omega), c_{\gamma}(\mathbf{A}_{\gamma}, \mathbf{v}) = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}. \end{array} \right. \qquad (19)$$

We remark that

$$\text{curl} (\mu^{-1} \text{curl} \mathbf{A}_{\gamma}) + \gamma \mathbf{A}_{\gamma} = \mathbf{J} \text{ in } \Omega, \qquad (20)$$

so in general $\mathbf{A}_{\gamma} \neq \mathbf{A}$.

This perturbed variational formulation corresponds to a model that has been considered in [35–37], in a simplified geometrical setting, that is in a simply connected domain, with a connected boundary. Also, the assumptions on the permeability $\mu$ are more restrictive in those Refs, namely it is a real-valued symmetric tensor field.

Below, we study the difference between the exact solution $\mathbf{A}$ and the approximate solution $\mathbf{A}_{\gamma}$. The next proposition is proved in [35, 36] in the "simpler" configuration, this is the reason why we provide a sketched proof.
**Proposition 5** For all $\gamma > 0$, $A_\gamma \in K_N^-(\Omega)$. Furthermore, one has the estimates

$$ \begin{aligned} \|\mathrm{curl}(A_\gamma - A)\|_{L^2(\Omega)} &\le C'\gamma \|J\|_{L^2(\Omega)} \quad \text{and} \\ \|A_\gamma - A\|_{L^2(\Omega)} &\le C' C_W \gamma \|J\|_{L^2(\Omega)}, \end{aligned} \quad (21) $$

with $C' = (C_W)^3 (\mu_+)^4 (\mu_-)^{-2}$.

*Proof* Given $q \in H_{\partial\Omega}^1(\Omega)$, one can choose the test-field $v = \nabla q$ in (19). By assumption on the data, cf. (2), we find $(A_\gamma, \nabla q)_{L^2(\Omega)} = 0$. Hence, $A_\gamma \in K_N^-(\Omega)$.

In addition, according to the Weber inequality (16), it follows that $\|A_\gamma\|_{L^2(\Omega)} \le C_W \|\mathrm{curl} A_\gamma\|_{L^2(\Omega)}$. Thanks to $0 \in \Theta_{\mu^{-1}}$, then using the definition of $c_\gamma$ and finally taking $v = A_\gamma$ in (19), we find that

$$ \begin{aligned} \mu_-(\mu_+)^{-2} \|\mathrm{curl} A_\gamma\|_{L^2(\Omega)}^2 &\le \Re \left( (\mu^{-1} \mathrm{curl} A_\gamma, \mathrm{curl} A_\gamma)_{L^2(\Omega)} \right) \\ &\le \Re (c_\gamma(A_\gamma, A_\gamma)) \\ &\le \|J\|_{L^2(\Omega)} \|A_\gamma\|_{L^2(\Omega)} \\ &\le C_W \|J\|_{L^2(\Omega)} \|\mathrm{curl} A_\gamma\|_{L^2(\Omega)}. \end{aligned} $$

Using the above, it follows that $\|\mathrm{curl} A_\gamma\|_{L^2(\Omega)} \le C_W (\mu_+)^2 (\mu_-)^{-1} \|J\|_{L^2(\Omega)}$, and $\|A_\gamma\|_{L^2(\Omega)} \le (C_W)^2 (\mu_+)^2 (\mu_-)^{-1} \|J\|_{L^2(\Omega)}$.

Comparing (9a) to (19), we find that

$$ \forall v \in H_0(\mathrm{curl}; \Omega), \quad (\mu^{-1}\mathrm{curl}(A_\gamma - A), \mathrm{curl} v)_{L^2(\Omega)} + \gamma(A_\gamma, v)_{L^2(\Omega)} = 0. $$

Using $v = (A_\gamma - A) \in K_N^-(\Omega)$ above, we now find

$$ \begin{aligned} \mu_-(\mu_+)^{-2} \|\mathrm{curl}(A_\gamma - A)\|_{L^2(\Omega)}^2 &\le \Re \left( (\mu^{-1} \mathrm{curl}(A_\gamma - A), \mathrm{curl}(A_\gamma - A))_{L^2(\Omega)} \right) \\ &\le \gamma \|A_\gamma\|_{L^2(\Omega)} \|A_\gamma - A\|_{L^2(\Omega)} \\ &\le C_W \gamma \|A_\gamma\|_{L^2(\Omega)} \|\mathrm{curl}(A_\gamma - A)\|_{L^2(\Omega)}, \end{aligned} $$

so that

$$ \begin{aligned} \|\mathrm{curl}(A_\gamma - A)\|_{L^2(\Omega)} &\le C_W (\mu_+)^2 (\mu_-)^{-1} \gamma \|A_\gamma\|_{L^2(\Omega)} \\ &\le (C_W)^3 (\mu_+)^4 (\mu_-)^{-2} \gamma \|J\|_{L^2(\Omega)}, \end{aligned} $$

i.e. the first inequality in (21). And, according to (16), we conclude that

$$ \|A_\gamma - A\|_{L^2(\Omega)} \le (C_W)^4 (\mu_+)^4 (\mu_-)^{-2} \gamma \|J\|_{L^2(\Omega)}, $$

which is the second inequality in (21). □

## 5.5 Approximation by edge finite elements

Below, we propose to discretize the perturbed problem (19) using edge finite elements. The discrete solution will prove to be a suitable approximation of the vector potential $A$, under appropriate conditions on the perturbation parameter $\gamma$. Recall that we are looking for an approximation of $B_J = \mathrm{curl} A$, so we focus on finding a suitable approximation, and suitable approximation results, for $\mathrm{curl} A$. In this respect, compared to [35–37], the process is simplified, since one is looking for approximation results for both $A$ and $\mathrm{curl} A$ in those Refs.
14 Solving magnetostatic systems

Denoting by $h_K$ the diameter of $K$, each mesh is indexed by the meshsize $h = \max_K h_K$. We choose the Nédélec’s first family of edge finite elements [4, 5] to define finite dimensional subspaces $(V_h)_h$ of $H_0(\text{curl}, \Omega)$. Below, we consider first-order finite elements. In principle, the results can be extended to higher-order finite elements without any difficulty. So, let $\mathcal{R}_1(K)$ be the vector space of polynomials on $K$ defined by

$$ \mathcal{R}_1(K) = \{v \in \mathbf{P}_1(K) : v(x) = a + b \times x, a, b \in \mathbb{C}^3\}. $$

For given $h$, one introduces

$$ V_h = \{v_h \in H_0(\text{curl}, \Omega) : v_{h|K} \in \mathcal{R}_1(K), \forall K \in \mathcal{T}_h \}. $$

We introduce $\Pi_h^{\text{curl}}$, the classical global Nédélec interpolant in $H_0(\text{curl}, \Omega)$ with values in $V_h$ [4, 5]. According Lemma 5.1 of [38]

**Proposition 6** Assume that $v \in \mathbf{PH}^s(\Omega)$ and $\text{curl } v \in \mathbf{PH}^{s'}(\Omega)$ for some $s > 1/2$, $s' > 0$. Then one can define $\Pi_h^{\text{curl}}v$ and, in addition, one has the interpolation result:

$$ \|v - \Pi_h^{\text{curl}} v\|_{H(\text{curl}; \Omega)} \lesssim h^{\min(s,s',1)} \{\|v\|_{\mathbf{PH}^s(\Omega)} + \|\text{curl } v\|_{\mathbf{PH}^{s'}(\Omega)}\}. \quad (22) $$

Since we are interested in the approximation of the curl of $A$, we now refine (22). We introduce $\Pi_h^{\text{div}}$ the classical global Raviart-Thomas interpolation operator in $H_0(\text{div}; \Omega)$ with values in $W_h$ [4, 39], where $(W_h)_h$ are designed with the help of $H(\text{div}; \Omega)$-conforming, first-order finite element spaces. Here,

$$ W_h = \{w_h \in H_0(\text{div}; \Omega) : w_{h|K} \in \mathcal{D}_1(K), \forall K \in \mathcal{T}_h \}. $$

with $\mathcal{D}_1(K)$ the vector space of polynomials on $K$ defined by

$$ \mathcal{D}_1(K) = \{v \in \mathbf{P}_1(K) : v(x) = a + bx, a \in \mathbb{C}^3, b \in \mathbb{C}\}. $$

Using the results of Chapter 5 in [5], we recall that

$$ \forall v \in H_0(\text{curl}; \Omega) \text{ s.t. } \Pi_h^{\text{curl}} v \text{ exists, } \Pi_h^{\text{div}}(\text{curl } v) = \text{curl}(\Pi_h^{\text{curl}} v). $$

Then, applying Lemma 3.3 of [40] (interpolation error for Raviart-Thomas discretization) together with the previous property, we have

**Proposition 7** Assume that $v \in H_0(\text{curl}; \Omega)$ and $\text{curl } v \in \mathbf{PH}^{s'}(\Omega)$ for some $s' > 0$. Then if one can define $\Pi_h^{\text{curl}}v$, one has the interpolation result:

$$ \|\text{curl}(v - \Pi_h^{\text{curl}} v)\|_{L^2(\Omega)} \lesssim h^{\min(s',1)} \|\text{curl } v\|_{\mathbf{PH}^{s'}(\Omega)}. \quad (23) $$
Given a perturbation parameter $\gamma > 0$ and a meshsize $h$, the discrete variational formulation of the perturbed problem (19) is then

$$ \begin{cases} \text{Find } \mathbf{A}_{\gamma}^{h} \in \mathbf{V}_{h} \text{ such that} \\ \forall \mathbf{v}_{h} \in \mathbf{V}_{h}, c_{\gamma}(\mathbf{A}_{\gamma}^{h}, \mathbf{v}_{h}) = (\mathbf{J}, \mathbf{v}_{h})_{L^{2}(\Omega)}. \end{cases} \qquad (24) $$

## 5.6 Numerical analysis

To carry out the numerical analysis, we rely on a variant of the Céa lemma.

**Proposition 8** For all $\gamma > 0$ and all $h$, one has the estimate

$$ \|\mathrm{curl}(\mathbf{A}_{\gamma} - \mathbf{A}_{\gamma}^{h})\|_{L^2(\Omega)} \le \inf_{\mathbf{v}_h \in \mathbf{V}_h} \left\{ \gamma^{1/2} \frac{\mu_+}{\sqrt{2}(\mu_-)^{1/2}} \|\mathbf{A}_{\gamma} - \mathbf{v}_h\|_{L^2(\Omega)} + \frac{\mu_+^2}{\mu_-^2} \|\mathrm{curl}(\mathbf{A}_{\gamma} - \mathbf{v}_h)\|_{L^2(\Omega)} \right\}. \qquad (25) $$

*Proof* For any $\mathbf{v}_h \in \mathbf{V}_h$ and any $\eta > 0$ (for Young's inequality), it holds that

$$ \begin{align*} \mu_- (\mu_+)^{-2} \|\mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h)\|_{L^2(\Omega)} &\le \Re \left( (\mu^{-1} \mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h), \mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h))_{L^2(\Omega)} \right) \\ &= \Re \left( c_\gamma (\mathbf{A}_\gamma - \mathbf{A}_\gamma^h, \mathbf{A}_\gamma - \mathbf{A}_\gamma^h) \right) - \gamma \|\mathbf{A}_\gamma - \mathbf{A}_\gamma^h\|_{L^2(\Omega)}^2 \\ \text{(19) and (24) with } \mathbf{A}_\gamma^h - \mathbf{v}_h \quad &= \Re \left( c_\gamma (\mathbf{A}_\gamma - \mathbf{A}_\gamma^h, \mathbf{A}_\gamma - \mathbf{v}_h) \right) - \gamma \|\mathbf{A}_\gamma - \mathbf{A}_\gamma^h\|_{L^2(\Omega)}^2 \\ &= \Re \left( \gamma (\mathbf{A}_\gamma - \mathbf{A}_\gamma^h, \mathbf{A}_\gamma - \mathbf{v}_h)_{L^2(\Omega)} \right) - \gamma \|\mathbf{A}_\gamma - \mathbf{A}_\gamma^h\|_{L^2(\Omega)}^2 \\ &\quad + \Re \left( (\mu^{-1} \mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h), \mathrm{curl}(\mathbf{A}_\gamma - \mathbf{v}_h))_{L^2(\Omega)} \right) \\ (ab \le a^2 + \frac{1}{4}b^2) \quad &\le \frac{\gamma}{4} \|\mathbf{A}_\gamma - \mathbf{v}_h\|_{L^2(\Omega)}^2 \\ &\quad + \Re \left( (\mu^{-1} \mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h), \mathrm{curl}(\mathbf{A}_\gamma - \mathbf{v}_h))_{L^2(\Omega)} \right) \\ &\le \frac{\gamma}{4} \|\mathbf{A}_\gamma - \mathbf{v}_h\|_{L^2(\Omega)}^2 \\ &\quad + \frac{1}{\mu_-} \|\mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h)\|_{L^2(\Omega)} \|\mathrm{curl}(\mathbf{A}_\gamma - \mathbf{v}_h)\|_{L^2(\Omega)} \\ (ab \le \frac{\eta}{2}a^2 + \frac{1}{2\eta}b^2)) \quad &\le \frac{\gamma}{4} \|\mathbf{A}_\gamma - \mathbf{v}_h\|_{L^2(\Omega)}^2 + \frac{1}{\mu_-} \frac{1}{2\eta} \|\mathrm{curl}(\mathbf{A}_\gamma - \mathbf{v}_h)\|_{L^2(\Omega)}^2 \\ &\quad + \frac{1}{\mu_-} \frac{\eta}{2} \|\mathrm{curl}(\mathbf{A}_\gamma - \mathbf{A}_\gamma^h)\|_{L^2(\Omega)}^2. \end{align*} $$

Then, choosing $\eta = (\mu_-)^2(\mu_+)^{-2}$ above yields

$$ \begin{align*} \|\mathrm{curl}(\mathbf{A}_{\gamma} - \mathbf{A}_{\gamma}^{h})\|_{L^2(\Omega)}^2 &\le \gamma \frac{\mu_{+}^2}{2\mu_{-}} \|\mathbf{A}_{\gamma} - \mathbf{v}_{h}\|_{L^2(\Omega)}^2 + \frac{\mu_{+}^4}{\mu_{-}^4} \|\mathrm{curl}(\mathbf{A}_{\gamma} - \mathbf{v}_{h})\|_{L^2(\Omega)}^2 \\ &\le \left( \gamma^{1/2} \frac{\mu_{+}}{\sqrt{2}(\mu_{-})^{1/2}} \|\mathbf{A}_{\gamma} - \mathbf{v}_{h}\|_{L^2(\Omega)} \right)^2 \end{align*} $$
16 Solving magnetostatic systems

$$+ \frac{\mu_+^2}{\mu_-^2} \| \mathbf{curl} (\mathbf{A}_\gamma - \mathbf{v}_h) \|_{L^2(\Omega)}\bigg)^2,$$

which leads to the bound (25).

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

Remark 2 If the permittivity $\mu$ is a hermitian tensor field, then $c_\gamma$ defines an inner product. Denoting $\| \cdot \|_{\gamma} : \mathbf{v} \mapsto (\gamma \|\mathbf{v}\|_{L^2(\Omega)}^2 + (\mu^{-1} \mathbf{curl} \mathbf{v}, \mathbf{curl} \mathbf{v})_{L^2(\Omega)})^{1/2}$ the associated norm, one has the "straightforward" estimate for all $\gamma > 0$ and all $h$,

$$\| \mathbf{A}_{\gamma} - \mathbf{A}_{\gamma}^{h} \|_{\gamma} \leq \inf_{\mathbf{v}_{h} \in \mathbf{V}_{h}} \| \mathbf{A}_{\gamma} - \mathbf{v}_{h} \|_{\gamma}, \quad (26)$$

and it holds that $\|\mathbf{curl}(\mathbf{A}_{\gamma} - \mathbf{A}_{\gamma}^{h})\|_{L^2(\Omega)} \leq \mu_{+} \inf_{\mathbf{v}_{h} \in \mathbf{V}_{h}} \|\mathbf{A}_{\gamma} - \mathbf{v}_{h}\|_{\gamma}$.

Using (21) and (25), we find

$$\|\mathbf{curl}(\mathbf{A} - \mathbf{A}_{\gamma}^{h})\|_{L^2(\Omega)} \lesssim \inf_{\mathbf{v}_h \in \mathbf{V}_h} \left\{ \gamma^{1/2} \|\mathbf{A}_{\gamma} - \mathbf{v}_h\|_{L^2(\Omega)} + \|\mathbf{curl}(\mathbf{A}_{\gamma} - \mathbf{v}_h)\|_{L^2(\Omega)} \right\} + \gamma \|\mathbf{J}\|_{L^2(\Omega)}. \quad (27)$$

To bound the infimum, we would like to use $\mathbf{v}_h = \Pi_h^{\text{curl}} \mathbf{A}_\gamma$. This is possible if $\mathbf{A}_\gamma$ is sufficiently smooth, in the sense of proposition 6.

Forgetting for the moment the topic of discretization, let us now recall some abstract results regarding a priori smoothness of $\mathbf{A}_\gamma$. First, a classical shift theorem in the domain $\Omega$. We call $\sigma_{Dir} \in ]0, 1]$ the limit regularity exponent for the Laplace problem with Dirichlet boundary condition

Find $u \in H_0^1(\Omega)$ such that $-\Delta u = f$ in $\Omega$

with data $f \in L^2(\Omega)$. According to [41], we know that $u \in H^{3/2}(\Omega)$, with $\|u\|_{H^{3/2}(\Omega)} \lesssim \|f\|_{L^2(\Omega)}$; hence $\sigma_{Dir} \geq 1/2$. In addition, cf. [42, 43]:

* if $\Omega$ is convex or if its boundary is of class $C^2$, then $u \in H^2(\Omega)$, with $\|u\|_{H^2(\Omega)} \lesssim \|f\|_{L^2(\Omega)}$, and in this case $\sigma_{Dir} = 1$;

* if $\Omega$ is a non-convex polyhedron, then there exists $\sigma_{Dir} \in ]1/2, 1[$ such that

- $\forall f \in L^2(\Omega)$, $u \in \mathcal{H}^{1+\sigma_{Dir}}(\Omega) = \bigcap_{0 \le s < \sigma_{Dir}} H^{1+s}(\Omega)$;

- $\exists f \in L^2(\Omega)$ such that $u \notin H^{1+\sigma_{Dir}}(\Omega)$;

- for each $s \in ]0, \sigma_{Dir}$, one has $\|u\|_{H^{1+s}(\Omega)} \lesssim \|f\|_{L^2(\Omega)}$.

We also recall the Birman-Solomyak decomposition, see theorem 4.1 in [44], also called the regular-gradient splitting, see lemma 2.4 in [45]. Any element $u \in H_0(\mathbf{curl}; \Omega)$ can be split as

$$\begin{aligned} u = u_{reg} + \nabla q_u \quad \text{in } \Omega, \text{ with } u_{reg} \in H^1(\Omega), q_u \in H_0^1(\Omega), \text{ and} \\ \|u_{reg}\|_{H^1(\Omega)} + \|q_u\|_{H^1(\Omega)} \lesssim \|u\|_{H(\mathbf{curl}; \Omega)}. \end{aligned} \quad (28)$$

Let us now apply those results to our solution $\mathbf{A}_\gamma$.
**Proposition 9** Let $A_γ$ be the solution to (19). Let $s = 1$ if $σ_{Dir} = 1$, and $s \in ]1/2, σ_{Dir}[$ else. Then, $A_γ \in H^s(\Omega)$, and $||A_γ||_{H^s(\Omega)} \lesssim s ||J||_{L^2(\Omega)}$.

*Proof* Since $A_γ \in K_N^-(\Omega) \subset H_0(\text{curl}; \Omega)$, with the help of the regular-gradient splitting (28), we find that

$$A_γ = A_{reg} + \nabla q_A \text{ in } \Omega, \text{ with } A_{reg} \in H^1(\Omega), q_A \in H_0^1(\Omega), \text{ and} \\ ||A_{reg}||_{H^1(\Omega)} + ||q_A||_{H^1(\Omega)} \lesssim ||A_γ||_{H(\text{curl}; \Omega)}.$$

By construction, $Δq_A = -\text{div } A_{reg}$ in $\Omega$, with $||\text{div } A_{reg}||_{L^2(\Omega)} \le \sqrt{3}||A_{reg}||_{H^1(\Omega)}$. According to the shift theorem, using the definition of the limit regularity exponent $σ_{Dir}$, we find that $\nabla q_A \in H^s(\Omega)$, with $||\nabla q_A||_{H^s(\Omega)} \lesssim s ||A_{reg}||_{H^1(\Omega)}$. Hence, $A_γ \in H^s(\Omega)$, with the estimate $||A_γ||_{H^s(\Omega)} \lesssim s ||A_γ||_{H(\text{curl}; \Omega)}$, and the first Weber inequality yields $||A_γ||_{H^s(\Omega)} \lesssim s ||\text{curl } A_γ||_{L^2(\Omega)}$. Finally, we recall that $||\text{curl } A_γ||_{L^2(\Omega)} \lesssim ||J||_{L^2(\Omega)}$ (proof of proposition 5), so the conclusion follows. □

Then, if one changes the Dirichlet boundary condition to a Neumann boundary condition, one has similar results. The results that we invoke below now involve the tensor $μ$ (in the Dirichlet case, see above, we considered a unit tensor). There holds a second Weber inequality (cf. theorem 2.11 in [18]), namely there exists $C'_W > 0$ such that, for all $v \in H(\text{curl}; \Omega) \cap H_0(\text{div } μ; \Omega)$, one has

$$||v||_{L^2(\Omega)} \le C'_W \left\{ ||\text{curl } v||_{L^2(\Omega)} + ||\text{div } \mu v||_{L^2(\Omega)} + \sum_{i=1, I} |\langle \mu v \cdot n, 1 \rangle_{H^{1/2}(\Sigma_i)}| \right\}.$$

Next, if one considers that $\Omega$ is a polyhedron, one has a second shift theorem, which we state directly for elements of $H(\text{curl}; \Omega) \cap H_0(\text{div } μ; \Omega)$, cf. [46, 47]. It is derived thanks to another regular-gradient splitting, which is omitted. Let $v \in H(\text{curl}; \Omega) \cap H_0(\text{div } μ; \Omega)$. There exists $σ_{Neu}(μ) \in ]0, 1]$ such that

$$H(\text{curl}; \Omega) \cap H_0(\text{div } μ; \Omega) \subset \underline{PH}^{\sigma_{Neu}(μ)}(\Omega) = \bigcap_{0 \le s' < \sigma_{Neu}(μ)} \underline{PH}^{s'}(\Omega), (29)$$

with continuous imbedding. Applying the two results together, one concludes that, for every $s' \in ]0, σ_{Neu}(μ)]$, one has

$$||v||_{PH^{s'}(\Omega)} \lesssim s' \left\{ ||\text{curl } v||_{L^2(\Omega)} + ||\text{div } \mu v||_{L^2(\Omega)} + \sum_{i=1, I} |\langle \mu v \cdot n, 1 \rangle_{H^{1/2}(\Sigma_i)}| \right\}.$$

If $σ_{Neu}(μ) = 1$, the result also holds for $s' = 1$.

We refer for instance to section 7 in [46], section 4.5 in [48], section 3.5 in [47] or section 3.3 in [18] and References therein for a detailed discussion, and possible values of the limit regularity exponent $σ_{Neu}(μ)$. In particular, if $∂Ω$
18 Solving magnetostatic systems

is of class $\mathcal{C}^2$ and if $\mu \in (\mathcal{C}^1(\bar{\Omega}))^{3 \times 3}$, then $\sigma_{Neu}(\mu) = 1$.
Regarding the curl of $A_\gamma$, we find the results below.

**Proposition 10** Let $A_\gamma$ be the solution to (19). Let $\mathbf{s}' = 1$ if $\sigma_{Neu}(\mu) = 1$, and $\mathbf{s}' \in [0, \sigma_{Neu}(\mu)]$ [else. Then, curl $A_\gamma \in PH^{s'}(\Omega)$, and $\|\text{curl } A_\gamma\|_{PH^{s'}(\Omega)} \lesssim_{\mathbf{s}'} (1+\gamma) \|\mathbf{J}\|_{L^2(\Omega)}$.

*Proof* The field $H_\gamma = \mu^{-1}\text{curl } A_\gamma$ is such that $H_\gamma \in H(\text{curl}; \Omega)$ with curl $H_\gamma = \mathbf{J} - \gamma A_\gamma$ in $\Omega$ according to (20). Also, $H_\gamma \in H_0(\text{div } \mu; \Omega)$, with $\text{div}(\mu H_\gamma) = 0$ in $\Omega$. Using the second shift theorem for $H_\gamma$, we find $H_\gamma \in PH^{s'}(\Omega)$. Next, because $\mu H_\gamma \in \text{curl } [H_0(\text{curl}; \Omega)]$, one has always $\langle \mu H_\gamma \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Sigma_i)} = 0$ for $i=1, I$ (cf. remark 3.5.2. in [8]). So we have $\|H_\gamma\|_{PH^{s'}(\Omega)} \lesssim_{\mathbf{s}'} \|\text{curl } H_\gamma\|_{L^2(\Omega)}$.
Going back to curl $H_\gamma = \mathbf{J} - \gamma A_\gamma$ in $\Omega$, we infer

$$ \|H_{\gamma}\|_{PH^{s'}(\Omega)} \lesssim_{\mathbf{s}'} \{ \|J\|_{L^2(\Omega)} + \gamma \|A_{\gamma}\|_{L^2(\Omega)} \}. $$

Recall that the partition $(\Omega_p^\mu)_{p=1,P}$ is chosen so that $\mu|_{\Omega_p^\mu} \in (W^{1,\infty}(\Omega_p^\mu))^{3 \times 3}$ for $p=1, P$. Because multiplying by $\mu|_{\Omega_p^\mu}$ is stable in $H^{s'}(\Omega_p^\mu)$ for $p=1, P$ [49], we find that curl $A_\gamma = \mu H_\gamma \in PH^{s'}(\Omega)$, with

$$ \|\text{curl } A_{\gamma}\|_{PH^{s'}(\Omega)} \lesssim_{\mathbf{s}'} \{ \|\mathbf{J}\|_{L^2(\Omega)} + \gamma \|A_{\gamma}\|_{L^2(\Omega)} \}. $$

To conclude, we recall that $\|A_\gamma\|_{L^2(\Omega)} \lesssim \|\mathbf{J}\|_{L^2(\Omega)}$ (proof of proposition 5). $\square$

As announced, one may use $v_h = \Pi_h^{\text{curl}} A_\gamma$ in (27).

**Theorem 11** For $\gamma > 0$ and $h > 0$, let $A_\gamma^h$ be the solution to (24).
Then, for $\mathbf{s} = 1$ if $\sigma_{Dir} = 1$, and $\mathbf{s} \in [1/2, \sigma_{Dir}]$ [else, and for $\mathbf{s}' = 1$ if $\sigma_{Neu}(\mu) = 1$, and $\mathbf{s}' \in [0, \sigma_{Neu}(\mu)]$ [else, one has the error estimate

$$ \|\text{curl}(\mathbf{A} - \mathbf{A}_{\gamma}^{h})\|_{L^2(\Omega)} \lesssim_{\mathbf{s}, \mathbf{s}'} (\gamma + \gamma^{1/2}(1+\gamma)h^{\min(\mathbf{s}, \mathbf{s}')} + (1+\gamma)h^{\mathbf{s}'}) \|\mathbf{J}\|_{L^2(\Omega)}. \quad (30) $$

*Proof* Let us bound the infimum in (27). Combining the interpolation errors (22) and (23) with the last two propositions, we find

$$ \begin{aligned} \inf_{v_h \in V_h} & \{ \gamma^{1/2} \|A_\gamma - v_h\|_{L^2(\Omega)} + \|\text{curl}(A_\gamma - v_h)\|_{L^2(\Omega)} \} \\ & \lesssim \gamma^{1/2} \|A_\gamma - \Pi_h^{\text{curl}} A_\gamma\|_{L^2(\Omega)} + \|\text{curl}(A_\gamma - \Pi_h^{\text{curl}} A_\gamma)\|_{L^2(\Omega)} \\ & \lesssim_{\mathbf{s}, \mathbf{s}'} \gamma^{1/2} (1+\gamma)h^{\min(\mathbf{s}, \mathbf{s}')} \|\mathbf{J}\|_{L^2(\Omega)} + (1+\gamma)h^{\mathbf{s}'} \|\mathbf{J}\|_{L^2(\Omega)}. \end{aligned} $$

Using (27) leads to the claim. $\square$
To conclude this part, we now consider that $\gamma$ is a function of the meshsize $h$ (with strictly positive values), that is $\gamma = \gamma(h)$, and we solve

$$ \begin{cases} \text{Find } \mathbf{A}_h \in \mathbf{V}_h \text{ such that} \\ \forall \mathbf{v}_h \in \mathbf{V}_h, c_{\gamma(h)}(\mathbf{A}_h, \mathbf{v}_h) = (\mathbf{J}, \mathbf{v}_h)_{L^2(\Omega)}. \end{cases} \quad (31) $$

**Corollary 12** Assume that $h \in [0, 1]$. Let $\gamma = \gamma(h)$ be such that $\gamma(h) \lesssim h^{\sigma_{Neu}(\mu)}$. For $s' = 1$ if $\sigma_{Neu}(\mu) = 1$, and $s' \in [0, \sigma_{Neu}(\mu)]$ [else, one has the error estimate

$$ \|\mathrm{curl}(\mathbf{A} - \mathbf{A}_h)\|_{L^2(\Omega)} \lesssim_{s'} h^{s'} \|\mathbf{J}\|_{L^2(\Omega)}. \quad (32) $$

*Proof* Since in theorem 11 one has $s > 1/2$, it holds $h^{\min(s,s')} \le h^{\min(1/2,s')}$. Hence from (30) one gets the bound:

$$ \|\mathrm{curl}(\mathbf{A} - \mathbf{A}_{\gamma}^{h})\|_{L^2(\Omega)} \lesssim_{s'} (\gamma + \gamma^{1/2}(1+\gamma)h^{\min(1/2,s')}) + (1+\gamma)h^{s'} \|\mathbf{J}\|_{L^2(\Omega)}. $$

Let us find the dominant term between the parentheses. If $s' < 1/2$,

$$ \gamma(h) + \gamma(h)^{1/2}(1 + \gamma(h))h^{\min(1/2, s')} \lesssim h^{\sigma_{Neu}(\mu)} + h^{\sigma_{Neu}(\mu)/2 + s'} = O(h^{s'}). $$

If $s' \ge 1/2$, taking into account that $\sigma_{Neu}(\mu)/2 + 1/2 > \sigma_{Neu}(\mu)$, one has

$$ \gamma(h) + \gamma(h)^{1/2}(1 + \gamma(h))h^{\min(1/2, s')} \lesssim h^{\sigma_{Neu}(\mu)} + h^{\sigma_{Neu}(\mu)/2 + 1/2} = O(h^{\sigma_{Neu}(\mu)}). $$

Since $\gamma(h)h^{s'} = o(h^{s'})$, the dominant term is $h^{s'}$. $\square$

Hence, choosing $\gamma(h) \lesssim h^{\sigma_{Neu}(\mu)}$ and going back to the magnetic field $\mathbf{B}_J$, it is guaranteed that

$$ \|\mathbf{B}_J - \mathrm{curl} \, \mathbf{A}_h\|_{H(\mathrm{div}; \Omega)} \lesssim_{s'} h^{s'} \|\mathbf{J}\|_{L^2(\Omega)}. $$

According to the definition of $s'$, one concludes that the a priori error is equal to $O(h)$ if $\sigma_{Neu}(\mu) = 1$, and almost like $O(h^{\sigma_{Neu}(\mu)})$ if $\sigma_{Neu}(\mu) < 1$, since the result is $O(h^{s'})$ for all $s' < \sigma_{Neu}(\mu)$ in the latter case.

Because the form $c_{\gamma(h)}$ is coercive, the linear system equivalent to (31) is expected to be simpler to solve than the one resulting from the discretization of the saddle-point variational formulation (10) or (14). Moreover, according to remark 1, it is advised to take the "largest" possible values of $\gamma(h)$ to improve the conditioning of the matrix. In other words, one should take $\gamma(h) \approx h^{\sigma_{Neu}(\mu)}$, see corollary 12.

One can find a number of numerical experiments in [35–37] to support the claims, and to show the robustness of the approach.
# 6 A conforming variational formulation in a weighted space

We review here an approach that has been introduced at the turn of the millenium by Costabel and Dauge, see [6]. It relies on taking the divergence explicitly into account in the measure of the vector potentials. Below, we focus on the original version of Costabel and Dauge. See [50–56] for applications in electromagnetism in a homogeneous material. It is especially of interest when the domain $\Omega$ is a non-convex polyhedron$^3$, so we focus on this geometrical configuration from now on. This approach is aimed at designing a mathematical framework that can be efficiently approximated numerically by continuous vector finite elements.$^4$ In most of the above cited References, the problem is solved in a simplified geometrical setting, and, the assumptions on the permeability $\mu$ are more restrictive.

## 6.1 Reminders about the WRM

When $\Omega$ is non-convex, following [6], the divergence of the fields is evaluated in a weighted $L^2$ space. Denote by $E$ the non-empty set of reentrant edges of $\partial\Omega$. Let $d$ be the distance to $E$: $d(x) = \text{dist}(x, \cup_{e\in E}\bar{e})$. Consider $w_\alpha$ a smooth non-negative function of $x$, that depends on a real parameter $\alpha$. The (simplified) weight $w_\alpha$ is chosen to behave locally as $d^\alpha$ in the neighborhood of reentrant edges and corners, and is bounded above and below by a strictly positive constant outside a fixed neighborhood of $E$. Let:

$$L_{\alpha}^{2}(\Omega) = \{q \in L_{\text{loc}}^{2}(\Omega) \text{ such that } w_{\alpha} q \in L^{2}(\Omega)\},$$

endowed with the natural norm $\|\cdot\|_{L_{\alpha}^2(\Omega)}: q \mapsto \|w_{\alpha} q\|_{L^2(\Omega)}$, and associated scalar product $(\cdot, \cdot)_{L_{\alpha}^2(\Omega)}$. Then, we introduce a subspace of $H_0(\text{curl}, \Omega)$:

$$X_N^\alpha(\Omega) = \{v \in H_0(\text{curl}, \Omega) \text{ such that } \text{div } v \in L_\alpha^2(\Omega)\}.$$

It is endowed with the natural norm

$$v \mapsto \left( \|v\|_{H(\text{curl}, \Omega)}^2 + \|\text{div } v\|_{L_{\alpha}^2(\Omega)}^2 \right)^{1/2}.$$

$$X_N(\Omega) = \{v \in H_0(\text{curl}; \Omega) \text{ s.t. } \text{div } v \in L^2(\Omega)\}$$

that is equivalent to the magnetostatic system (4), simply by taking the weight $w_\alpha$ introduced afterwards equal to 1 everywhere. It is endowed with the natural norm $v \mapsto (\|v\|_{H(\text{curl}; \Omega)}^2 + \|\text{div } v\|_{L^2(\Omega)}^2)^{1/2}$. Then, "classical" Lagrange finite elements can be used to approximate this variational formulation. We refer e.g. to Appendix B in [54] for details and elements of numerical analysis. Convergence stems from the fact that, when $\Omega$ is a convex polyhedron, $X_N(\Omega) \cap H^1(\Omega)$ is dense in $X_N(\Omega)$. See also section 8.2.2.B in [6], with another choice of measure of the divergence, similar but not identical to the one we present here.

<footnote>$^3$Non-convexity always occurs when the boundary $\partial\Omega$ is not connected, that is $K \ge 1$.</footnote>

<footnote>$^4$In a convex polyhedron $\Omega$, one can build a variational formulation in the "natural" function space</footnote>
Introducing the value $\alpha_{min} = 1 - \sigma_{Dir} \in ]0, 1/2[, $ one has that for all $\alpha \in ]\alpha_{min}, 1[, \mathbf{X}_{N}^{\alpha}(\Omega) \cap \mathbf{H}^{1}(\Omega)$ is dense in $\mathbf{X}_{N}^{\alpha}(\Omega)$ according to theorem 5.1 of [6]. This result is paramount to obtain the basic approximability property for a conforming finite element method in $\mathbf{X}_{N}^{\alpha}(\Omega)$.<fnref n="5" /> So, we assume from now on that

$$ \alpha \in ]\alpha_{min}, 1[. $$

Among the noticeable properties proven in [6], we observe that $H_0^1(\Omega) \subset L_{-\alpha}^2(\Omega)$ (with dense embedding), see bottom of page 249 of loco citato. For $k=1, K$, let $q_k$ be characterized by

$$ \text{Find } q_k \in H_{\partial\Omega}^1(\Omega) \text{ such that } -\Delta q_k = 0 \text{ in } \Omega, q_k|_{\Gamma_m} = \delta_{km} \text{ for } m=1, K. \quad (33) $$

By construction, one has the orthogonal decomposition

$$ H_{\partial\Omega}^{1}(\Omega) = H_{0}^{1}(\Omega) \oplus \operatorname{span}_{\ell=1,K}(q_{\ell}). \quad (34) $$

One can check that 1 belongs to $L_{-\alpha}^2(\Omega)$ so, for $k=1, K$, we also have that $q_k \in L_{-\alpha}^2(\Omega)$. As a consequence, the fluxes $((\mathbf{v} \cdot \mathbf{n}, 1)_{H^{1/2}(\Gamma_{\ell})})_{\ell=1,K}$ are meaningful for all $\mathbf{v} \in \mathbf{X}_{N}^{\alpha}(\Omega)$, and equal to $((\mathbf{v} \cdot \mathbf{n}, q_{\ell})_{H^{1/2}(\partial\Omega)})_{\ell=1,K}$.

Then, following for instance [8, Proof of Theorem 3.4.3] and using the embedding of $L_{\alpha}^2(\Omega)$ in $H^{-1}(\Omega)$, one can prove that

$$ \|\cdot\|_{\mathbf{X}_{N}^{\alpha}(\Omega)} : \mathbf{v} \mapsto \left( \|\mathrm{curl} \mathbf{v}\|_{L^2(\Omega)}^2 + \|\mathrm{div} \mathbf{v}\|_{L^2_{\alpha}(\Omega)}^2 + \sum_{\ell=1,K} |\langle \mathbf{v} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})}|^2 \right)^{1/2} $$

defines a norm, that is equivalent to the natural norm, in $\mathbf{X}_{N}^{\alpha}(\Omega)$. We let $((\cdot, \cdot))_{\mathbf{X}_{N}^{\alpha}(\Omega)}$ be the associated scalar product.

Going back to the magnetic system (4), it is obvious that its divergence-free solution $\mathbf{A}$ automatically belongs to $\mathbf{X}_{N}^{\alpha}(\Omega)$. Let $\mathbf{v} \in \mathbf{X}_{N}^{\alpha}(\Omega)$. Then, adding the contributions

$$ \begin{aligned} (\mu^{-1}\mathrm{curl} \mathbf{A}, \mathrm{curl} \mathbf{v})_{L^2(\Omega)} &= (\mathrm{curl} \mu^{-1}\mathrm{curl} \mathbf{A}, \mathbf{v})_{L^2(\Omega)} = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}; \\ (\mathrm{div} \mathbf{A}, \mathrm{div} \mathbf{v})_{L^2_{\alpha}(\Omega)} &= 0; \quad \sum_{\ell=1,K} \langle \mathbf{A} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})} \overline{\langle \mathbf{v} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})}} = 0; \end{aligned} $$

we find that $\mathbf{A}$ is governed by the variational formulation

$$ \left\{ \begin{array}{ll} \text{Find } \mathbf{A} \in \mathbf{X}_{N}^{\alpha}(\Omega) \text{ such that} & \\ \forall \mathbf{v} \in \mathbf{X}_{N}^{\alpha}(\Omega), \quad (\mu^{-1}\mathrm{curl} \mathbf{A}, \mathrm{curl} \mathbf{v})_{L^2(\Omega)} + (\mathrm{div} \mathbf{A}, \mathrm{div} \mathbf{v})_{L^2_{\alpha}(\Omega)} & (35) \\ \phantom{\forall \mathbf{v} \in \mathbf{X}_{N}^{\alpha}(\Omega), \quad (\mu^{-1}\mathrm{curl} \mathbf{A}, \mathrm{curl} \mathbf{v})_{L^2(\Omega)} + (\mathrm{div} \mathbf{A}, \mathrm{div} \mathbf{v})_{L^2_{\alpha}(\Omega)} } + \sum_{\ell=1,K} \langle \mathbf{A} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})} \overline{\langle \mathbf{v} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})}} = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}. & \end{array} \right. $$

<footnote n="5">⁵On the other hand, when $\Omega$ is a non-convex polyhedron, according e.g. to [57], $\mathbf{X}_N(\Omega) \cap \mathbf{H}^1(\Omega)$ is not dense in $\mathbf{X}_N(\Omega)$. Hence, a conforming finite element method in $\mathbf{X}_N(\Omega)$ fails to satisfy the basic approximability property, so it can not be used to approximate those fields that do not belong to $\mathbf{X}_N(\Omega) \cap \mathbf{H}^1(\Omega)$.</footnote>
## 22 Solving magnetostatic systems

In the spirit of [50], this is an augmented variational formulation.

Due to the assumption $0 \in \Theta_{\mu^{-1}}$ (otherwise, one can proceed similarly to footnote² page 10), this variational formulation (35) falls within the framework of the Lax-Milgram theorem. As a matter of fact, the sesquilinear form

$$a_{\alpha}(\mathbf{u}, \mathbf{v}) = (\mu^{-1}\mathbf{curl}\mathbf{u}, \mathbf{curl}\mathbf{v})_{L^2(\Omega)} + (\mathrm{div}\,\mathbf{u}, \mathrm{div}\,\mathbf{v})_{L^2_{\alpha}(\Omega)} \\ + \sum_{\ell=1,K} \langle \mathbf{u} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})} \overline{\langle \mathbf{v} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})}}$$

is obviously continuous, and coercive, on $X_N^\alpha(\Omega) \times X_N^\alpha(\Omega)$.

On the other hand, starting from the augmented variational formulation (35), one can prove that its solution $\mathbf{A}$ is governed by (4). To start with, we know that $\mathbf{A} \in H_0(\mathbf{curl}, \Omega)$. Next, we take appropriate test functions in (35). Below, we make use of a triple orthogonality property, which allows to deal with the three terms defining the form $a_\alpha$ separately.

Note that $\mathbf{v}_\beta = \sum_{k=1,K} \beta_k \nabla q_k$, where $(\beta_k)_{k=1,K} \in \mathbb{C}^K$ and $(q_k)_{k=1,K}$ is given by (33), belongs to $X_N^\alpha(\Omega)$. Using it as a test function in (35), we find that

$$\sum_{\ell=1,K} \langle \mathbf{A} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})} \overline{\langle \mathbf{v}_{\beta} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})}} = 0$$

after integrating by parts the right-hand side, thanks to assumption (2). On the other hand, the mapping

$$(\beta_k)_{k=1,K} \mapsto (\langle \mathbf{v}_{\beta} \cdot \mathbf{n}, 1 \rangle_{H^{1/2}(\Gamma_{\ell})})_{k=1,K}$$

is onto. This classical result stems from the fact that the capacitance matrix, with entries $((\nabla q_{\ell}, \nabla q_k)_{L^2(\Omega)})_{1 \le k, \ell \le K}$, is invertible (see e.g. corollary 3.3.8 in [8]). As a consequence, the fluxes $((\mathbf{A} \cdot \mathbf{n}, 1)_{H^{1/2}(\Gamma_k)})_{k=1,K}$ all vanish. Going back to (35), we find that $\mathbf{A}$ is such that, for all $\mathbf{v} \in X_N^\alpha(\Omega)$,

$$(\mu^{-1}\mathbf{curl}\mathbf{A}, \mathbf{curl}\mathbf{v})_{L^2(\Omega)} + (\mathrm{div}\,\mathbf{A}, \mathrm{div}\,\mathbf{v})_{L^2_{\alpha}(\Omega)} = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}.$$

By definition, $\mathrm{div}\,\mathbf{A} \in L^2_{\alpha}(\Omega) \subset H^{-1}(\Omega)$, so there exists a unique $q_A \in H_0^1(\Omega)$ such that $-\Delta q_A = \mathrm{div}\,\mathbf{A}$. The test function $\mathbf{v}_A = \nabla q_A$ can be used above, which now yields

$$\|\mathrm{div}\,\mathbf{A}\|_{L^2_{\alpha}(\Omega)}^2 = 0,$$

after integrating by parts the right-hand side, thanks again to assumption (2). Therefore, $\mathbf{A}$ is such that, for all $\mathbf{v} \in X_N^\alpha(\Omega)$,

$$(\mu^{-1}\mathbf{curl}\mathbf{A}, \mathbf{curl}\mathbf{v})_{L^2(\Omega)} = (\mathbf{J}, \mathbf{v})_{L^2(\Omega)}.$$

Finally, taking $\mathbf{v} \in D(\Omega)$, we conclude that $\mathbf{curl}(\mu^{-1}\mathbf{curl}\mathbf{A}) = \mathbf{J}$, and $\mathbf{A}$ solves the magnetic system (4) as claimed.
**Theorem 13** Let the data **J** fulfill assumption (2). Then, the vector potential **A** solves (4) if, and only if, **A** solves (35).

Note that the exact solution **A** is recovered exactly for all values of $\alpha \in ]\alpha_{min}, 1[$.

## 6.2 Approximation by continuous vector finite elements

Among others, the discretization of problem (4) via the solution of the augmented variational formulation (35) was originally studied in [6], together with the numerical analysis of the resulting method. We report below the main steps of the analysis which has been carried out there. It starts with the study of singularities, cf. sections 6 and 7 in loc. cit..

*Remark 3* Within the taxonomy of [58, section 6], one can prove that **A** contains only singularities of type 1 and type 2, and no singularity of type 3, because $\mathbf{J} \in \mathbf{L}^2(\Omega)$ and $\mathrm{div} \, \mathbf{A} = 0$ (see loc. cit., page 259).

A solution **A** to (4) can be split into regular and gradient parts as follows:

$$ \mathbf{A} = \mathbf{A}_{\mathrm{reg}} + \nabla \phi_A, \quad \text{with} \quad \left\{ \begin{array}{l} \mathbf{A}_{\mathrm{reg}} \in \underline{\mathbf{P}}\underline{\mathbf{H}}^{1+\sigma_{\mathrm{Neu}}(\mu)}(\Omega) \cap \mathbf{H}_0(\mathrm{curl}; \Omega) \\ \phi_A \in \underline{\mathbf{H}}^{1+\sigma_{\mathrm{Dir}}}(\Omega) \cap \mathbf{H}_{\partial\Omega}^1(\Omega) \end{array} \right. . \quad (36) $$

In addition, there exists $C_\alpha > 0$ (independent of **A**) such that

$$ \|\mathbf{A}\|_{\mathbf{X}_N^\alpha(\Omega)} \leq C_\alpha (\|\mathbf{A}_{\mathrm{reg}}\|_{\mathbf{H}^1(\Omega)} + \|\phi_A\|_{V_\alpha^2(\Omega)}), $$

where $V_\alpha^2(\Omega)$ is a weighted Sobolev space à la Nazarov-Plamenevski (cf. loc. cit. sections 4 and 5). Next, one uses this regular-gradient splitting (36) to obtain error estimates; precisely, the error on $\mathbf{A}_{\mathrm{reg}}$ in $\mathbf{H}^1(\Omega)$, and the error on $\phi_A$ in $V_\alpha^2(\Omega)$.

Given $k \ge 1$, we choose the continuous Lagrange elements of order $k$ to define finite dimensional subspaces $(\mathbf{X}_{h,k})_h$ of $\mathbf{X}_N^\alpha(\Omega)$. Given $h$, one introduces

$$ \mathbf{X}_{h,k} = \{ \mathbf{v}_h \in \mathbf{H}^1(\Omega) \cap \mathbf{H}_0(\mathrm{curl}; \Omega) : \mathbf{v}_{h|K} \in \mathcal{P}_k(K), \forall K \in \mathcal{T}_h \}. $$

The discrete variational formulation of the augmented variational formulation (35) is then

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{A}_{\alpha}^{h} \in \mathbf{X}_{h,k} \text{ such that} \\ \forall \mathbf{v}_h \in \mathbf{X}_{h,k}, a_{\alpha}(\mathbf{A}_{\alpha}^{h}, \mathbf{v}_h) = (\mathbf{J}, \mathbf{v}_h)_{\mathbf{L}^2(\Omega)}. \end{array} \right. \quad (37) $$

Observe that the discrete solution $\mathbf{A}_{\alpha}^{h}$ depends on both $h$ and $\alpha$.
24 Solving magnetostatic systems

Classically, $(\mathbf{X}_{h,k})_h$ exhibits approximability properties for the regular part:

$$ egin{cases} 
\text{For } \mathbf{s}' = 1 \text{ if } \sigma_{\text{Neu}}(\mu) = 1, \text{ and all } \mathbf{s}' \in ]0, \sigma_{\text{Neu}}(\mu)] \text{ else} \\
\forall \mathbf{v} \in \underline{\mathbf{P}\mathbf{H}}^{1+\sigma_{\text{Neu}}(\mu)}(\Omega) \cap \mathbf{H}_0(\text{curl}; \Omega), \exists c_{v,\mathbf{s}'}, \forall h, \exists \mathbf{v}_h \in \mathbf{X}_{h,k}, \\
\|\mathbf{v} - \mathbf{v}_h\|_{\mathbf{H}^1(\Omega)} \le c_{v,\mathbf{s}'}h^{\mathbf{s}'}. 
\end{cases} 
\quad (38) $$

In addition, one must have good approximability properties for the gradient part: namely, one needs that there exists a family of spaces $(\Phi_h)_h$ of scalar fields such that

$$ \nabla[\Phi_h] \subset \mathbf{X}_{h,k}, \forall h, \quad (39) $$

with<fnref n="6" />

$$ \left\{ \begin{array}{l} \forall \varepsilon > 0, \quad \forall \phi \in \underline{H}^{1+\sigma_{\text{Dir}}}(\Omega) \cap H_{\partial\Omega}^{1}(\Omega), \quad \exists c_{\phi,\varepsilon}, \forall h, \quad \exists \phi_h \in \Phi_h, \\ \qquad \|\phi - \phi_h\|_{V_{\alpha}^2(\Omega)} \le c_{\phi,\varepsilon} h^{\alpha - \alpha_{\text{min}} - \varepsilon}. \end{array} \right. \quad (40) $$

A natural question is to check in which case the assumption (39)-(40) holds. Below, we report some results in 2D and 3D domains. The issue was first discussed by Costabel and Dauge in [6]. They observed that, in a 2D polygonal domain, one can choose Hsieh-Clough-Tocher finite elements [59] to define the discrete spaces $(\Phi_h)_h$. Indeed, the discrete fields belong to $C^1(\bar{\Omega})$ and are piecewise $P_3$ on a subdivision of triangles, so their gradients are naturally in $C^0(\bar{\Omega})$ and are piecewise $P_2$. Moreover, one can choose those fields to vanish on the boundary, to enforce the condition on the tangential trace. This ensures that the assumption (39)-(40) holds in a 2D domain at least for $k=2$.

Browsing the classical literature, we note that, according to [60, section 46], there exist similar finite elements in 2D for $k \ge 3$, and also in 3D: the former can be found in [61], and the latter in [62]. More recent results can be found in [63–65] and Refs. therein, including the case $k=1$, and are proved by considering again appropriate subdivisions of tetrahedra and/or triangles.

## 6.3 Numerical analysis of the augmented variational formulation

The numerical analysis is straightforward, because one can apply Céa's lemma. Using the approximability properties of section 6.2, we conclude that the error estimates below hold.

**Theorem 14** For $\alpha \in ]\alpha_{\text{min}}, 1[$, let $\tau_\alpha = \min(\sigma_{\text{Neu}}(\mu), \sigma_{\text{Dir}} + (\alpha - 1))$. For $h > 0$, let $\mathbf{A}_\alpha^h$ be the solution to (37). Then, for all $\varepsilon > 0$, there exists $c_{\alpha, \varepsilon}$ such that for all $h$:

$$ \|\mathbf{A} - \mathbf{A}_\alpha^h\|_{\mathbf{X}_N^\alpha(\Omega)} \le c_{\alpha, \varepsilon} h^{\tau_\alpha - \varepsilon}, \quad (41) $$

<footnote n="6">⁶Recall that $\alpha_{\text{min}} = 1 - \sigma_{\text{Dir}} \in ]0, 1/2[$, so the right-hand side of (40) is equal to $c_{\phi, \varepsilon} h^{\sigma_{\text{Dir}}+(\alpha-1)-\varepsilon}$.</footnote>
Going back to the magnetic field $B_J$, that it is now guaranteed that

$$ \|B_J - \text{curl } A_{\alpha}^h\|_{H(\text{div}; \Omega)} \le c_{J,\varepsilon} h^{\tau_{\alpha}-\varepsilon}. $$

Within the framework of the WRM, we can take any $\alpha \in ]\alpha_{min}, 1[$ so, if one picks $\alpha$ "close to 1", the a priori error is almost like $O(h^{\tau_{max}})$, where $\tau_{max} = \min(\sigma_{Neu}(\mu), \sigma_{Dir})$.

We note that the linear system equivalent to (37) is expected to be easy to solve because the form $a_{\alpha}$ is coercive. In 2D and 3D domains, one can find numerical experiments [51, 53, 66, 67] to support the claims, and to show the robustness of the approach. In particular, in [66, 67], a very coarse approximation of the weight $w_{\alpha}$ is chosen, namely a discrete weight which is equal to 0 in simplices located near the reentrant edges, and equal to 1 in all other simplices.

# 7 Concluding remarks

Going back to the finite dimensional part of the magnetic induction made of curl- and div-free fields, we note that, by definition, the solution $\dot{p}$ of (6) is such that $\nabla\dot{p}$ belongs to $H(\text{curl}; \Omega) \cap H_0(\text{div}\,\mu; \Omega)$. Thanks to the a priori regularity result (29), one gets that the a priori error between the exact and the discrete solutions is equal to $O(h)$ if $\sigma_{Neu}(\mu) = 1$, and behaves almost like $O(h^{\sigma_{Neu}(\mu)})$ if $\sigma_{Neu}(\mu) < 1$. And one has the same a priori error for $B_c = \mu\nabla\dot{p}$ and the corresponding discrete field. More precisely, the plain formulation in $P_{zmv}(\Omega)$ for solving numerically (6) is discretized with conforming scalar Lagrange finite elements, and the error is in $L^2(\Omega)$-norm. While if ones uses a mixed approach, discretized with Raviart-Thomas finite elements, the error is in $H(\text{div}; \Omega)$-norm.

To study the variational formulations for the vector potential part, we used a double orthogonality property for the formulation in $H_0(\text{curl}; \Omega)$, respectively a triple orthogonality property for the formulation in $X_N^{\alpha}(\Omega)$, where we considered the latter in a non-convex polyhedron $\Omega$. If $\sigma_{Neu}(\mu) \le \sigma_{Dir}$, we recover the same orders of convergence. For instance, this is automatically true in a 2D polygonal domain for smooth $\mu$, since $\sigma_{Neu}(\mu) = \sigma_{Neu}(1)$ in that case, and it always holds that $\sigma_{Neu}(1) = \sigma_{Dir}$ in 2D. On the other hand, if $\sigma_{Neu}(\mu) > \sigma_{Dir}$, the expected order of convergence is better for the edge finite element discretization than for the vector Lagrange finite elements. This can be explained by the fact that, with the WRM, the norm is sharper as one also controls the divergence of the vector potential $A$, which is of no relevance to compute an approximation of the magnetic field. Finally, in a convex polyhedron, adapting the results of section 8.2.2.B in [6], we observe that, in principle, the same order of convergence is recovered for both methods, namely almost like $O(h^{\sigma_{Neu}(\mu)})$.
**Dedication**

This article is dedicated to Professor Alain Bossavit on the occasion of his
80th birthday.

**Conflict of Interest**

The authors declare that they have no conflict of interest.

**References**

[1] Alonso Rodriguez, A., Bertolazzi, E., Valli, A.: The curl–div system: theory and finite element approximation. Maxwell’s equations-analysis and numerics. Radon Ser. Comput. Appl. Math., vol. 24, pp. 1–43. De Gruyter, Berlin, Germany (2019)

[2] Ciarlet Jr, P.: T-coercivity: Application to the discretization of Helmholtz-like problems. Computers Math. Applic. **64**, 22–34 (2012)

[3] Chesnel, L., Ciarlet Jr, P.: T-coercivity and continuous Galerkin methods: application to transmission problems with sign changing coefficients. Numer. Math. **124**, 1–29 (2013)

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

[5] Monk, P.: Finite Element Methods for Maxwell’s Equations. Oxford University Press, Oxford, Great Britain (2003)

[6] Costabel, M., Dauge, M.: Weighted Regularization of Maxwell’s Equations in polyhedral domains. SIAM J. Numer. Anal. **93**, 239–277 (2002)

[7] Kirsch, A., Hettlich, F.: Mathematical Theory of Time-harmonic Maxwell’s Equations. Springer, Cham, Switzerland (2016)

[8] Assous, F., Ciarlet Jr, P., Labrunie, S.: Mathematical Foundations of Computational Electromagnetism. Springer, Cham, Switzerland (2018)

[9] Gross, P.W., Kotiuga, P.R.: Electromagnetic Theory and Computation: a Topological Approach. MSRI Publications Series, **48**. Cambridge University Press, Cambridge (2004)

[10] Bossavit, A.: Computational Electromagnetism. Academic Press, New York, USA (1998)

[11] Bermúdez, A., Gómez, D., Salgado, P.: Mathematical Models and Numerical Simulation in Electromagnetism. Springer, Cham, Switzerland (2014)
[12] Bermúdez, A., Rodriguez, R., Salgado, P.: A finite element method for the magnetostatic problem in terms of scalar potentials. SIAM J. Numer. Anal. **46**, 1338–1363 (2008)

[13] Bermúdez, A., López-Rodriguez, B., Rodriguez, R., Salgado, P.: Equivalence between two finite element methods for the eddy current problem. C. R. Acad. Sci. Paris, Ser. I **348**, 769–774 (2010)

[14] Tang, Z., Le Menach, Y., Creusé, E., Nicaise, S., Piriou, F., Nemitz, N.: Residual and equilibrated error estimators for magnetostatic problems solved by Finite Element Method. IEEE Trans. Magn. **49**, 5715–5723 (2013)

[15] Bermúdez, A., Piñiero, M., Rodriguez, R., Salgado, P.: Analysis of an ungauged $T, \phi$–$\phi$ formulation of the eddy current problem with currents and voltage excitations. Math. Mod. Num. Anal. **51**, 2487–2509 (2017)

[16] Bermúdez, A., Piñiero, M., Salgado, P.: Mathematical and numerical analysis of a transient magnetic model with voltage drop excitations. Computers Math. Applic. **76**, 2710–2722 (2018)

[17] Creusé, E., Dular, P., Nicaise, S.: About the gauge conditions arising in Finite Element magnetostatic problems **77**, 1563–1582 (2019)

[18] Chicaud, D., Ciarlet Jr, P., Modave, A.: Analysis of variational formulations and low-regularity solutions for time-harmonic electromagnetic problems in complex anisotropic media. SIAM J. Math. Anal. **53**, 2691–2717 (2021)

[19] Chicaud, D., Ciarlet Jr, P.: Analysis of time-harmonic Maxwell impedance problems in anisotropic media. SIAM J. Math. Anal. **55**, 1969–2000 (2023)

[20] Feistauer, M., Ženišek, A.: Finite element solution of nonlinear elliptic problems. Numer. Math. **50**, 451–475 (1987)

[21] Dello Russo, A., Alonso, A.: Finite element approximation of Maxwell eigenproblems on curved Lipschitz polyhedral domains. Appl. Numer. Math. **59**, 1796–1822 (2009)

[22] Bramble, J.H., Pasciak, J.E.: A new approximation technique for div-curl systems. Math. Comp. **73**, 1739–1762 (2003)

[23] Duan, H., Jia, F., Lin, P., Tan, R.C.E.: The local $L^2$ projected $C^0$ finite element method for Maxwell problem. SIAM J. Numer. Anal. **47**, 1274–1303 (2009)
[24] Bonito, A., Guermond, J.-L.: Approximation of the eigenvalue problem for the time harmonic Maxwell system by continuous Lagrange finite elements. Math. Comp. **80**, 1887–1910 (2011)

[25] Badia, S., Codina, R.: A nodal-based finite element approximation of the Maxwell problem suitable for singular solutions. SIAM J. Numer. Anal. **50**, 398–417 (2012)

[26] Duan, H., Du, Z., Liu, W., Zhang, S.: New mixed element for Maxwell equations. SIAM J. Numer. Anal. **57**, 320–354 (2019)

[27] Ern, A., Guermond, J.-L.: Finite Elements II. Springer, Cham, Switzerland (2021)

[28] Jamelot, E., Ciarlet Jr, P.: Fast non-overlapping Schwarz domain decomposition methods for solving the neutron diffusion equation. J. Comput. Phys. **241**, 445–463 (2013)

[29] Ciarlet Jr, P., Jamelot, E., Kpadonou, F.D.: Domain decomposition methods for the diffusion equation with low-regularity solution. Computers Math. Applic. **74**, 2369–2384 (2017)

[30] Giret, L.: Non-Conforming Domain Decomposition for the Multigroup Neutron SPN Equation. PhD thesis, Université Paris-Saclay (2018)

[31] Barré, M., Ciarlet Jr, P.: The T-coercivity approach for mixed problems (2022). https://hal.inria.fr/hal-03820910v1

[32] Jamelot, E.: Improved stability estimates for solving Stokes problem with Fortin-Soulie finite elements. Technical Report HAL (2023)

[33] Ciarlet Jr, P.: Lecture notes on Maxwell's equations and their approximation (in French). Master's degree Analysis, Modelling and Simulation from Paris-Saclay University and Institut Polytechnique de Paris (2022). https://hal.inria.fr/hal-03153780v3

[34] Weber, C.: A local compactness theorem for Maxwell's equations. Math. Meth. Appl. Sci. **2**, 12–25 (1980)

[35] Reitzinger, S., Schöberl, J.: An algebraic multigrid method for finite element discretizations with edge elements. Numer. Linear Algebra Appl. **9**, 223–238 (2002)

[36] Duan, H., Li, S., Tan, R.C.E., Zheng, W.: A delta-regularization finite element method for a double curl problem with divergence-free constraint. SIAM J. Numer. Anal. **50**, 3208–3230 (2012)

[37] Ciarlet, Jr, P., Wu, H., Zou, J.: Edge element methods for Maxwell's
equations with strong convergence for Gauss' laws. SIAM J. Numer. Anal. **52**, 779-807 (2014)

[38] Bermúdez, A., Rodriguez, R., Salgado, P.: Numerical treatment of realistic boundary conditions for the eddy current problem in an electrode via Lagrange multipliers. Math. Comp. **74**, 123–151 (2005)

[39] Raviart, P.-A., Thomas, J.-M.: A mixed finite element method for second order elliptic problems. In: Mathematical Aspects of Finite Element Methods. Lecture Notes in Mathematics, vol. 606, pp. 292–315. Springer, Berlin, Germany (1977)

[40] Bermudez, A., Gamallo, P., Nogueiras, M.R., Rodriguez, R.: Approximation of a structural acoustic vibration problem by hexahedral finite elements. IMA J. Numer. Anal. **26**, 391–421 (2006)

[41] Jerison, D.S., Kenig, C.E.: The Neumann problem on Lipschitz domains. Bull. Amer. Math. Soc. **4**, 203–207 (1981)

[42] Grisvard, P.: Elliptic Problems in Nonsmooth Domains. Pitman, Boston, London, Melbourne (1985)

[43] Dauge, M.: Elliptic Boundary Value Problems on Corner Domains. Springer, Berlin, Germany (1988)

[44] Birman, M.S., Solomyak, M.Z.: Maxwell operator in regions with nonsmooth boundaries. Sib. Math. J. **28**, 12–24 (1987)

[45] Hiptmair, R.: Finite elements in computational electromagnetics. Acta Numerica **50**, 237–339 (2002)

[46] Costabel, M., Dauge, M., Nicaise, S.: Singularities of Maxwell interface problems. Modél. Math. Anal. Numér. **33**, 627–649 (1999)

[47] Ciarlet Jr, P.: On the approximation of electromagnetic fields by edge finite elements. Part 1: sharp interpolation results for low-regularity fields. Computers Math. Applic. **71**, 85–104 (2016)

[48] Costabel, M., Dauge, M., Nicaise, S.: Corner singularities and analytic regularity for linear elliptic systems. Part I: Smooth domains (2010). https://hal.archives-ouvertes.fr/hal-00453934v2

[49] Bonito, A., Guermond, J.-L., Luddens, F.: Regularity of the Maxwell equations in heterogeneous media and Lipschitz domains. J. Math. Anal. Appl. **408**, 498–512 (2013)

[50] Ciarlet Jr, P.: Augmented formulations for solving Maxwell equations. Comp. Meth. Appl. Mech. and Eng. **194**, 559–586 (2005)
30 *Solving magnetostatic systems*

[51] Jamelot, E.: Continuous Galerkin methods for solving Maxwell equations (in French). PhD thesis, Ecole Polytechnique (2005)

[52] Ciarlet Jr, P., Jamelot, E.: Continuous Galerkin methods for solving the time-dependent Maxwell equations in 3d geometries. J. Comput. Phys. **226**, 1122–1135 (2007)

[53] Buffa, A., Ciarlet Jr, P., Jamelot, E.: Solving electromagnetic eigenvalue problems in polyhedral domains. Numer. Math. **113**, 497–518 (2009)

[54] Ciarlet Jr, P., Labrunie, S.: Numerical analysis of the generalized Maxwell equations (with an elliptic correction) for charged particle simulations. Math. Models Meth. App. Sci. **19**, 1959–1994 (2009)

[55] Ciarlet, Jr, P., Scheid, C.: Electrowetting of a 3D drop: numerical modelling with electrostatic vector fields. Math. Mod. Num. Anal. **44**, 647–670 (2010)

[56] Otin, R.: Regularized Maxwell equations and nodal finite elements for electromagnetic field computations. Electromagnetics **30**, 190–204 (2010)

[57] Costabel, M.: A coercive bilinear form for Maxwell's equations. J. Math. Anal. and Appl. **157**, 527–541 (1991)

[58] Costabel, M., Dauge, M.: Singularities of electromagnetic fields in polyhedral domains. Arch. Rational Mech. Anal. **151**, 221–276 (2000)

[59] Clough, R.W., Tocher, J.L.: Finite element stiffness matrices for analysis of plates in bending. In: Conference on Matrix Methods in Structural Mechanics, Wright Patterson A.F.B, Dayton, USA, pp. 515–545 (1965)

[60] Ciarlet, P.G.: Basic error estimates for elliptic problems. Handbook of Numerical Analysis, vol. II, pp. 17–351. North Holland, Amsterdam, The Netherlands (1991)

[61] Percell, P.: On cubic and quartic Clough-Tocher finite elements. SIAM J. Numer. Anal. **13**, 100–103 (1976)

[62] Alfeld, P.: A trivariate Clough-Tocher scheme for tetrahedral data. Comput. Aided Geom. Design **1**, 169–181 (1984)

[63] Sorokina, T., Worsey, A.J.: A multivariate Powell-Sabin interpolant. Adv Comput Math **29**, 71–89 (2008)

[64] Badia, S., Codina, R.: A combined nodal continuous-discontinuous finite element formulation for the Maxwell problem. Applied Mathematics and Computation **218**, 4276–4294 (2011)
[65] Eddargani, S.: Approximation by splines functions on Powell-Sabin tri-angulations. PhD thesis, Universidad de Granada and Université Hassan 1er (2021)

[66] Otin, R.: Regularized Maxwell equations and nodal finite elements for electromagnetic field computations in frequency domain. PhD thesis, Polytechnic University of Catalonia (2011)

[67] Otin, R.: ERMES: A nodal-based finite element code for electromagnetic simulations in frequency domain. Computer Physics Communications 184, 2588-2595 (2013)