An arbitrary-order method for magnetostatics on polyhedral meshes based on a discrete de Rham sequence

Highlights

• Exact global discrete de Rham sequence of arbitrary order.

• First scheme for the mixed formulation on arbitrary polyhedral cells.

• Unconditional stability and well-posedness of the discrete problem.

• Favourable convergence rates of $h^{k+1}$ for both magnetic field and vector potential.

• Fewer unknowns and simpler interpolators than Finite Element methods on hexahedra.

Abstract

In this work we develop a discretisation method for the mixed formulation of the magnetostatic problem supporting arbitrary orders and polyhedral meshes. The method is based on a global discrete de Rham (DDR) sequence, obtained by patching the local spaces constructed in [22] by enforcing the single-valuedness of the components attached to the boundary of each element. The first main contribution of this paper is a proof of exactness relations for this global DDR sequence, obtained leveraging the exactness of the corresponding local sequence and a topological assembly of the mesh valid for domains that do not enclose any void. The second main contribution is the formulation and well-posedness analysis of the method, which includes the proof of uniform Poincaré inequalities for the discrete divergence and curl operators. The convergence rate in the natural energy norm is numerically evaluated on standard and polyhedral meshes. When the DDR sequence of degree $k\ge 0$ is used, the error converges as $h^{k+1}$, with h denoting the mesh size.

Introduction

In this work we develop a discretisation method for the mixed formulation of the magnetostatic problem supporting arbitrary orders and polyhedral meshes. The stability of the method hinges on a global version of the discrete de Rham (DDR) sequence of [22].

Let $\mathrm{\Omega }\subset {\mathbb{R}}^3$ be an open connected polyhedral domain that does not enclose any void (that is, its second Betti number is zero), with boundary ∂Ω and unit outward normal n. Denote by $\mathbf{H}(\mathbf{curl};\mathrm{\Omega })$ the space of vector-valued functions over Ω that are square-integrable along with their curl and by $\mathbf{H}(\mathrm{div};\mathrm{\Omega })$ the space of vector-valued functions over Ω that are square-integrable along with their divergence. Let $\mathit{J}\in \mathbf{curl}\mathbf{H}(\mathbf{curl};\mathrm{\Omega })$ and $\mathit{g}\in {\mathrm{L}}^2{(\partial \mathrm{\Omega })}^3$ denote, respectively, the free current density and boundary datum. We consider the following problem (see, e.g., [1, Section 4.5.3]): Find $(\mathit{H},\mathit{A})\in \mathbf{H}(\mathbf{curl};\mathrm{\Omega })\times \mathbf{H}(\mathrm{div};\mathrm{\Omega })$ such that

$$a(\mathit{H},\mathit{\zeta })-b(\mathit{\zeta },\mathit{A})=-\underset{\partial \mathrm{\Omega }}{\int }\mathit{g}\cdot \mathit{\zeta }\forall \mathit{\zeta }\in \mathbf{H}(\mathbf{curl};\mathrm{\Omega }),b(\mathit{H},\mathit{v})+c(\mathit{A},\mathit{v})=\underset{\mathrm{\Omega }}{\int }\mathit{J}\cdot \mathit{v}\forall \mathit{v}\in \mathbf{H}(\mathrm{div};\mathrm{\Omega }),$$

with bilinear forms $a:\mathbf{H}(\mathbf{curl};\mathrm{\Omega })\times \mathbf{H}(\mathbf{curl};\mathrm{\Omega })\to \mathbb{R}$, $b:\mathbf{H}(\mathbf{curl};\mathrm{\Omega })\times \mathbf{H}(\mathrm{div};\mathrm{\Omega })\to \mathbb{R}$, and $c:\mathbf{H}(\mathrm{div};\mathrm{\Omega })\times \mathbf{H}(\mathrm{div};\mathrm{\Omega })\to \mathbb{R}$ such that, for all $\mathit{\upsilon },\mathit{\zeta }\in \mathbf{H}(\mathbf{curl};\mathrm{\Omega })$ and all $\mathit{w},\mathit{v}\in \mathbf{H}(\mathrm{div};\mathrm{\Omega })$,

$$a(\mathit{\upsilon },\mathit{\zeta }):=\underset{\mathrm{\Omega }}{\int }\mu \mathit{\upsilon }\cdot \mathit{\zeta },b(\mathit{\zeta },\mathit{v}):=\underset{\mathrm{\Omega }}{\int }\mathit{v}\cdot \mathbf{curl}\mathit{\zeta },c(\mathit{w},\mathit{v}):=\underset{\mathrm{\Omega }}{\int }\mathrm{div}\mathit{w}\mathrm{div}\mathit{v}.$$

