# Finite element methods with matching and nonmatching meshes for Maxwell equations with disc

SIAM J. NUMER. ANAL. Vol. 37, No. 5, pp. 1542¨C1570

c 2000 Society for Industrial and Applied Mathematics

FINITE ELEMENT METHODS WITH MATCHING AND NONMATCHING MESHES FOR MAXWELL EQUATIONS WITH DISCONTINUOUS COEFFICIENTS?

ZHIMING CHEN? , QIANG DU? , AND JUN ZOU¡ì Abstract. We investigate the ?nite element methods for solving time-dependent Maxwell equations with discontinuous coe?cients in general three-dimensional Lipschitz polyhedral domains. Both matching and nonmatching ?nite element meshes on the interfaces are considered, and optimal error estimates for both cases are obtained. The analysis of the latter case is based on an abstract framework for nested saddle point problems, along with a characterization of the trace space for H (curl ; D), a new extension theorem for H (curl ; D) functions in any Lipschitz domain D, and a novel compactness argument for deriving discrete inf-sup conditions. Key words. ?nite element method, Maxwell equations, interface problems, nonmatching meshes, error estimates, saddle point formulation, extension theorem AMS subject classi?cations. 65M60, 65N30, 35Q60, 46E35 PII. S0036142998349977

1. Introduction. This paper is concerned with the following Maxwell equations in a dielectric medium: (1.1) (1.2) (1.3) (1.4) ¦Å ?t E ? curl H =J in ? ¡Á (0, T ), ? ?t H + curl E =0 in ? ¡Á (0, T ), div(¦Å E) = ¦Ñ in ? ¡Á (0, T ), div(? H) =0 in ? ¡Á (0, T ).

Here ? ? R3 is a simply-connected Lipschitz polyhedral domain with connected boundary which is occupied by the dielectric material. E and H are the electric and magnetic ?elds, and J and ¦Ñ are the current density and charge density. We assume that the permeability parameter ? and the permittivity parameter ¦Å of the medium are discontinuous across an interface ¦£ ? ?, where ¦£ is the boundary of a simply? 1 ? ? and ?2 = ? \ ? ? 1 . ?2 is also connected Lipschitz polyhedral domain ?1 with ? assumed to be simply-connected which, in turn, implies that ¦£ is connected. Without loss of generality we consider only the case with ¦Å and ? being two piecewise constant functions in the domain ?, namely, ¦Å= ¦Å1 ¦Å2 in in ?1 , ?2 , ?= ?1 ?2 in in ?1 , ?2 ,

? Received by the editors January 18, 1999; accepted for publication (in revised form) September 30, 1999; published electronically May 4, 2000. http://www.siam.org/journals/sinum/37-5/34997.html ? Institute of Mathematics, Academia Sinica, Beijing 100080, People¡¯s Republic of China (zmchen@math03.math.ac.cn). The work of this author was partially supported by the China NSF under grant 19771080, China National Key Project ¡°Large Scale Scienti?c and Engineering Computing.¡± ? Department of Mathematics, Hong Kong University of Science and Technology, Clear Water Bay, Hong Kong, and the Department of Mathematics, Iowa State University, Ames, IA 50011 (madu@ust.hk). The work of this author was partially supported by the USNSF under grant DMS9796208 and by a grant from HKRGC. ¡ì Department of Mathematics, The Chinese University of Hong Kong, Shatin, N.T., Hong Kong (zou@math.cuhk.edu.hk). The work of this author was partially supported by Hong Kong RGC grants CUHK 338/96E and CUHK 4004/98P.

1542

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1543

and ¦Åi , ?i (i = 1, 2) are positive constants. It is known that the electric and magnetic ?elds E and H must satisfy the following jump conditions across the interface ¦£ : (1.5) (1.6) [E ¡Á n] = 0, [H ¡Á n] = J¦£ , [¦Å E ¡¤ n] = ¦Ñ¦£ , [? H ¡¤ n] = 0 ,

where n is the unit outward normal to ? ?1 and ¦Ñ¦£ and J¦£ are the surface charge and current density. Throughout the paper, the jump of any function A across the interface ¦£ is de?ned as [A] := A2 |¦£ ? A1 |¦£ with Ai = A|?i , i = 1, 2. We supplement Maxwell equations (1.1)¨C(1.6) with the initial conditions E(x, 0) = E0 (x), H(x, 0) = H0 (x), x¡Ê?.

Instead of solving the fully coupled Maxwell system (1.1)¨C(1.6), we are interested in ?nding the electric and magnetic ?elds separately. To do so, we ?rst eliminate the magnetic ?eld H in (1.1)¨C(1.2) to obtain the electric ?eld equations, (1.7) (1.8) ¦Å ?tt E + curl (??1 curl E) =?t J div (¦Å E) =¦Ñ in in ? ¡Á (0, T ), ? ¡Á (0, T ),

with the following interface and boundary conditions: (1.9) (1.10) [E ¡Á n] = 0, [¦ÅE ¡¤ n] = ¦Ñ¦£ , E ¡Á n = 0 on ? ? ¡Á (0, T ). [??1 curl E ¡Á n] = ??t J¦£ ,

