DeepAI

# Constraint-preserving hybrid finite element methods for Maxwell's equations

Maxwell's equations describe the evolution of electromagnetic fields, together with constraints on the divergence of the magnetic and electric flux densities. These constraints correspond to fundamental physical laws: the nonexistence of magnetic monopoles and the conservation of charge, respectively. However, one or both of these constraints may be violated when one applies a finite element method to discretize in space. This is a well-known and longstanding problem in computational electromagnetics. We use domain decomposition to construct a family of primal hybrid finite element methods for Maxwell's equations, where the Lagrange multipliers are shown to correspond to a numerical trace of the magnetic field and a numerical flux of the electric flux density. Expressing the charge-conservation constraint in terms of this numerical flux, we show that both constraints are strongly preserved. As a special case, these methods include a hybridized version of Nédélec's method, implying that it preserves the constraints more strongly than previously recognized. These constraint-preserving properties are illustrated using numerical experiments in both the time domain and frequency domain. Additionally, we observe a superconvergence phenomenon, where hybrid post-processing yields an improved estimate of the magnetic field.

• 6 publications
• 6 publications
10/05/2022

### Structure preserving transport stabilized compatible finite element methods for magnetohydrodynamics

We present compatible finite element space discretizations for the ideal...
05/03/2021

### A structure-preserving finite element method for compressible ideal and resistive MHD

We construct a structure-preserving finite element method and time-stepp...
09/17/2020

### Broadband Finite-Element Impedance Computation for Parasitic Extraction

Parasitic extraction is a powerful tool in the design process of electro...
03/23/2020

### Charge-conserving hybrid methods for the Yang-Mills equations

The Yang-Mills equations generalize Maxwell's equations to nonabelian ga...
12/08/2020

### A Finite Element Method for MHD that Preserves Energy, Cross-Helicity, Magnetic Helicity, Incompressibility, and div B = 0

We construct a structure-preserving finite element method and time-stepp...
11/23/2021

### A Conservative Finite Element Solver for MHD Kinematics equations: Vector Potential method and Constraint Preconditioning

A new conservative finite element solver for the three-dimensional stead...
11/16/2022

### Stimulation of soy seeds using environmentally friendly magnetic and electric fields

The study analyzes the impact of constant and alternating magnetic field...

## 1. Introduction

Maxwell’s equations consist of two vector evolution equations, together with two scalar constraint equations,

and , where is magnetic flux density, is electric flux density, and is charge density. These constraints are automatically preserved by the evolution, so given initial conditions satisfying the constraints, one can simply evolve forward in time without needing to “enforce” the constraints in any way.

However, if one applies a finite element method in space, then the semidiscretized evolution equations no longer necessarily preserve these constraints, at least not strongly. Nédélec [29] showed that, if one uses curl-conforming edge elements for the electric field and divergence-conforming face elements for , then the semidiscretized equations preserve strongly. On the other hand, holds only in the Galerkin sense (i.e., when both sides are integrated against certain continuous, piecewise-polynomial test functions). Observing this, Christiansen and Winther [15] commented that strong preservation of both divergence constraints “appears to be necessary for many applications in electromagnetics,” and Houston et al. [20] call this “one of the main difficulties in the numerical solution of Maxwell’s equations.” For this reason, alternative approaches have been developed that enforce the constraints strongly—for instance, using Lagrange multipliers [4, 13]—instead of attempting to preserve them automatically but weakly, as Nédélec’s method does. In cases where , another idea is to use divergence-free elements to construct nonconforming methods [8, 10] or discontinuous Galerkin methods [16, 20, 9].

In this paper, we attack the problem of constraint preservation from a different perspective. We perform domain decomposition of the Lagrangian (i.e., primal) variational principle for Maxwell’s equations, in terms of the vector potential and scalar potential , using Lagrange multipliers and to enforce inter-element continuity and boundary conditions. These Lagrange multipliers are shown to correspond to boundary traces of the magnetic field and electric flux density . After using gauge symmetry to fix , we show that the evolution of automatically preserves the constraints and . Finally, we semidiscretize this domain-decomposed variational principle, obtaining primal hybrid finite element methods that preserve this formulation of the constraints in a strong sense. As a special case, we give a hybridized formulation of Nédélec’s method, implying that it preserves the constraints in a stronger sense than previously recognized.