The solution of problem (1) satisfies almost everywhere

$$\mu \mathit{H}-\mathbf{curl}\mathit{A}=\mathbf{0}\text{in }\mathrm{\Omega },$$

$$\mathbf{curl}\mathit{H}=\mathit{J}\text{in }\mathrm{\Omega },$$

$$\mathrm{div}\mathit{A}=0\text{in }\mathrm{\Omega },$$

$$\mathit{A}\times \mathit{n}=\mathit{g}\text{on }\partial \mathrm{\Omega }.$$

In the context of magnetostatics, the interpretation of the above relations is as follows: equation (2a) expresses the magnetic field H in terms of the vector potential A, with $\mu :\mathrm{\Omega }\to \mathbb{R}$ denoting the permeability such that $\bar{\mu }\ge \mu \ge \underset{\_}{\mu }>0$ almost everywhere in Ω, where $\bar{\mu }$ and $\underset{\_}{\mu }$ are constant numbers; equation (2b) is Ampère's law relating the magnetic field to the current density; (2c) is the so-called Coulomb's gauge which, together with the boundary condition (2d), ensures the uniqueness of the vector potential (notice that, since the second Betti number of Ω is zero, the space of 2-harmonic forms is trivial).

The well-posedness of problem (1) hinges on the fact that, under the above assumptions on the domain, the image of the curl operator coincides with the kernel of the divergence operator, and the latter is surjective in ${\mathrm{L}}^2(\mathrm{\Omega })$; cf., e.g., the discussion in [22, Section 2]. These relations correspond to the exactness of the rightmost portion of the de Rham sequence

where $i_{\mathrm{\Omega }}$ is the operator that maps a real value to a constant function over Ω and ${\mathrm{H}}^1(\mathrm{\Omega })$ the space of scalar-valued functions over Ω that are square-integrable along with their gradient. The design of stable numerical approximations of problem (1) requires to mimic these exactness properties at the discrete level.

Other formulations of the magnetostatic problem are possible, most notably the H- and A-methods presented in [24]. The H-method consists in eliminating A from (2a), (2b), (2c), (2d) by expressing the fact that $\mu \mathit{H}\in \mathbf{H}(\mathbf{curl};\mathrm{\Omega })$ (see (2a)) to write $\mathrm{div}(\mu \mathit{H})=0$, and by retaining (2b); the divergence constraint on $\mu \mathit{H}$ is enforced in the weak formulation by the introduction of a Lagrange multiplier p (scalar potential, which appears through its gradient). The A-method uses (2a) to eliminate H from (2b), leading to the second order $\mathbf{curl}({\mu }^{-1}\mathbf{curl}\mathit{A})=\mathit{J}$ equation on A. Here too, the constraint $\mathrm{div}\mathit{A}=0$ is imposed in the weak formulation via a scalar potential that appears through its gradient. In either of these methods, the well-posedness of the weak formulation hinges on the exactness of the leftmost portion of the de Rham sequence (3), specifically the relation $\mathrm{Ker}\mathbf{curl}\subset \mathrm{Im}\mathbf{grad}$, which holds only when the first Betti number of Ω is zero. The most convenient formulation to use as a starting point for the numerical approximation (magnetic field-vector potential as in (1), or H- or A-methods as in [24]) thus depends on the topology of the domain: for solid toroidal or cylindrical domains (that have non-zero first Betti number but zero second Betti number), e.g., problem (1) is naturally well-posed, whereas the formulations of [24] require an additional condition to account for the existence of non-trivial 1-harmonic forms. From the standpoint of numerical approximation, this necessitates the (possibly expensive) computation of cohomology generators (see, e.g., [23], [27]). Clearly, for domains that enclose voids (non-zero second Betti number) but do not have tunnels (zero first Betti number), the situation is reversed, and an additional condition enforcing the ${\mathrm{L}}^2$-orthogonality of the vector potential A to 2-harmonic forms is required to ensure the well-posedness of (1); see, e.g., [1, Section 4.5.4]. We note in passing that, although this paper focuses on the formulation (1), the framework we develop provides an entire discrete De Rham sequence that could equally be used on the H- and A-methods.