Similarly, we obtain the magnetic ?eld equations from (1.1)¨C(1.6), (1.11) (1.12) ? ?tt H + curl (¦Å?1 curl H) = ? curl (¦Å?1 J) in ? ¡Á (0, T ) , div (? H) = 0 in ? ¡Á (0, T ) ,

with the following interface and boundary conditions: (1.13) [H ¡Á n] = J¦£ , [? H ¡¤ n] = 0, [¦Å?1 curl H ¡Á n] = ?[¦Å?1 J ¡Á n] , (1.14) curl H ¡Á n = ?J ¡Á n on ? ? ¡Á (0, T ). Due to the practical interests, there has been a great deal of work on the numerical approximation to time-dependent Maxwell equations and also on the convergence analysis of numerical schemes for stationary Maxwell equations and related models; see, for example, [1, 3, 10, 13, 18, 20, 25] and the references therein. On the contrary, not too much work was available in the literature concerning the convergence analysis or error estimates for the fully discrete numerical methods for the time-dependent Maxwell systems. For some recent work in this aspect, we refer readers to [12, 19] and [24] for time-dependent Maxwell systems with continuous coe?cients. In this paper we will study a fully discrete ?nite element method for the timedependent Maxwell equations with discontinuous coe?cients. There have been numerous studies on the use of ?nite element methods for elliptic equations having discontinuous coe?cients (see, for example, [5, 11, 14, 7]). For the study presented

1544

ZHIMING CHEN, QIANG DU, AND JUN ZOU

here, we will concentrate on the electric ?eld equations (1.7)¨C(1.10). There is no essential di?culty for the extension of our numerical analysis here to the case for magnetic ?eld equation (1.11)¨C(1.14). We will discretize the system (1.7)¨C(1.10) using an implicit scheme in time and a N? ed? elec¡¯s edge element scheme in space (cf. [21, 22]). In particular, we will investigate ?nite element methods with both matching and nonmatching grids on the interface ¦£. As we will see later, the convergence analysis on the nonmatching case seems to be much more technical than the matching case. For example, a major obstacle is how to ?nd some appropriate ?nite element spaces which satisfy suitable inf-sup conditions. The paper is organized as follows. In section 2, we study the stationary model elliptic problem, (1.15) (1.16) curl (¦Á curl A) + ¦Ã¦Â A = f div (¦Â A) = g in ?, in ?,

with the following interface conditions and boundary conditions: (1.17) (1.18) [A ¡Á n] = 0, A¡Án=0 [¦Â A ¡¤ n] = g¦£ , on ? ?. [¦Á curl A ¡Á n] = f¦£ on ¦£,

Here ¦Á and ¦Â are piecewise constants in ?; i.e., ¦Á = ¦Ái , ¦Â = ¦Âi in ?i for i = 1, 2, and ¦Ái , ¦Âi are positive constants. For technical reasons, ¦Ã is a chosen constant which may be zero in some cases (as in section 3) or strictly positive in some other cases (as in section 4). Such a lower-order term is required to ensure that the corresponding bilinear form induces a full H (curl ; ?)-norm in the latter cases. Finite element methods with matching and nonmatching meshes on the interface will be considered in sections 3 and 4. The results for the above stationary system (1.15)¨C(1.16) are needed in the convergence analysis of ?nite element methods for the time-dependent equations (1.7)¨C(1.10), which will be discussed in section 5. A number of interesting theoretical results concerning the weak formulation and the ?nite element approximation will be stated there while some technical proofs for an extension theorem and an abstract framework for the convergence analysis will be provided in an appendix at the end of the paper. A few conclusion remarks will be furnished in section 6. 2. Stationary model problem. For the convenience of presentation, we ?rst give some notation that will be used throughout the paper: H (curl ; ?) = {v ¡Ê L2 (?)3 ; curl v ¡Ê L2 (?)3 } , H ¦Á (curl ; ?) = {v ¡Ê H ¦Á (?)3 ; curl v ¡Ê H ¦Á (?)3 } (¦Á>0), H0 (curl ; ?) = {v ¡Ê H (curl ; ?); v ¡Á n = 0 on ? ?} . The spaces H (curl ; ?) and H ¦Á (curl ; ?) are equipped with the norms

2 ||v||curl ,? = ||v||2 0,? + ||curl v||0,? 1 /2

, .

2 ||v||¦Á,curl ,? = ||v||2 ¦Á,? + ||curl v||¦Á,?

1 /2

Here and in what follows, ¡¤ 0,? denotes the L2 (?)3 -norm (or the L2 (?)-norm for scalar functions) and ¡¤ ¦Á,? and | ¡¤ |¦Á,? denote the norm and the seminorm of the Sobolev space H ¦Á (?)3 (or H ¦Á (?) for scalar functions). Similar de?nitions are adopted

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1545