The paper is organized as follows:

• In Section 2, we review Maxwell’s equations, the Lagrangian variational principle, and semidiscretization using edge elements.

• In Section 3, we domain decompose the Lagrangian variational principle, relate solutions to the classical (non-domain-decomposed) formulation of Maxwell’s equations, and study the domain-decomposed version of the constraints and their preservation.

• In Section 4, we consider primal hybrid finite element methods for semidiscretizing the domain-decomposed evolution equations, showing that constraints are preserved in a strong sense.

• Finally, in Section 5 we conduct numerical experiments demonstrating the behavior of the hybridized Nédélec method. In addition to the constraints being preserved to machine precision, these results illustrate a superconvergence phenomenon for the post-processed magnetic field , similar to that observed for other hybridized mixed methods (cf. Arnold and Brezzi [2], Brezzi et al. [11]).

## 2. Maxwell’s equations

### 2.1. Maxwell’s equations

We begin by reviewing the classical formulation of Maxwell’s equations, first in terms of the electric and magnetic fields and flux densities, and then in terms of the vector and scalar potentials. We postpone the discussion of regularity until the introduction of the weak formulation, in Section 2.2

; for the moment, everything may be assumed to be smooth.

#### 2.1.1. Standard formulation

In their most familiar form, Maxwell’s equations consist of the vector evolution equations,

 (1a) ˙B =−curlE, (1b) ˙D+J =curlH,

together with the scalar constraint equations,

 (2a) divB =0, (2b) divD =ρ.

Here, and denote the electric field and magnetic field, and denote the electric flux density and magnetic flux density, and

are the electric permittivity and magnetic permeability tensors, and

and are current density and charge density, respectively. We use the “dot” notation to denote partial differentiation with respect to time.

The evolution equations (1) automatically preserve the constraints (2). Indeed, taking the divergence of (1a) implies , so (2a) is preserved. Similarly, taking the divergence of (1b) implies , so (2b) is preserved if and only if and satisfy , which is the law of conservation of charge. We refer to (2b) as the charge-conservation constraint, since it is equivalent to this condition.

#### 2.1.2. Formulation in terms of potentials

Alternatively, Maxwell’s equations may be expressed in terms of a vector field , called the vector potential, and a scalar field , called the scalar potential. Given and , we define the electric field and magnetic flux density by

Note that (1a) and (2a) are automatically satisfied, so we may restrict our attention entirely to the single evolution equation (1b), which we have already seen preserves (2b).

However, Maxwell’s equations do not uniquely determine the evolution of . Observe that if is any time-dependent scalar field, then the transformation leaves , , , unchanged. Such transformations are called gauge transformations, and the invariance of Maxwell’s equations under gauge transformations is called gauge symmetry. In particular, any solution may be transformed into one of the form by taking to be a solution of . Therefore, we may restrict our attention to solutions with .

###### Remark .

This procedure of restricting to particular solutions, which are related to a general solution by some gauge transformation, is called gauge fixing. The choice , called temporal gauge, is the most convenient for our purposes, but there are other choices as well. Note that there is still some remaining gauge symmetry, even after performing temporal gauge fixing: we may transform for any constant in time.

After temporal gauge fixing, we can write (1b) as either a first-order system in , ,

 ˙A=−ϵ−1D,˙D+J=curl(μ−1curlA),

or as a second-order equation in alone,

 −∂t(ϵ˙A)+J=curl(μ−1curlA).

In the special case where and are simply positive constants with (as in vacuum, with units chosen so that the speed of light is ) and , the latter equation just becomes

 ¨A+curlcurlA=0.

Taking the Fourier transform with respect to time (the so-called

frequency domain or time-harmonic

approach), this latter equation transforms into the eigenvalue problem for the

operator.

### 2.2. Weak formulation