In the context of Finite Element (FE) approximations, the exactness property discussed above is achieved by a nontrivial choice of finite-dimensional subspaces of $\mathbf{H}(\mathbf{curl};\mathrm{\Omega })$ and $\mathbf{H}(\mathrm{div};\mathrm{\Omega })$; cf. [1] for a comprehensive introduction to this topic. Finite Elements can, however, display severe practical limitations: the construction of the finite-dimensional spaces hinges upon conforming meshes with elements of simple geometric shape; the number of degrees of freedom on hexahedral elements can become very large when increasing the polynomial degree (see, e.g., [22, Table 2]); the definition of unisolvent degrees of freedom can be tricky for high-order versions (cf., e.g., [6] and references therein). A recent generalisation of FE methods is provided by the Isogeometric Analysis (IGA), which is designed to facilitate exchanges with Computer Assisted Design software. In this framework, spline spaces and projection operators that verify a de Rham diagram have been developed in [10], placing on solid grounds the IGA method for electromagnetism originally proposed in [11]; see also [12] for further developments. Low-order frameworks that involve exact discrete sequences of spaces and operators on general polyhedral meshes have been developed over the last years, among which we cite Mimetic Finite Differences (see [5] and references therein), the Discrete Geometric Approach (see, e.g., [17]), or Compatible Discrete Operators (see [8], [9] and also [7]). The above polyhedral frameworks are tightly related to the lowest-order version of the DDR sequence corresponding to $k=0$, and have already found their way in commercial codes. Compatible Discrete Operators, along with the more recent Hybrid High-Order (HHO) technology [20], are available in widely used simulators such as Code_Aster (https://www.code-aster.org) or Code_Saturne (https://www.code-saturne.org). More recently, a de Rham sequence of arbitrary-order virtual spaces has been proposed in [4]; see also the related works [2], [3] concerning the Virtual Element approximation of the magnetostatic problem based on the formulation of [24]. Virtual spaces are spanned by functions whose expression is not available at each point, therefore projections on polynomial spaces are used in practice. For this reason, the exactness of the virtual sequence cannot be directly exploited to prove the stability of numerical schemes. Fully nonconforming approximations of the magnetostatic problem have also been explored, where stability is ensured by additional penalty terms. We mention here, in particular, the Discontinuous Galerkin method of [26], the Hybridizable Discontinuous Galerkin methods of [15], [16], [25], and, on polyhedral meshes, the HHO methods of [13], [14].

The approach proposed in this work relies on a global DDR sequence involving spaces of (polynomial) discrete unknowns and discrete counterparts of vector operators acting thereon. This sequence is obtained by patching the local spaces constructed in [22] by enforcing the single-valuedness of the components of these spaces on the element boundaries. Its exactness, which constitutes the first key contribution of this work, is proved in Theorem 3 below, leveraging the exactness of the local sequence and a topological assembly of the mesh valid for domains that do not enclose any void. The second contribution of this work is the development and well-posedness analysis of a discretisation method for the mixed formulation (1) of the magnetostatic problem (the first, to our knowledge, supporting polyhedral meshes). The discrete problem is formulated in terms of the spaces and operators appearing in the global DDR sequence, along with discrete counterparts of ${\mathrm{L}}^2$-products. For this reason, the stability and well-posedness (Theorem 10 and Corollary 11) of the discrete problem are direct consequences of the exactness of the DDR sequence, together with uniform Poincaré inequalities for the discrete divergence and curl operators. Besides supporting polyhedral meshes and arbitrary orders, the proposed method has fewer unknowns than (non-serendipity) Finite Elements on hexahedra (cf. Remark 8 below) and allows for great freedom in the practical implementation of the polynomial spaces that lie at its core. The convergence rate of the method is numerically evaluated on a set of standard and polyhedral refined mesh families. When the DDR sequence of degree $k\ge 0$ is used, the error in the natural energy norm associated with the problem behaves as $h^{k+1}$, with h denoting the mesh size.

The rest of this paper is organised as follows. In Section 2 we introduce the setting, recalling the appropriate notion of polyhedral mesh, the definitions of vector operators on faces, and those of the local polynomial spaces. In Section 3 we define the global DDR sequence and prove the required exactness relations. Section 4 contains the statement of the discrete problem along with a theoretical well-posedness result and a numerical assessment of its convergence rate. In Section 5 we discuss the practical implementation. Finally, Appendix A contains the proof of the stability result, which is based on uniform discrete Poincaré inequalities. The paper is structured so as to offer two levels of reading. In particular, the implementation aspects in Section 5 and the proofs in Appendix A are rather technical, and the reader mainly interested in the formulation of the proposed numerical scheme can skip them at first reading.