for ?1 and ?2 . If no confusion is caused, we will often omit ?, ?1 , and ?2 in the subscripts of the norms. The constant C will always represent a generic constant independent of the time step ¦Ó and the mesh size h, unless otherwise speci?ed. We now discuss the variational formulation and the ?nite element approximation for the stationary model equations (1.15)¨C(1.18). As usual, we introduce a Lagrangian multiplier p to relax the divergence condition (1.16); hence, we consider the following system: (2.1) (2.2) (2.3) (2.4) (2.5) curl (¦Á curl A) + ¦Ã¦Â A ? ¦Â ?p = f in ?, div (¦Â A) = g in ?, [A ¡Á n] = 0, [¦Â A ¡¤ n] = g¦£ on ¦£, [p] = 0, [¦Á curl A ¡Á n] = f¦£ on ¦£, p = 0, A ¡Á n = 0 on ? ?.

Here ¦Ã is a constant which may be taken to be zero unless otherwise stated. Note that if the source terms f and g satisfy div f = ¦Ãg, then we know the solution p of the above system is equal to zero by taking the divergence of (2.1) and using the boundary condition p = 0 on ? ?. This justi?es the introduction of the Lagrangian multiplier p. To establish an appropriate variational formulation for the system (2.1)¨C(2.5), we need a new trace space of H0 (curl ; ?) on the interface ¦£. We next introduce this trace space and present some of its properties for the later use. We know that any v ¡Ê H0 (curl ; ?) has a tangential trace v ¡Á n in H ?1/2 (¦£)3 , de?ned by (2.6) or v ¡Á n, ?

¦£

v ¡Á n, ?

¦£

=

?1

v ¡¤ curl ? dx ?

?1

curl v ¡¤ ?dx ? ? ¡Ê H 1 (?1 )3

=

?2

curl v ¡¤ ?dx ?

?2

v ¡¤ curl ?dx ?? ¡Ê H 1 (?2 )3 ¡É H0 (curl ; ?2 ),

where ¡¤, ¡¤ ¦£ denotes the duality pairing between H ?1/2 (¦£)3 and H 1/2 (¦£)3 (or the duality pairing between H ?1/2 (¦£) and H 1/2 (¦£) for scalar functions) and H0 (curl ; ?2 ) = {v ¡Ê H (curl ; ?2 ) ; v ¡Á n = 0 on ? ?}. However, this characterization of the trace v ¡Á n is rather inconvenient in applications since this trace mapping from H0 (curl ; ?) to H ?1/2 (¦£)3 is not onto. To overcome this di?culty, we introduce the following trace space on ¦£ for functions in H0 (curl ; ?): T (¦£) = {s ¡Ê H ?1/2 (¦£)3 ; ? v ¡Ê H0 (curl ; ?) such that v ¡Á n = s on ¦£}. It is not di?cult to see that T (¦£) is a Banach space under the norm: (2.7) s

T (¦£)

= inf { v

curl ,? ; v

¡Ê H0 (curl ; ?) and v ¡Á n = s on ¦£}.

We note that for domains with smooth boundary, T (¦£) is well known, and it is often denoted by H ?1/2 (div , ¦£) (cf. [8]). For domains with a Lipschitz boundary, we next provide a new characterization of this space.

1546

ZHIMING CHEN, QIANG DU, AND JUN ZOU

For convenience, we denote X1 = H (curl ; ?1 ) and X2 = H0 (curl ; ?2 ). Note that the left-hand side of (2.6) does not make sense if ? belongs to X1 , but the right-hand side of (2.6) is well-de?ned for all v and ? in X1 ; thus, for any s ¡Ê T (¦£), instead of (2.6) we can de?ne (2.8) s, ?

1 ,¦£

=

?1

v ¡¤ curl ? dx ?

?1

curl v ¡¤ ? dx

? ? ¡Ê X1 ,

where v ¡Ê H0 (curl , ?) such that v ¡Á n = s on ¦£. Using the Green¡¯s formula (2.6) and the density of H 1 (?1 ) in X1 , we know that (2.8) is independent of the choice of v ¡Ê H (curl ; ?) such that v ¡Á n = s on ¦£. Thus (2.8) is well-de?ned. Similarly, we can de?ne (2.9) s, ?

2 ,¦£

=

?2

curl v ¡¤ ?dx ?

?2

v ¡¤ curl ?dx

?? ¡Ê X2 ,

where v ¡Ê H0 (curl ; ?) such that s = v ¡Á n on ¦£. It is clear from the de?nitions (2.7)¨C(2.9) that (2.10) s, ?

1 ,¦£

? s, ?

2 ,¦£

¡Ü s

T (¦£)

?

X1 ¡ÁX2

? ? ¡Ê X1 ¡Á X2 ,

which implies that s, ¡¤ 1,¦£ ? s, ¡¤ 2,¦£ de?nes a continuous linear functional in X1 ¡ÁX2 . Hence by Riesz representation theorem, there exists a Q ¡Ê X1 ¡Á X2 satisfying (2.11) (curl Q ¡¤ curl ? + Q ¡¤ ?) dx = s, ?

?1 ¡È?2 1 ,¦£

? s, ?

2 ,¦£

? ? ¡Ê X1 ¡Á X2

and the following identity holds: (2.12) Q

X1 ¡ÁX2

=

?¡ÊX1 ¡ÁX2

sup

s, ?

1 ,¦£

? s, ?

X1 ¡ÁX2

2 ,¦£

?

.