We next discuss the weak formulation of Maxwell’s equations, first using a Lagrangian variational principle in terms of the potentials and , and then fixing the temporal gauge to arrive at a weak formulation in terms of alone.

#### 2.2.1. Function spaces and regularity

Let be a bounded Lipschitz domain, and define the function spaces

 H1(Ω) \coloneqq{u∈L2(Ω):gradu∈L2(Ω,R3)}, H(curl;Ω) \coloneqq{u∈L2(Ω,R3):curlu∈L2(Ω,R3)}, H(div;Ω) \coloneqq{u∈L2(Ω;R3):divu∈L2(Ω)}.

We also define the following subspaces, with boundary conditions imposed:

 ˚H1(Ω) \coloneqq{u∈H1(Ω):u|∂Ω=0}, ˚H(curl;Ω) \coloneqq{u∈H(curl;Ω):u×n|∂Ω=0}, ˚H(div;Ω)

Here, denotes the outer unit normal to , and restrictions to are interpreted in the trace sense.

Let be a curve in and be a curve in . It follows that is a curve in , that is a curve in , and that (1a) and (2a) hold strongly in . We also assume that both and are , symmetric, and uniformly elliptic. In particular, this implies that and are both curves in . Henceforth, we restrict our attention to such that is in fact a curve in .

Finally, let the current density be a given curve in and the charge density be a given curve in , satisfying the charge conservation condition .

#### 2.2.2. The Lagrangian and Euler–Lagrange equations

For as above, define the Lagrangian

 L(A,φ,˙A,˙φ)\coloneqq∫Ω(12E⋅D−12B⋅H+A⋅J−φρ).

The Euler–Lagrange equations are

 (3a) ∫Ω(A′⋅(˙D+J)−curlA′⋅H) =0, ∀A′ ∈˚H(curl;Ω), (3b) ∫Ω(gradφ′⋅D+φ′ρ) =0, ∀φ′ ∈˚H1(Ω),

which are weak expressions of (1b) and (2b), respectively.

These Euler–Lagrange equations imply that solutions have additional regularity properties. Since is in , we have that is in . Likewise, since is in , we have that is in . Hence, solutions to this weak problem are in fact strong solutions of Maxwell’s equations.

###### Remark .

When and are constant in time, the electric and magnetic fields have precisely the same regularity assumed by Monk [27, eqs. (7)–(8)], namely: is in and in , while is in and in .

As in Section 2.1, this formulation is symmetric with respect to gauge transformations , where is now an arbitrary curve in . Fixing the temporal gauge , the Lagrangian becomes

 L(A,˙A)=∫Ω(12E⋅D−12B⋅H+A⋅J),

and the Euler–Lagrange equations are just (3a). This again implies that is in , so (1b) holds strongly. By the same argument as in Section 2.1, this automatically preserves the charge-conservation constraint.

###### Remark .

Preservation of the charge-conservation constraint may also be seen as a consequence of the remaining gauge symmetry , mentioned in Section 2.1.2, where is constant in time. This is a particular instance of Noether’s theorem, which relates symmetries to conservation laws. See Marsden and Ratiu [24, Section 1.6] for an account of the case, as well as the discussion in Christiansen and Winther [15].

### 2.3. Galerkin semidiscretization using Nédélec elements

The use of finite elements in computational electromagnetics is a broad topic with a long history, and we do not attempt to give a full account here. We refer the reader to the texts by Monk [28] and Jin [21], as well as the excellent survey article by Hiptmair [19], which relates these methods to the more recent theory of finite element spaces of differential forms. In this section, we briefly review the semidiscretization of Maxwell’s equations using the elements of Nédélec [29, 30], an approach that was subsequently analyzed in a series of papers by Monk [25, 26, 27].

Galerkin semidiscretization of the variational problem (3a) restricts the trial and test functions to some finite-dimensional subspace , resulting in a finite-dimensional system of ODEs. That is, we seek a curve such that

 (4) ∫Ω(A′h⋅(˙Dh+J)−curlA′h⋅Hh)=0,∀A′h∈V1h,

