A mass-lumped mixed finite element method for Maxwell's equations
AA mass-lumped mixed finite element method forMaxwell’s equations
Herbert Egger and Bogdan Radu
Abstract
A novel mass-lumping strategy for a mixed finite element approximationof Maxwell’s equations is proposed. On structured orthogonal grids the resultingmethod coincides with the spatial discretization of the Yee scheme. The proposedmethod, however, generalizes naturally to unstructured grids and anisotropic mate-rials and thus yields a variational extension of the Yee scheme for these situations.
We consider the propagation of electromagnetic radiation through a linear non-dispersive and non-conducting medium described by Maxwells equations ε∂ t E = curl H , (1) µ∂ t H = − curl E . (2)Here E , H denote the electric and magnetic field intensities and ε , µ are the symmet-ric and positive definite permittivity and permeability tensors. For ease of notation,we assume that E × n = ε ∂ t e = C (cid:48) h , (3)M µ ∂ t h = − C e . (4) H. EggerDept. of Mathematics, TU Darmstadt, Dolivostraße 15, 64293 Darmstadt, Germany e-mail:[email protected]. RaduGraduate School for Computational Engineering, TU Darmstadt, Dolivostraße 15, 64293 Darm-stadt, Germany e-mail: [email protected] 1 a r X i v : . [ m a t h . NA ] O c t Herbert Egger and Bogdan Radu
Due to the particular structure of the system, the stability of such discretizationschemes can easily be ensured by the simple algebraic conditions(i) C (cid:48) = C (cid:62) ,(ii) M ε , M µ are symmetric and positive definite.In order to enable an efficient solution of (3)–(4) by explicit time-stepping methods,one additionally has to assume that(iii) M − ε , M − µ can be applied efficiently.The finite difference approximation of (1)–(2) on staggered orthogonal grids yieldsapproximations of the form (3)–(4) satisfying (i)–(iii) with diagonal M ε , M µ [12].Moreover, the entries e i , h j in the solution vectors yield second order approxima-tions for the line integrals of E , H along edges of the primal and duals grid, re-spectively [3, 10]. An extension to unstructured grids and anisotropic coefficientsis possible [2, 9], but these approaches rely on the use of two sets of grids whichmakes a rigorous convergence analysis difficult.The finite element approximation of (1)–(2) yields systems of the form (3)–(4)satisfying (i)–(ii) automatically and a rigorous convergence analysis is possible inrather gerneral situations [7, 8]. Condition (iii) is, however, not valid in general, al-though the matrices M ε and M µ are usually sparse. This lack of efficiency can beovercome by mass-lumping, which aims at approximating M ε and M µ by diago-nal or block-diagonal matrices; [6, 4]. These approaches are usually based on anenrichment of the approximation spaces and appropriate quadrature; see [3].In this paper, we present a novel mass-lumping strategy for a mixed finite el-ement approximation of (1)–(2) that yields properties (i)–(iii) without such an in-crease of the system dimension. In special cases, the resulting scheme reduced tothe staggered-grid finite difference approximation of the Yee scheme. As a preliminary step, we consider a mass-lumped mixed finite element approxima-tion based on enriched approximation spaces and numerical quadrature. We seek forapproximations (cid:101) E h ( t ) ∈ (cid:101) V h , (cid:101) H h ( t ) ∈ (cid:101) Q h satisfying ( ε∂ t (cid:101) E h ( t ) , (cid:101) v h ) h = ( (cid:101) H h ( t ) , curl (cid:101) v h ) ∀ (cid:101) v h ∈ (cid:101) V h , (5) ( µ∂ t (cid:101) H h ( t ) , (cid:101) q h ) h , ∗ = − ( curl (cid:101) E h ( t ) , (cid:101) q h ) ∀ (cid:101) q h ∈ (cid:101) Q h , (6)for all t >
0. Here, (cid:101) V h ⊂ H ( curl; Ω ) and (cid:101) Q h ⊂ L ( Ω ) are appropriate finite di-mensional subspaces and ( a , b ) h , ( a , b ) h , ∗ are approximations for the scalar product ( a , b ) = (cid:82) Ω a ( x ) · b ( x ) dx to be defined below.We restrict our discussion in the sequel to problems where E = ( E x , E y , ) and H = ( , , H z ) with E x , E y , H z independent of z , which allows to represent the fieldsin two dimensions. The extension to three dimensions will be discussed in Section 5. mass-lumped mixed finite element method for Maxwell’s equations 3 Let T h = { T } be a conforming mesh of Ω consisting of triangles and parallelo-grams. Every element T ∈ T h is the image F T ( (cid:98) T ) of a reference triangle or referencesquare under an affine mapping F T ( (cid:98) x ) = a T + B T (cid:98) x with a T ∈ R and B T ∈ R × . Wedenote by h the maximal element diameter and assume uniform shape regularity.To every element T j , j = , . . . , n T of the mesh, we associate a basis function (cid:101) ψ j for the space (cid:101) Q h with (cid:101) ψ j | T k = δ jk . For every interior edge e i = T l ∩ T r , i = , . . . , n e (cid:98) φ , = (cid:0) − y + y − xy + x (cid:1) , (cid:98) φ , = (cid:0) y − y xy (cid:1) , (cid:98) φ , = (cid:0) − xy − x + x (cid:1) , (cid:98) φ , = (cid:0) − y + xyx − x (cid:1) , (cid:98) φ , = (cid:0) y − y xy − y (cid:1) , (cid:98) φ , = (cid:0) − y + y x + y − xy − (cid:1) , (cid:98) φ , = (cid:0) − x − y + xy + x − x (cid:1) , (cid:98) φ , = (cid:0) − xy + x − x + x (cid:1) . (cid:98) φ , = (cid:0) x (cid:1) , (cid:98) φ , = (cid:0) − y (cid:1) , (cid:98) φ , = (cid:0) − y − y (cid:1) , (cid:98) φ , = (cid:0) x + y − (cid:1) , (cid:98) φ , = (cid:0) − x − y (cid:1) , (cid:98) φ , = (cid:0) xx (cid:1) . Fig. 1
Degrees of freedom and basis functions for the unit triangle and unit square. The black dotsat the vertices represent the quadrature points for the quadrature formula introduced below. of the mesh, we further define two basis functions (cid:101) φ i , (cid:101) φ i + n e which are defined by (cid:101) φ i + (cid:96) · n e | T = B −(cid:62) T (cid:98) φ α , γ , (cid:96) = , , (7)on T ∈ { T l , T r } and vanish identically on all other elements. Here α ∈ { , . . . , (cid:98) n e } refers to the number of the edge e i on the reference element (cid:98) T and γ ∈ { , } dependson (cid:96) and the orientation of the edge e i . The functions (cid:98) φ α , γ are defined in Figure 1.We further set ( a , b ) h , ∗ = ( a , b ) and define ( a , b ) h = ∑ T ( a , b ) h , T with ( a , b ) h , T = | T | ∑ (cid:98) n p l = a ( F T ( (cid:98) x l )) · b ( F T ( (cid:98) x l )) w l , (8)where (cid:98) x l , l = , . . . , (cid:98) n p denote the quadrature points and w l = / (cid:98) n p the quadratureweights on the reference element as depicted in Figure 1.Using the bases defined above, all functions in (cid:101) V h and (cid:101) Q h can be represented as (cid:101) E h = ∑ i (cid:101) e i (cid:101) φ i + (cid:101) e i + n e (cid:101) φ i + n e and (cid:101) H h = ∑ j (cid:101) h j (cid:101) ψ j . (9)This allows to rewrite the variational problem (5)–(6) in algebraic form as (cid:101) M ε ∂ t (cid:101) e = (cid:101) C (cid:62) (cid:101) h (10) (cid:101) M µ ∂ t (cid:101) h = − (cid:101) C (cid:101) e (11) Herbert Egger and Bogdan Radu with matrices ( (cid:101) M ε ) i j = ( ε (cid:101) φ j , (cid:101) φ i ) h , ( (cid:101) M µ ) i j = ( µ (cid:101) ψ j , (cid:101) ψ i ) , and ( (cid:101) C ) i j = ( curl (cid:101) φ j , (cid:101) ψ i ) . Asa direct consequence of the particular choice of the basis functions, we obtain Lemma 1.
Let (cid:101) M ε , (cid:101) M µ , and (cid:101) C be defined as above. Then (i)–(iii) hold analogously.Proof. The properties (i)–(ii) follow directly form the construction. From the par-ticular choice of basis functions, one can deduce that (cid:101) M µ is diagonal and (cid:101) M ε isblock-diagonal; see [5, 11] for details. This implies conditions (iii). (cid:117)(cid:116) Let us mention that the quadrature rule satisfies ( a , b ) h , T = (cid:82) T a ( x ) · b ( x ) dx when a ( x ) · b ( x ) is affine linear. This ensures that the method (5)–(6) also has good approx-imation properties. By a slight adoption of the results given in [5], we obtain Lemma 2.
Let E , H be a smooth solution of (1)–(2) and let (cid:101) E h ( ) and (cid:101) H h ( ) bechosen appropriately. Then (cid:107) (cid:101) E h ( t ) − E ( t ) (cid:107) + (cid:107) (cid:101) H h ( t ) − H ( t ) (cid:107) ≤ Chfor all ≤ t ≤ T with C = C ( E , H , T ) . Moreover, (cid:107) (cid:101) H h ( t ) − π h H ( t ) (cid:107) ≤ Ch where π h H denotes the piecewise constant approximation of H on the mesh T h .Remark 1. For structured meshes and isotropic coefficients, one can observe secondorder convergence also for line integrals of the electric field along edges of themesh. In addition, second convergence can also obtained for unstructured meshesby a non-local post-processing strategy; see [5] for details.
The method of the previous section already yields a stable and efficient approxi-mation. We now show that one degree of freedom per edge can be saved withoutsacrificing the accuracy or efficiency of the method. To this end, we construct ap-proximations E h ( t ) ∈ V h , H h ( t ) ∈ Q h in spaces V h ⊂ (cid:101) V h and Q h = (cid:101) Q h . (cid:98) φ = (cid:0) x (cid:1) , (cid:98) φ = (cid:0) − y (cid:1) , (cid:98) φ = (cid:0) x − (cid:1) , (cid:98) φ = (cid:0) − y (cid:1) . (cid:98) φ = (cid:0) − yx (cid:1) , (cid:98) φ = (cid:0) − yx − (cid:1)(cid:98) φ = (cid:0) − yx (cid:1) . Fig. 2
Degrees of freedom and basis functions on the unit triangle and unit square.
We again define one basis function ψ j of Q h for every element T k by ψ j | T k = δ jk .To any edge e i = T l ∩ T r , we now associate one single basis function φ i defined by mass-lumped mixed finite element method for Maxwell’s equations 5 φ i = (cid:101) φ i + (cid:101) φ i + n e . (12)Using the construction of (cid:101) φ i , one can give an equivalent definition of φ i via φ i | T = B −(cid:62) T (cid:98) φ α , T ∩ e i (cid:54) = /0 , (13)with basis functions (cid:98) φ α = (cid:98) φ α , + (cid:98) φ α , defined on the reference element in Figure 2.Let us note that the space V h coincides with the Nedelec space of lowest order [1].Any function E h ∈ V h and H h ∈ Q h can then be expanded as E h = ∑ i e i φ i and H h = ∑ j h j ψ j . (14)As a consequence of (12), any E h ∈ V h can be interpreted as function (cid:101) E h ∈ (cid:101) V h by E h = ∑ i e i φ i = ∑ i e i ( (cid:101) φ i + (cid:101) φ i + n e ) = ∑ i e i (cid:101) φ i + e i (cid:101) φ i + n e = (cid:101) E h . (15)The coordinates of (cid:101) E h and E h are thus simply connected by (cid:101) e i = (cid:101) e i + n e = e i . Viceversa, we can associate to any function (cid:101) E h ∈ (cid:101) V h a function E h = Π h (cid:101) E h ∈ V h bydefining its coordinates as e i = ( (cid:101) e i + (cid:101) e i + n e ) . In linear algebra notation, this reads e = P (cid:101) e (16)with projection matrix P defined by P i j = if j = i or j = i + n e , and P i j = ( M µ ) i j = ( µψ j , ψ i ) ,C i j = C (cid:48) ji = ( curl φ j , ψ i ) , and M − ε = P (cid:101) M − ε P (cid:62) ,where (cid:101) M ε is defined as in the previ-ous sections. This construction has the folllowing properties. Lemma 3.
Let M µ , C , C (cid:48) , and M − ε be defined as above, and set M ε = ( M − ε ) − .Then the conditions (i)–(iii) are satisfied.Proof. Condition (i) follows by construction. The matrix M µ is diagonal and pos-itive definite and therefore M − µ has the same properties. This verifies (ii) and (iii)for the matrix M µ . Since P is sparse and has fully rank and (cid:101) M − ε is block diagonal,symmetric, and positive definite, one can see that also M − ε is sparse, symmetric,and positive-definite. This verifies conditions (ii) and (iii) for M ε . (cid:117)(cid:116) In the following, we investigate more closely the relation of the system (3)–(4)with matrices as defined above and the system (10)–(11) discussed in the previoussection. We start with an auxiliary result.
Lemma 4.
Let C , P , and (cid:101) C be defined as above. Then one has (cid:101) C = CP .Proof. The result follows directly from the definition of the basis functions. (cid:117)(cid:116)
As a direct consequence, we can reveal the following close connection betweenthe methods (3)–(4) and (10)–(11) discussed in the preceding sections.
Lemma 5.
Let (cid:101) e ( t ) , (cid:101) h ( t ) be a solution of (10)–(11). Then e ( t ) = P (cid:101) e ( t ) , h ( t ) = (cid:101) h ( t ) solves (3)–(4) with matrices M ε , M µ , and C as defined above. Herbert Egger and Bogdan Radu
Proof.
From equation (10), the definition of e , h , and Lemma 4, we deduce that ∂ t e = P ∂ t (cid:101) e = P (cid:101) M − ε (cid:101) C (cid:62) (cid:101) h = P (cid:101) M − ε P (cid:62) C (cid:62) (cid:101) h = M − ε C (cid:62) h . This verifies the validity of equation (3). Using equation (11), we obtainM µ ∂ t h = (cid:101) M µ ∂ t (cid:101) h = − (cid:101) C (cid:101) e = − CP (cid:101) e = − C e , which verifies the validity of equation (4). Finally, using the discrete stability of theprojection completes the proof. (cid:117)(cid:116) Remark 2.
The vectors e ( t ) , h ( t ) computed via (3)–(4) with the above choice of ma-trices correspond to finite element approximations E h ( t ) ∈ V h , H h ( t ) ∈ Q h . There-fore, the procedure described above can be interpreted as a mixed finite elementmethod with mass-lumping based on the approximation spaces V h and Q h .As an immediate consequence of Lemma 5 and the approximation result ofLemma 2, we now obtain the following assertions. Lemma 6.
Let e ( t ) , h ( t ) denote the solutions of (3)–(4) with appropriate initial con-ditions and set E h ( t ) = ∑ i e i ( t ) φ i , H h ( t ) = ∑ j h j ( t ) ψ j . Then (cid:107) E h ( t ) − E ( t ) (cid:107) + (cid:107) H h ( t ) − H ( t ) (cid:107) ≤ Chfor all < t ≤ T . In addition, (cid:107) π h H ( t ) − H h ( t ) (cid:107) ≤ Ch where π h H denotes thepiecewise constant approximation of H on the mesh T h . By some elementary computations, one can verify the following observation.
Lemma 7.
Let T h be a uniform mesh consisting of orthogonal quadrilaterals T ofthe same size. Furthermore, let ε and µ be positive constants. Then the matrices M ε , M µ , and C , defined above coincide with those obtained by the finite differenceapproximation on staggered grids; see [3] for the two dimensional version. The method proposed in this section therefore can be understood as a variationalgeneralization of the Yee scheme. In the two dimensional setting, one degree offreedom e i is required for every edge, and one value h j for every element. Consider the domain Ω = ( − , ) \ { ( x , y ) : ( x − . ) + y ≤ . } , which is splitby an interior boundary into Ω = Ω ∪ Ω ; see Figure 3 for a sketch. We set ε = Ω , ε = Ω and µ = Ω , and consider a plane wave that enters the domainfrom the left boundary. The wave gets slowed down and refracted, when enteringthe domain Ω , and reflected at the circle ∂ Ω , where we enforce a perfect electricboundary conditions. Convergence rates for the numerical solution are depicted inTable 4 and a few snapshots of the solution are depicted Figure 5. mass-lumped mixed finite element method for Maxwell’s equations 7 ∂ Ω ∂ Ω Ω Ω ( − , − ) ( , − )( − , ) ( , ) Fig. 3
Geometry. h DOF ||| E h − π h E h ∗ ||| eoc ||| π h ( H h − π h H h ∗ ) ||| eoc2 − . . − . .
46 0 . . − . .
19 0 . . − . .
08 0 . . Fig. 4
Errors and estimated order of convergence (eoc)with respect to a fine solution ( E h ∗ , H h ∗ ) for h ∗ = − . Thetotal number of degrees of freedom (DOF) is also given. Fig. 5
Snapshots of the post-processed pressure fields (cid:101) p h for time t = . , . , . , . , . Before we conclude, let us briefly discuss an alternative formulation and the exten-sion to three dimensions and higher order approximations.
Remark 3.
Eliminating h from (3)–(4) leads to a second order equationM ε ∂ tt e + K µ − e = e alone, with K µ − = C (cid:48) M − µ C. A sufficient condition forthe stability of the scheme (17) is(iv) M ε and K µ − are symmetric and positive definite, respectively, semi-definite,and for an efficient numerical integration of (17), one now requires that(v) M − ε and K µ − can be applied efficiently.The conditions (iv) and (v) can be seen to be a direct consequence of the conditions(i)–(iii), and the special form K µ − = C (cid:48) M − µ C of the matrix K µ − . Remark 4.
Using the definition of the matrices M µ , C, and C (cid:48) = C (cid:62) given in the pre-vious section, one can verify that K µ − is given by ( K µ − ) i j = ( µ − curl φ j , curl φ i ) .Thus K µ − can be assembled without constructing C or M µ explicitly. Moreover,the conditions (iv) and (v) for K µ − are satisfied automatically. The essential ingre-dient for a mass-lumped mixed finite element approximation of (1)–(2) thus is theconstruction of a positive definite and sparse matrix M − ε . Herbert Egger and Bogdan Radu
Remark 5.
The construction of the approximation M ε discussed in Section 3 imme-diately generalizes to three space dimensions. Like in the two dimensional case, twobasis functions (cid:101) φ i , (cid:101) φ i + n e of the space (cid:101) V h are defined for every edge e i of the mesh andthe approximation ( · , · ) h is defined via numerical quadrature by the vertex rule. Thelumped mass matrix given by ( (cid:101) M ε ) i j = ( ε (cid:101) φ j , (cid:101) φ i ) h then is again block-diagonal. Asbefore, the basis functions for the space V h are then defined by φ i = (cid:101) φ i + (cid:101) φ i + n e andthe inverse mass matrix for the reduced space is again given by M − ε = P (cid:101) M − ε P (cid:62) with projection matrix P of the same form as in two dimensions. Acknowledgements
The authors are grateful for support by the German Research Foundation (DFG)via grants TRR 146, TRR 154, and Eg-331/1-1 and through grant GSC 233 of the“Excellence Initiative” of the German Federal and State Governments.
References
1. Boffi, D., Brezzi, F., Fortin, M.: Mixed finite element methods and applications,