1 (?i )3 Now applying the Green¡¯s formula (2.6) to (2.11) with test functions ? ¡Ê H0 and ? ¡Ê H 1 (?i )3 , respectively, we obtain

curl curl Q + Q = 0 in ?1 ¡È ?2 , curl Q ¡Á n = 0 on ? ?, curl Q ¡Á n = s on ¦£ where the last two relations are understood in the sense of H ?1/2 (? ?)3 and H ?1/2 (¦£)3 , respectively. This yields curl Q ¡Ê H0 (curl ; ?) and (2.13) s

T (¦£)

¡Ü curl Q

curl ,?

= Q

X1 ¡ÁX2 .

Combining (2.10)¨C(2.13) we have proved the following lemma. Lemma 2.1. For any s ¡Ê T (¦£), we have the equality s

T (¦£)

=

?¡ÊX1 ¡ÁX2

sup

s, ?

1 ,¦£

? s, ?

X1 ¡ÁX2

2 ,¦£

?

.

A direct consequence of this lemma is that T (¦£) is indeed a Hilbert space. In fact, if

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1547

si ¡Ê T (¦£) and Qi ¡Ê X1 ¡Á X2 (i = 1, 2) is the corresponding solution of (2.11), then we can de?ne an inner product ((¡¤, ¡¤))¦£ as follows: ((s1 , s2 ))¦£ =

?1 ¡È?2 1 /2

(curl Q1 ¡¤ curl Q2 + Q1 ¡¤ Q2 )dx .

Clearly, s T (¦£) = ((s, s))¦£ . In practice, Lemma 2.1 is rather inconvenient as it uses information from both ?1 and ?2 to de?ne the norm on T (¦£). To overcome the inconvenience, we note that (2.14) |s |1,¦£ = sup ?¡ÊX1 s, ? 1,¦£ , ? X1 |s |2,¦£ = sup ?¡ÊX2 s, ? 2,¦£ ? X2

are also norms of T (¦£). It is clear that |s |1,¦£ ¡Ü s

T (¦£) ,

|s |2,¦£ ¡Ü s

T (¦£)

?s ¡Ê T (¦£).

Both | ¡¤ |1,¦£ and | ¡¤ |2,¦£ are, in fact, equivalent norms on T (¦£) to ¡¤ T (¦£) . To show this, we need the following extension theorem which may be of independent interest. Lemma 2.2 (extension theorem). Assume U is a bounded domain in R3 with a Lipschitz boundary ?U . Let U ?? D. Then there exists a bounded linear operator E : H (curl ; U ) ¡ú H (curl ; R3 ) such that E v = v on U and supp(E v) ? D ? v ¡Ê H (curl ; U ). Because the proof of Lemma 2.2 has not been seen in the literature and the full proof of the result is somewhat involved, we will describe it in more detail in the appendix, section A1. Lemma 2.3. The norms | ¡¤ |1,¦£ and | ¡¤ |2,¦£ are equivalent to ¡¤ T (¦£) . Proof. We need only to prove that |s |1,¦£ ¡Ü C |s |2,¦£ ?s ¡Ê T (¦£). For any ?1 ¡Ê X1 , we use Lemma 2.2 to get a function ? ¡Ê H0 (curl ; ?) such that ? = ?1 on ?1 and ? Note that s, ? Thus, we have

1 ,¦£ X2

¡Ü ?

curl ;?

¡Ü C ? ?1

X1 .

=

s, ?

2 ,¦£

for any ? ¡Ê H0 (curl ; ?) by the Green¡¯s formula. s, ? 2,¦£ = sup ?1 ¡ÊX1 ?1 X1 s, ? 2,¦£ s, ? 2,¦£ = |s |2,¦£ . ¡Ü C ? sup ? curl ;? ?¡ÊX2 ? X2

|s |1,¦£ = sup ?1 ¡ÊX1 ¡Ü C?

s, ?1 ?1

1 ,¦£ X1

sup ?¡ÊH0 (curl ;?)

This completes the proof. 3. Finite element methods with a matching grid. In this section we discuss the ?nite element method based on a matching ?nite element grid on the interface ¦£; i.e., the restrictions from two ?nite element triangulations in ?1 and ?2 match with each other on ¦£. So both triangulations from ?1 and ?2 are combined into a

1548

ZHIMING CHEN, QIANG DU, AND JUN ZOU

standard triangulation of the whole domain ?. In this case, we set ¦Ã = 0 in the system 1 (?) (conforming across (2.1)¨C(2.5) and we require that A ¡Ê H (curl ; ?) and p ¡Ê H0 the interface). This leads to the following weak formulation for the system (2.1)¨C(2.5) with ¦Ã = 0. 1 Find (A, p) ¡Ê H0 (curl ; ?) ¡Á H0 (?) such that (3.1) (3.2) a(A, B) + b(B, p) =

?

f ¡¤ Bdx ? f¦£ , B g qdx + g¦£ , q

? ¦£

2 ,¦£

? B ¡Ê H0 (curl ; ?),

b(A, q ) =

2 3

? q ¡Ê H1 0 (?),