where , , , and . The discrete versions of (1a) and (2a),

 ˙Bh=−curlEh,divBh=0,

follow immediately. In fact, both hold strongly in , by the same argument as in Section 2.2.1, since and . On the other hand, we cannot conclude that is in , nor that is in , since (4) only holds for test functions in and not all of .

Consequently, the charge-conservation constraint (2b) is only preserved in the following, much weaker sense. Let be a finite-dimensional subspace such that . Then, for all , taking in (4) and applying gives

Hence, if the initial conditions satisfy , for all , then this is preserved by the flow of (4).

In particular, suppose now that is polyhedral, and that is a triangulation of by -simplices (i.e., tetrahedra) . We may take to be the space of continuous degree- piecewise polynomials on vanishing on , corresponding to standard Lagrange finite elements. For , we may take either degree- Nédélec edge elements of the first kind [29] or degree- Nédélec edge elements of the second kind [30]

with vanishing degrees of freedom on

. These are spaces of piecewise-polynomial vector fields in with tangential (but not necessarily normal) continuity between neighboring simplices. These choices ensure that , so the weak charge-conservation argument above holds.

Note, however, that only says that holds in an “averaged” sense, since (unlike in the infinite-dimensional case) nonzero cannot be taken to have arbitrarily small support. We cannot even conclude that the constraint holds in the sense that , since the indicator function is discontinuous and therefore not an admissible test function. (Christiansen and Winther [15] give a compactness argument for why this weak form of the constraint “might be just as good” as the strong form, in the limit as ; see also Christiansen [14].) This motivates our proposed hybrid approach, based on domain decomposition, for which piecewise-constants are admissible test functions.

###### Remark .

The method above describes the evolution of . Equivalently, one may evolve and by augmenting (4) with . This is the original approach described by Nédélec [29], where is given by face elements on .

## 3. Domain decomposition

In this section, we introduce an alternative variational formulation for Maxwell’s equations, based on domain decomposition. Specifically, we decompose the problem on into a collection of problems on , weakly enforcing internal continuity and external boundary conditions using Lagrange multipliers. This is similar in spirit to the standard approach to domain decomposition for Poisson’s equation, cf. Brezzi and Fortin [12]. We show that the Lagrange multipliers enforcing these conditions on and correspond to the traces of and , respectively, and we show that the latter satisfies an appropriate version of the charge-conservation constraint.

### 3.1. Function spaces

We begin by introducing the following discontinuous function spaces, which are larger than the spaces used in the previous variational formulation:

 H1(Th) \coloneqq{u∈L2(Ω):u|K∈H1(K), for all K∈Th}, H(curl;Th) \coloneqq{u∈L2(Ω,R3):u|K∈H(curl;K), for all K∈Th}, H(div;Th) \coloneqq{u∈L2(Ω,R3):u|K∈H(div;K), for all K∈Th}.

Brezzi and Fortin [12, Proposition III.1.1] show that

 ˚H1(Ω)={u∈H1(Th):∑K∈Th∫∂Kuλ⋅n=0, for all % λ∈H(div;Ω)}.

That is, is the subspace of where internal continuity and external boundary conditions are enforced by Lagrange multipliers . Likewise, [12, Proposition III.1.2] shows that

 H(div;Ω)={u∈H(div;Th):∑K∈Th∫∂Kuλ⋅n=0, for all λ∈˚H1(Ω)}.

Using a similar argument, we now prove the corresponding result for the spaces.

.

###### Proof.

If , then for any , we have

 ∑K∈Th∫∂K(u×λ)⋅n =∑K∈Th∫K(curlu⋅λ−u⋅curlλ) =∫Ω(curlu⋅λ−u⋅curlλ) =∫∂Ω(u×λ)⋅n =0,

so the forward inclusion holds. To get the reverse inclusion , suppose that satisfies the condition above, and let . Then, integrating by parts, we have

 =∣∣∣∑K∈Th∫Kcurlu⋅λ−∑K∈Th∫∂K(u×λ)⋅n∣∣∣ =∣∣∣∑K∈Th∫Kcurlu⋅λ∣∣∣ ≤(∑K∈Th∥curlu∥2L2(K,R3))1/2∥λ∥L2(Ω,R3),