where we assume f ¡Ê L (?) , g ¡Ê L2 (?), f¦£ ¡Ê T (¦£) and g¦£ ¡Ê H ?1/2 (¦£). The two bilinear forms a(¡¤, ¡¤) and b(¡¤, ¡¤) are de?ned by a(B, D) =

?

¦Á curl B ¡¤ curl Ddx ? B, D ¡Ê X, ¦Â B ¡¤ ?qdx

?

b(B, q ) = ?

? B ¡Ê X, q ¡Ê Q,

1 (?). with X = H0 (curl ; ?) and Q = H0 The following compactness result is due to Weber [28]. Lemma 3.1. Let ¦Â be a piecewise constant function in ? and {En } be a sequence in L2 (?)3 satisfying

En

0 ,?

¡Ü C,

curl En

0 ,?

¡Ü C,

div ¦Â En

0 ,?

¡Ü C,

En ¡Á n|? ? = 0

2

for some constant C ; then {En } has a convergent subsequence in L (?)3 . With the above lemma, we may get the following theorem. Theorem 3.2. There exists a unique solution (A, p) ¡Ê X ¡Á Q to (3.1)¨C(3.2). The proof is omitted as it following from a standard argument (see, for example, [21]) based on the Babuska¨CBrezzi theory. We next discuss the ?nite element discretization of (3.1)¨C(3.2). We note that for the matching grid case, studies of ?nite element discretization for the Maxwell systems have been given in [2, 20]. The main error estimates given in this section are, in fact, similar to the existing ones in the literature. We include the discussion here for the sake of completeness, and some of the technical results are also useful for the later discussion on nonmatching grids. Let T h be a shape regular triangulation of ? which matches with the interface ¦£, i.e., K ¡É¦£ = ?

? ?

? K ¡Ê T h.

Here K is the interior of K . Let Sh be the standard continuous piecewise linear ed? elec H (curl , ?)-conforming edge element space ?nite element space and Vh be the N? de?ned by (3.3) Vh = {vh ¡Ê H (curl ; ?); vh = aK + bK ¡Á x on K ? K ¡Ê T h} ,

where aK and bK are two constant vectors. It is known (cf. [21]) that any function vh ¡Ê Vh is uniquely determined by the degrees of freedom in the set ME (v) of the moments on each element K ¡Ê T h , which is given by ME (v) =

e

v ¡¤ ¦Ó ds ;

e is an edge of K

.

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1549

Here ¦Ó is the unit vector along the edge e. For any K ¡Ê T h and p > 2, we denote Yp (K ) = {v ¡Ê Lp (K )3 ; curl v ¡Ê Lp (K )3 ; and v ¡Á n ¡Ê Lp (?K )3 }. We know from [2, Lemma 4.7] that the integrals required in the de?nition of ME (v) make sense for any v ¡Ê Yp (K ). Thus for any v ¡Ê H s (?)3 with curl v ¡Ê Lp (?)3 , where s > 1/2 and p > 2, we can de?ne an interpolation ¦Ðh v such that ¦Ðh v ¡Ê Vh , and ¦Ðh v has the same degrees of freedom as v on each K in T h . Following the same argument as the one used in [12], we can prove the following interpolation error estimate. Lemma 3.3. For 1/2 < s ¡Ü 1, we have v ? ¦Ðh v

0,K

+ curl (v ? ¦Ðh v)

0,K

¡Ü C hs K v

s,curl ,K

? v ¡Ê H s (curl ; K ),

where hK is the diameter of K ¡Ê T h . Furthermore, under the assumptions on the domain ? given in section 1, the interpolation operator ¦Ðh has the following property [21]. 1 Lemma 3.4. Let v = ?p with p ¡Ê H0 (?). Then, if v is regular enough to ensure 1 the existence of ¦Ðh v, we have ¦Ðh v = ?ph for some ph ¡Ê Sh ¡É H0 (?). Now de?ne two ?nite element subspaces of Vh and Sh : Xh = Vh ¡É H0 (curl ; ?),

1 Qh = Sh ¡É H0 (?).

Then the ?nite element approximation of (3.1)¨C(3.2) can be formulated as follows. Find (Ah , ph ) ¡Ê Xh ¡Á Qh such that (3.4) (3.5) a(Ah , Bh ) + b(Bh , ph ) =

?

f ¡¤ Bh dx ? f¦£ , Bh g qh dx + g¦£ , qh

? ¦£

2 ,¦£

? Bh ¡Ê X h ,

b(Ah , qh ) =

? qh ¡Ê Qh .

The following theorem is the main result of this section. Theorem 3.5. Assume the solution (A, p) of problems (3.1)¨C(3.2) has the following regularity: for some 1/2 < s ¡Ü 1, A ¡Ê H s (curl ; ?i ), p ¡Ê H 1+s (?i ), i = 1, 2.

Then there exists a unique solution (Ah , ph ) ¡Ê Xh ¡Á Qh to the discrete problem (3.4)¨C (3.5), and we have the error estimate

2

(3.6)

A ? Ah

curl ,? + p ? ph

1 ,?

¡Ü C hs

i=1

{ A

s,curl ,?i

+ p

1+s,?i }.