where the last line uses the triangle and Cauchy–Schwarz inequalities. It follows that , so . This implies that for all . Hence, in the trace sense, which completes the proof. ∎

###### Remark .

A variant of this result is stated in Boffi et al. [7, Proposition 2.1.3], where is taken to be in rather than . However, the version given here is more natural for the purposes of the hybrid methods discussed in Section 4.

### 3.2. Domain decomposition of the Lagrangian variational principle

We now introduce a new Lagrangian for Maxwell’s equations, which allows the potentials to live in the discontinuous function spaces defined in the previous section, enforcing continuity and boundary conditions using Lagrange multipliers.

Let and , and introduce the Lagrange multipliers and . We adopt the notation, often seen in the literature on discontinuous Galerkin and hybrid methods, of placing hats over variables that act like weak traces/fluxes. As before, suppose that is and that is , such that is . Furthermore, suppose that and are both . Define the Lagrangian

The Euler–Lagrange equations are then

 (5a) ∫K(A′⋅(˙D+J)−curlA′⋅H)+∫∂K(A′×ˆH)⋅n =0, ∀A′ ∈H(curl;K), (5b) ∫K(gradφ′⋅D+φ′ρ)−∫∂Kφ′ˆD⋅n =0, ∀φ′ ∈H1(K), (5c) ∑K∈Th∫∂K(A×ˆH′)⋅n =0, ∀ˆH′ ∈H(curl;Ω), (5d) ∑K∈Th∫∂KφˆD′⋅n =0, ∀ˆD′ ∈H(div;Ω),

where (5a) and (5b) hold for all . We now relate this to the classical variational form of Maxwell’s equations, stated in (3).

###### Proposition .

is a solution to (5) if and only if is a solution to (3) with and . In particular, if is a solution to (3), then is a solution to (5).

###### Proof.

Suppose is a solution to (5). By Section 3.1, (5c) implies , so taking and summing (5a) over , the integrals over cancel, yielding (3a). As previously stated, (3a) implies , so substituting this into (5a) gives

 ∫∂K(A′×ˆH)⋅n=∫K(curlA′⋅H−A′⋅curlH)=∫∂K(A′×H)⋅n,∀A′∈H(curl;K),

so . Similarly, (5d) implies , so taking and summing (5b) over yields (3b). This implies , and substituting into (5b) gives .

Conversely, suppose is a solution to (3). Since and , it follows that (5c) and (5d) hold. Furthermore, (3) implies that and , so (5a) and (5b) hold with and . In particular, we could take and . ∎

###### Remark .

Note that, in addition to (5b) implying that , we also see by taking that satisfies the conservation law , for all .

### 3.3. Temporal gauge fixing and the charge-conservation constraint

As in Section 2.2, if is a solution to (5), then so is for any curve . Therefore, we perform temporal gauge fixing by taking . This yields the gauge-fixed Lagrangian

 L(A,ˆH,˙A,˙ˆH)=∑K∈Th[∫K(12E⋅D−12B⋅H+A⋅J)+∫∂K(A×ˆH)⋅n],

whose Euler–Lagrange equations are simply (5a) and (5c). Of course, (5d) is satisfied trivially, since . The next result shows that the charge-conservation constraint (5b) is automatically preserved, for an appropriately-defined .

###### Proposition .

Let be a solution to (5a) and (5c). Suppose initial values for satisfy (5b), and let be the solution to . Then is a solution to (5).

###### Proof.

As we have already mentioned, trivially satisfies (5d), so it suffices to show that (5b) holds. Let be arbitrary. Taking in (5a) and integrating by parts gives

so if (5b) holds at the initial time, then it holds for all time. ∎

###### Remark .

As in Section 3.2, taking implies . Furthermore, if the initial conditions also satisfy , then we have for all time, since . Finally, if , and if the initial conditions for equal those for , then we recover .