To show the theorem, we need the following embedding result proved in [2] for any polyhedral domains. Lemma 3.6. Let XN (?) and XT (?) be the following subspaces of L2 (?)3 : XN (?) = {v ¡Ê L2 (?)3 ; curl v ¡Ê L2 (?)3 , div v ¡Ê L2 (?), v ¡Á n = 0 on ? ?}, XT (?) = {v ¡Ê L2 (?)3 ; curl v ¡Ê L2 (?)3 , div v ¡Ê L2 (?), v ¡¤ n = 0 on ? ?}. If the domain ? is a Lipschitz polyhedral domain in R3 , then there exists a real number r > 1/2 such that XT (?) and XN (?) are continuously embedded into H r (?)3 .

1550

ZHIMING CHEN, QIANG DU, AND JUN ZOU

Proof of Theorem 3.5. Let Nh = {Bh ¡Ê Xh ; b(Bh , qh ) = 0 ? qh ¡Ê Qh }. Then Theorem 3.5 follows immediately from Lemma 3.3, from the coercivity, and from the inf-sup condition below by applying the standard Babuska¨CBrezzi theory (for example, cf. [6]): (3.7) (3.8) sup

Bh ¡ÊXh

a(Bh , Bh ) ¡Ý C0 Bh b(Bh , qh ) ¡Ý C1 qh Bh curl ,?

2 curl ; ? 1 ,?

? Bh ¡Ê N h ,

? qh ¡Ê Qh .

The inf-sup condition (3.8) is a consequence of the Poinc? are inequality by taking Bh = ?qh ¡Ê Xh for any ?xed qh ¡Ê Qh . The coercivity (3.7) was proved in [21] under the assumption that the ?nite element triangulation T h is quasi-uniform. We now 1 give a proof assuming only that T h is shape regular. For any Bh ¡Ê Nh , let ? ¡Ê H0 (?) be the solution of the following problem: ?? ¡¤ ?vdx =

? ? 1 Bh ¡¤ ?vdx ? v ¡Ê H0 (?) .

Set w = Bh ? ??; then it is obvious to see (3.9) curl w = curl Bh , div w = 0 in ?; w ¡Á n = 0 on ? ? .

Hence w ¡Ê XN (?) and so w ¡Ê H r (?) for some r > 1/2 by Lemma 3.6. Now since curl w = curl Bh ¡Ê L¡Þ (?)3 , we know that ¦Ðh w is well-de?ned. Note that ¦Ðh Bh = Bh ; we obtain by Lemma 3.4 that (3.10) Bh = ¦Ðh w + ? ? h for some ?h ¡Ê Qh .

On the other hand, by using Lemma 3.3, (3.9), and the local inverse estimate, we have w ? ¦Ðh w

0,K

¡Ü Chr K ( w r,K + curl w r,K ) = Chr K ( w r,K + curl Bh r,K ) ¡Ü C w r,K + C curl Bh 0,K ,

which yields, by using Lemma 3.6 and (3.9), ¦Ðh w

0 ,?

¡ÜC( w ¡ÜC( w

0 ,?

+ w r,? + curl Bh 0,? ) 0,? + curl w 0,? + curl Bh

0 ,? )

0 ,? )

(3.11)

¡Ü C ( w 0,? + curl Bh ¡Ü C curl Bh 0,? .

Here we have used Lemma 3.1 to conclude that w

0 ,?

¡Ü C curl w

0 ,?

= C curl Bh

0 ,? .

Now the theorem follows by multiplying (3.10) by ¦Â Bh , then using (3.11) and the fact that Bh ¡Ê Nh . 4. Finite element method with a nonmatching grid. The ?nite element method discussed in section 3 for solving the system (2.1)¨C(2.5) is based on a matching ?nite element mesh on the interface ¦£. This imposes a serious restriction, especially

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1551

in three dimensions, on the triangulations in ?1 and ?2 : both must match with each other on ¦£. We are now going to relax this restriction and consider a nonmatching mesh on the interface that allows the two triangulations in ?1 and ?2 to be generated independently. This advantage, however, brings some di?culty to the convergence analysis since the resulting ?nite element spaces will be nonconforming for both the unknown A and Lagrangian multiplier p. Similarly to the treatment of the divergence condition, we will also deal with the constraints [A ¡Á n] = 0 and [p] = 0 on ¦£ by a Lagrangian multiplier approach. As we shall see, the introduction of these new multipliers leads to a nested saddle point problem, for which we need to ?rst generalize the standard saddle point theory (for example, cf. [6, 17]). 4.1. Abstract framework. Let X, Q, M be three real Hilbert spaces and a : X ¡Á X ¡ú R, b : X ¡Á Q ¡ú R, c : Q ¡Á M ¡ú R be three continuous bilinear functionals. Given f ¡Ê X ¡ä , g ¡Ê Q¡ä , and ¦Ö ¡Ê M ¡ä , we consider the following problem: Find (u, p, ¦Ë) ¡Ê X ¡Á Q ¡Á M such that (4.1) (4.2) (4.3) a(u, v ) + b(v, p) = f, v b(u, q ) + c(q, ¦Ë) = g, q c(p, ?) = ¦Ö, ? ? v ¡Ê X, ? q ¡Ê Q, ? ? ¡Ê M.

In the above we have used the same notation ¡¤, ¡¤ to denote all three duality pairings. This saddle point problem and its discretization were considered earlier in [29]. The results to be shown below are more general than the ones in [29] in the sense that our statement on the coercivity and the inf-sup condition for both continuous and discrete cases is more precise and much less restrictive. We will give all the proofs in the appendix, section A2. We ?rst introduce two subspaces N1 ? Q and N2 ? X as follows: (4.4) (4.5) N1 = {q ¡Ê Q; c(q, ?) = 0 ? ? ¡Ê M }, N2 = {v ¡Ê X ; b(v, q ) = 0 ? q ¡Ê N1 }.

Note that in the de?nition of N2 , the test function q is required only to be in N1 , which naturally provides a relaxed condition on the spaces than that given in [29]. We have the following result on the existence and uniqueness of the solution for the system (4.1)¨C(4.3). Lemma 4.1. Assume that a(¡¤, ¡¤) is N2 -coercive, i.e., (4.6) a(v, v ) ¡Ý a0 v

2 X

? v ¡Ê N2 ,

and the following inf-sup conditions hold: (4.7) (4.8) inf sup b(v, q ) ¡Ý b0 , v X q Q c(q, ?) ¡Ý c0 q Q ? M

q ¡ÊN1 v ¡ÊX

?¡ÊM q ¡ÊQ

inf sup

for some positive constants a0 , b0 , c0 . Then the problem (4.1)¨C(4.3) has a unique solution (u, p, ¦Ë) in X ¡Á Q ¡Á M . Next let Xh ? X, Qh ? Q, and Mh ? M be three ?nite dimensional subspaces. We introduce the corresponding approximation of (4.1)¨C(4.3) as follows.

1552

ZHIMING CHEN, QIANG DU, AND JUN ZOU

Find (uh , ph , ¦Ëh ) ¡Ê Xh ¡Á Qh ¡Á Mh such that (4.9) (4.10) (4.11) a(uh , vh ) + b(vh , ph ) = f, vh b(uh , qh ) + c(qh , ¦Ëh ) = g, qh c(ph , ?h ) = ¦Ö, ?h ? vh ¡Ê Xh , ? qh ¡Ê Qh , ? ?h ¡Ê Mh .

Corresponding to (4.4)¨C(4.5) we set (4.12) (4.13) N1h = {qh ¡Ê Qh ; c(qh , ?h ) = 0 N2h = {vh ¡Ê Xh ; b(vh , qh ) = 0 ? ?h ¡Ê Mh }, ? qh ¡Ê N1h }.

We have the following result concerning the discrete problem (4.9)¨C(4.11). Lemma 4.2. Assume the following conditions hold. (i) a(¡¤, ¡¤) is N2h -coercive; i.e., there exists a constant a? > 0 such that (4.14) a(vh , vh ) ¡Ý a? vh

2 X

? vh ¡Ê N2h .

(ii) There exists a constant b? > 0 such that (4.15)

qh ¡ÊN1h vh ¡ÊXh

inf

sup

b(vh , qh ) ¡Ý b? . vh X qh Q

(iii) There exists a constant c? > 0 such that (4.16)

?h ¡ÊMh qh ¡ÊQh

inf

sup

c(qh , ?h ) ¡Ý c? . qh Q ?h M

Then the discrete problem (4.9)¨C(4.11) has a unique solution (uh , ph , ¦Ëh ) in the space Xh ¡Á Qh ¡Á Mh with the following error estimate, u ? uh ¡ÜC

X

+ p ? ph u ? vh

Q X

+ ¦Ë ? ¦Ëh + inf

qh ¡ÊQh

M

vh ¡ÊXh

inf

p ? qh

Q

+ inf

?h ¡ÊMh

¦Ë ? ?h

M

,

where the constant C depends only on a? , b? , c? , and on the operator norms a , b , and c of the bilinear functional a(¡¤, ¡¤), b(¡¤, ¡¤), and c(¡¤, ¡¤), respectively. 4.2. The variational formulation. In this subsection we introduce a new weak formulation for the system (2.1)¨C(2.5) with the constant coe?cient ¦Ã for the lowerorder term being strictly positive. Note that the results of section 3 for matching ?nite element grids are applicable to both stationary and time-dependent Maxwell equations even in the absence of the lower-order term. For nonmatching ?nite element grids, however, without this additional term, we have some technical di?culty in the veri?cation of the coercivity of the bilinear form corresponding to the term curl (¦Á curl A). Though such veri?cation has been made in the matching grid case; see Lemma 4.7. Let us now introduce the following spaces: X1 = H (curl , ?1 ), Q1 = H 1 (?1 ), Then set X = X1 ¡Á X2 , Q = Q1 ¡Á Q2 ¡Á T (¦£), M = H ?1/2 (¦£) X2 = {v ¡Ê H (curl , ?2 ); v ¡Á n = 0 on ? ?}, Q2 = {v ¡Ê H 1 (?2 ); v = 0 on ? ?}.

FINITE ELEMENT METHODS FOR MAXWELL EQUATIONS

1553

and de?ne three bilinear forms a : X ¡Á X ¡ú R, as follows:

2

a(A, B) =

?

¦Ã¦Â A ¡¤ B dx +

i=1 2 ?i

¦Ái curl Ai ¡¤ curl Bi dx ? A, B ¡Ê X, ? s, B1 ? B ¡Ê X, (q, s) ¡Ê Q,

b(B, (q, s)) = ?

i=1

¦Âi ?qi ¡¤ Bi dx + s, B2

?i ¦£

2 ,¦£

1 ,¦£

c((q, s), ?) = q1 ? q2 , ?

? (q, s) ¡Ê Q, ? ¡Ê M.

Now applying the standard technique of integration by parts leads immediately to the following weak formulation of problem (2.1)¨C(2.5). Problem (P). Given (f , g ) ¡Ê L2 (?)4 , ?nd (A, (p, t), ¦Ë) ¡Ê X ¡Á Q ¡Á M such that

2

(4.17) (4.18) (4.19)

a(A, B) + b(B, (p, t)) =

i=1 2 ?i

f ¡¤ Bi dx ? f¦£ , B2 g qi dx + g¦£ , q2

i=1 ?i

2 ,¦£

? B ¡Ê X,

b(A, (q, s)) + c((q, s), ¦Ë) = c((p, t), ?) = 0 ? ? ¡Ê M.

¦£

? (q, s) ¡Ê Q,

The above weak formulation is di?erent from the one used in section 3; it is more consistent with the ?nite element discretization on a nonmatching grid on the interface ¦£. Note that the above system can also be derived based on an optimal control formulation of the interface problem (1.15)¨C(1.18). Similar approaches have been studied for the Poisson-type equations in [14]. The above formulation can be used as a basis for the further development of nonoverlapping domain decomposition methods. First, we have the following result on the existence and uniqueness of solutions to the Problem (P). Theorem 4.3. There exists a unique solution (A, (p, t), ¦Ë) ¡Ê X ¡Á Q ¡Á M to the 1 system (4.17)¨C(4.19). Moreover, (A, p) ¡Ê H0 (curl , ?) ¡Á H0 (?) satis?es (2.1)¨C(2.2) in the sense of distribution, and the Lagrangian multiplier (t, ¦Ë) ¡Ê T (¦£) ¡Á H ?1/2 (¦£) satis?es the following relations: t = ¦Á1 curl A1 ¡Á n = ¦Á2 curl A2 ¡Á n ? f¦£ ¦Ë = ¦Â1 A1 ¡¤ n = ¦Â2 A2 ¡¤ n ? g¦£ in T (¦£), in H ?1/2 (¦£).

Proof. We only prove that the system (4.17)¨C(4.19) has a unique solution; the rest of the theorem can be easily obtained by an appropriate application of the Green¡¯s formula. For the proof, we need to verify the three conditions in Lemma 4.1. First, it is obvious that the bilinear form a : X ¡Á X ¡ú R is N2 -elliptic; i.e., (4.20) a(B, B) ¡Ý a0 B

2 X

? B ¡Ê N2

for some constant a0 > 0. Second, we show that there exists a constant b0 > 0 such that (4.21) sup

B¡ÊX

b(B, (q, s)) ¡Ý b0 (q, s) B X

Q

? (q, s) ¡Ê N1 ¡Á T (¦£).

1 (?). For any given s ¡Ê T (¦£), let Q ¡Ê H (curl , ?1 ) be the We note that N1 = H0 solution of the following problem:

(4.22)

?1

(curl Q ¡¤ curl B + Q ¡¤ B)dx = s, B

1 ,¦£

? B ¡Ê H (curl ; ?1 ).

1554 It is easy to see that (4.23)

ZHIMING CHEN, QIANG DU, AND JUN ZOU

|s |1,¦£ = Q

curl ,?1 .

1 (?) be the solution of the following elliptic interface problem: Now let u ¡Ê H0

(4.24)

?

¦Â ?u ¡¤ ?v dx =

?1

1 ¦Â1 Q ¡¤ ?v dx ? v ¡Ê H0 (?).

1 (?) which satis?es It is clear that (4.24) has a unique solution u ¡Ê H0

(4.25)

?u

0 ,?

¡ÜC Q

0 ,? 1 .

1 With these preparations, for any q ¡Ê H0 (?) and s ¡Ê T (¦£), we de?ne

? = B

??q1 + ?u1 ? Q ??q2 + ?u2

in ?1 , in ?2 ,

? ¡Ê X and where qi = q |?i , ui = u|?i , i = 1, 2. It is obvious that B ? B (4.26)

X

= ? ?q1 + ?u1 ? Q 0,?1 + curl Q ¡Ü C ( ?q 0,? + Q curl ,?1 )

0 ,? 1

+ ? ?q 2 + ? u 2

0 ,? 2

where we have used (4.25). Note that since q1 = q2 , u1 = u2 on ¦£, we have ?(q1 ? q2 ) ¡Á n = 0, ?(u1 ? u2 ) ¡Á n = 0 on ¦£.

Thus, by using (4.24) and (4.22), we are able to obtain ? , (q, s)) = ? b(B

?

¦Â ?q ¡¤ (??q + ?u)dx +

?1 1 ,¦£