Finally, we express this variational problem in the standard notation used for mixed and hybrid finite element methods, in terms of a pair of bilinear forms [12, Chapter II]. We will make use of this notation throughout the subsequent sections. Defining

 a :H(curl;Th)×H(curl;Th)→R, a(A,A′) \coloneqq∑K∈Th∫KcurlA′⋅μ−1curlA, b :H(curl;Th)×H(curl;Ω)→R, b(A′,ˆH) \coloneqq−∑K∈Th∫∂K(A′×ˆH)⋅n,

we seek and such that

 (6a) ⟨˙D+J,A′⟩ =a(A,A′)+b(A′,ˆH), ∀A′ ∈H(curl;Th), (6b) 0 =b(A,ˆH′), ∀ˆH′ ∈H(curl;Ω),

where is the inner product. Defining the map , , we see that (6) is equivalent to evolving by

 ⟨˙D+J,A′⟩=a(A,A′),∀A′∈kerB,

and subsequently solving for satisfying (6a). Since by Section 3.1, it follows that solves the non-domain-decomposed problem (3a).

## 4. Hybrid semidiscretization

We now perform Galerkin semidiscretization of the domain-decomposed variational problem with temporal gauge fixing, as introduced in the previous section. This results in a hybrid method for Maxwell’s equations, where “hybrid” means that the Lagrange multipliers and their test functions are both restricted to a subspace of . We then show that a suitably-defined satisfies the charge-conservation constraint in a strong sense, as opposed to the much weaker sense in which was seen to satisfy this constraint in Section 2.3. Finally, we discuss how certain choices of elements yield a hybridized version of Nédélec’s method, while others give nonconforming methods, and we remark on how this framework also applies to the 2-D Maxwell equations.

### 4.1. Semidiscretization of the variational problem

For each , let be a finite-dimensional subspace, so , and let . We seek and such that

 (7a) ∫K(A′h⋅(˙Dh+J)−curlA′h⋅Hh)+∫∂K(A′h×ˆHh)⋅n =0, ∀A′h ∈V1h(K), (7b) ∑K∈Th∫∂K(Ah×ˆH′h)⋅n =0, ∀ˆH′h ∈ˆV1h,

where (7a) holds for all . These are the semidiscretized versions of (5a) and (5c).

###### Remark .

Since (7b) only holds for test functions in , but not necessarily an arbitrary test function in , in general a solution will have . Hence, this method is generally not curl-conforming and is distinct from the conforming methods discussed in Section 2.3.

In terms of the bilinear forms and , this method may be written as

 (8a) ⟨˙Dh+J,A′h⟩ =a(Ah,A′h)+b(A′h,ˆHh), ∀A′h ∈V1h, (8b) 0 =b(Ah,ˆH′h), ∀ˆH′h ∈ˆV1h.

Defining the operator , , we see that (8) is equivalent to evolving by

 (9) ⟨˙Dh+J,A′h⟩=a(Ah,A′h),∀A′h∈kerBh,

and subsequently solving for satisfying (8a).

Since is finite-dimensional, we may apply Banach’s closed range theorem to deduce that is in the range of , so a solution exists, although generally not uniquely. A natural choice is to find the solution minimizing , which in a weak sense minimizes the distance between and . This existence-without-uniqueness is typical of hybrid methods, and one may formally resolve this by replacing by the quotient space (cf. Brezzi and Fortin [12, IV.1.3]). In practice, the evolution on specified by (9) is the essence of the method, and solving for may be seen as an optional post-processing step.

### 4.2. Preservation of the charge-conservation constraint

In order to discuss the charge-conservation constraint, we first suppose that are such that and for all . We consider whether the following discretization of (5b) is preserved,

for suitably defined.

###### Theorem 4.1.

Let be a solution to (7). Suppose initial values for satisfy (10), and let be the solution to . Then (10) holds for all time. In particular, . Moreover, if holds at the initial time, then it holds for all time.

###### Proof.

The proof is essentially similar to that of Section 3.3. Given , taking in (7a) and integrating by parts,