A Weak Galerkin Finite Element Method for the Maxwell Equations

# A Weak Galerkin Finite Element Method for the Maxwell Equations

Lin Mu Department of Mathematics, Michigan State University, East Lansing, MI 48824 (linmu@ msu.edu)    Junping Wang Division of Mathematical Sciences, National Science Foundation, Arlington, VA 22230, USA, jwang@nsf.gov. The research of Wang was supported by the NSF IR/D program, while working at the Foundation. However, any opinion, finding, and conclusions or recommendations expressed in this material are those of the author and do not necessarily reflect the views of the National Science Foundation.    Xiu Ye Department of Mathematics, University of Arkansas at Little Rock, Little Rock, AR 72204, USA, xxye@ualr.edu. This author was supported in part by the National Science Foundation under Grant No. DMS-1115097.    Shangyou Zhang Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA, (szhang@ udel.edu)
###### Abstract

This paper introduces a numerical scheme for time harmonic Maxwell’s equations by using weak Galerkin (WG) finite element methods. The WG finite element method is based on two operators: discrete weak curl and discrete weak gradient, with appropriately defined stabilizations that enforce a weak continuity of the approximating functions. This WG method is highly flexible by allowing the use of discontinuous approximating functions on arbitrary shape of polyhedra and, at the same time, is parameter free. Optimal-order of convergence is established for the weak Galerkin approximations in various discrete norms which are either -like or and -like. An effective implementation of the WG method is developed through variable reduction by following a Schur-complement approach, yielding a system of linear equations involving unknowns associated with element boundaries only. Numerical results are presented to confirm the theory of convergence.

w

eak Galerkin, finite element methods, weak curl, weak gradient, Maxwell equations, polyhedral meshes {AMS} Primary, 65N15, 65N30, 76D07; Secondary, 35B45, 35J50

## 1 Introduction

In this paper, we are concerned with new developments of numerical methods for the time-harmonic Maxwell equations in a heterogeneous medium . The model problem seeks unknown functions and satisfying

 (1) ∇×(μ∇×u)−ϵ∇p = f1inΩ, (2) ∇⋅(ϵu) = g1inΩ, (3) u×n = ϕon∂Ω, (4) p = 0on∂Ω,

where the coefficients and are the magnetic permeability and the electric permittivity of the medium, respectively.

A weak formulation for (1)-(4) seeks such that on and

 (5) (ν∇×u, ∇×v)−(v,∇p) = (f, v),∀v∈H0(curl;Ω) (6) (u,∇q) = −(g,q),∀q∈H10(Ω),

where , and .

The Maxwell equations have been studied extensively in literature by using various numerical methodologies including -conforming edge element approaches [1, 7, 8, 10, 11] and discontinuous Galerkin methods [2, 3, 4, 5, 12, 13]. Particularly in [6], a mixed DG formulation for the problem (1)-(4) was introduced and analyzed. In this DG formulation, both and are approximated by piecewise and functions if is a tetrahedron and by piecewise and if is a parallelepiped, where denotes the set of polynomials of total degree and the set of polynomials of degree in each variable.

The weak Galerkin (WG) finite element method refers to a general finite element technique for partial differential equations where differential operators are approximated as discrete distributions or discrete weak derivatives. The method was first introduced in [15, 16] for second order elliptic equations, and was later extended to other partial differential equations including the Stokes equations [17] and the biharmonic equation [9, 14]. The current research indicates that the concept of discrete weak differential operators offers a new paradigm in numerical methods for partial differential equations.

In this paper, we apply the idea of weak Galerkin to the problem (1)-(4). In essence, this procedure shall introduce a discrete curl operator, which shall be combined with the discrete weak gradient as introduced in [15] to yield a finite element scheme for the Maxwell equations. In this WG method, two types of weak functions are used: and , with and inside of each element and and on the boundary of the element. Error estimates of optimal order are established for the WG approximations in appropriate norms for the case of and , with . For the case of and , only numerical experiments are conducted to illustrate the performance of the corresponding WG finite element scheme; a theoretical study of this WG method is left to interested readers for an investigation.

The use of weak functions and weak derivatives makes the WG method highly flexible on the construction of finite element functions on partitions with arbitrary polygons or polyhedrons. Compared with the DG method in [6], our WG methods make use of additional variables and defined on the boundary of the elements. However, the variables and defined on each element can be eliminated through a local process/computation, yielding a system of linear equations involving only the variables and . Consequently, the WG method has much less number of globally coupled unknowns than DG methods. In addition, the weak Galerkin finite element method is parameter independent in its stability and convergence.

The paper is organized as follows. In Section 2, we introduce some basic notations. In Section 3, we discuss some discrete weak differential operators, particularly a discrete weak curl. Section 4 is devoted to a presentation of the weak Galerkin finite element scheme for the problem (5)-(6). In Section 5, we derive an error equation for the WG finite element approximation. In Section 6, we introduce two types of projection operators and derive some estimates for them. Sections 7 and 8 are devoted to an error analysis for the WG finite element approximations. In Section 9, we discuss an efficient implementation method by using variable reductions/elimination. Finally in Section 10, we present some numerical results that verify the theory established in the previous sections.

## 2 Preliminaries and Notations

Let be any open bounded domain with Lipschitz continuous boundary in . We use the standard definition for the Sobolev space and their associated inner products , norms , and seminorms for any . For example, for any integer , the seminorm is given by

 |v|s,D=⎛⎝∑|α|=s∫D|∂αv|2dD⎞⎠12

with the usual notation

 α=(α1,…,αd),|α|=α1+…+αd,∂α=3∏j=1∂αjxj.

The Sobolev norm is given by

 ∥v∥m,D=(m∑j=0|v|2j,D)12.

The space coincides with , for which the norm and the inner product are denoted by and , respectively. When , we shall drop the subscript in the norm and inner product notation.

The space is defined as the set of vector-valued functions on which, together with their curl, are square integrable; i.e.,

 H(curl;D)={v: v∈[L2(D)]3,∇×v∈[L2(D)]3}.

## 3 Weak Derivatives

The two differential operators used in (5) and (6) are curl and gradient operators. The goal of this section is to introduce an analogy of the curl and gradient operator, called weak curl and weak gradient operators, when the applied functions are discontinuous.

The concept of weak gradient and its discrete analogue was introduced in [15]. This subsection is presented for the sake of completeness of presentation.

Let be any polyhedral domain with boundary . A weak function on the region refers to a function such that and . The first component can be understood as the value of in , and the second component represents on the boundary of . Note that may not necessarily be related to the trace of on should a trace be well-defined. Denote by the space of weak functions on ; i.e.,

 (7) W(K):={v={v0,vb}: v0∈L2(K),vb∈L2(∂K)}.

The weak gradient operator is defined as follows.

###### Definition 3.1

(Weak Gradient) The dual of can be identified with itself by using the standard inner product as the action of linear functionals. With a similar interpretation, for any , the weak gradient of is defined as a linear functional in the dual space of whose action on each is given by

 (8) (∇wv,q)K:=−(v0,∇⋅q)K+⟨vb,q⋅n⟩∂K,

where is the outward normal direction to , is the inner product of and , and is the inner product of and in .

The Sobolev space can be embedded into the space by an inclusion map defined as follows

 iW(ϕ)={ϕ|K,ϕ|∂K},ϕ∈H1(K).

With the help of the inclusion map , the Sobolev space can be viewed as a subspace of by identifying each with .

Let be the set of polynomials on with degree no more than .

###### Definition 3.2

(Discrete Weak Gradient) The discrete weak gradient operator, denoted by , is defined as the unique polynomial satisfying the following equation

 (9) (∇w,r,Kv,q)K=−(v0,∇⋅q)K+⟨vb,q⋅n⟩∂K,∀q∈[Pr(K)]d.

### 3.2 Weak curl and discrete weak curl

To define weak curl, we require weak functions such that and . The first component can be understood as the value of in . The second component represents the value of on the boundary of .

Denote by the space of vector-valued weak functions on ; i.e.,

 (10) V(K)={v={v0,vb}: v0∈[L2(K)]3,vb×n∈[L2(∂K)]3}.

Then, we define a weak curl operator as follows.

###### Definition 3.3

(Weak Curl) The dual of can be identified with itself by using the standard inner product as the action of linear functionals. With a similar interpretation, for any , the weak curl of is defined as a linear functional in the dual space of whose action on each is given by

 (11) (∇w×v,φ)K:=(v0,∇×φ)K+⟨vb×n,φ⟩∂K,

where is the outward normal direction to , is the inner product of and , and is the inner product in .

The Sobolev space can be embedded into the space by an inclusion map defined as follows

 iV(ϕ)={ϕ|K,ϕ|∂K},ϕ∈[H1(K)]3.

Let be any polyhedral domain with boundary . For each face , let and be two assigned unit vectors on the face such that , and are orthogonal each other. Thus, we have . Define . Obviously,. Since the quantity of interest is not but , we will let in order to reduce the number of the unknowns.

###### Definition 3.4

(Discrete Weak Curl) For a given , a discrete weak curl operator, denoted by , is defined as the unique polynomial that satisfies the following equation

 (12) (∇w,r,K×v,φ)K:=(v0,∇×φ)K+⟨vb×n,φ⟩∂K,∀φ∈[Pr(K)]3.

## 4 Numerical Algorithms

Let be a partition of the domain with mesh size that consists of polyhedra of arbitrary shape. Assume that the partition is shape regular in the sense as defined in [16]; i.e. satisfies a set of conditions given in [16]. Denote by the set of all faces in , and let be the set of all interior faces.

Let be shared by two elements and . Let and be two tangential unit vectors on face . For , define a weak Galerkin finite element spaces associated with as

 (13) Vh={v={v0,vb=v1t1+v2t2}: v0|T∈[Pk(T)]3, v1,v2∈Pk(e), e⊂∂T},

and

 (14) Wh={w={w0,wb}: {w0,wb}|T∈Pk−1(T)×Pk(e), e⊂∂T, wb=0 on ∂Ω}.

We also introduce the following subspace of ,

 Vh,0={v={v0,vb}∈Vh,vb×n|e=0, e⊂∂Ω}.

The discrete weak gradient and the discrete weak curl on the finite element spaces and can be computed by using (9) and (12) on each element respectively; i.e.,

 (∇w,kv)|T = ∇w,k,T(v|T),∀v∈Wh (∇w,k−1×v)|T = ∇w,k−1,T×(v|T),∀v∈Vh.

For simplicity of notation, from now on we shall drop the subscript in and in for the discrete weak gradient and the discrete weak curl.

Corresponding to the bilinear forms in (5)-(6), we introduce the following bilinear forms:

 (ν∇w×v, ∇w×w)h = ∑T∈Th(ν∇w×v, ∇w×w)T (v, ∇wq)h = ∑T∈Th(∇wq, v)T.

Furthermore, we stabilize the first one by adding an appropriate stabilization term as follows:

 (15) a(v, w) = (ν∇w×v, ∇w×w)h+s1(v,w),

where

 (16) s1(v,w) = ∑T∈Thh−1⟨(v0−vb)×n,(w0−wb)×n⟩∂T.

For simplicity of notation, we introduce the following notation

 (17) b(v, q) = (v0,∇wq)h

and a second stabilization term

 (18) s2(p,q) = ∑T∈Thh⟨p0−pb,q0−qb⟩∂T.
###### Weak Galerkin Algorithm 1

Find and satisfying on and

 (19) a(uh, v)−b(v,ph) = (f,v0),∀ v={v0,vb}∈Vh,0, (20) b(uh,q)+s2(ph,q) = −(g,q0),∀ q={q0,qb}∈Wh,

where is an approximation of the boundary value in the polynomial space . For simplicity, one may take as the standard projection of the boundary value on each boundary segment.

{lemma}

The weak Galerkin finite element algorithm (19)-(20) has a unique solution.

{proof}

It suffices to show that zero is the only solution of (19)-(20) if and . To this end, assume that the homogeneous conditions are given. Take and in (19)-(20). By adding the two resulting equations, we obtain

 (ν∇w×uh, ∇w×uh)h + ∑T∈Thh−1⟨(u0−ub)×n, (u0−ub)×n⟩∂T + ∑T∈Thh⟨p0−pb, p0−pb⟩∂T=0,

which implies on each , and on . Note that the boundary condition implies on each . Then, it follows from (12) and the integration by parts that for any

 0 = (∇w×uh,v)T = (u0, ∇×v)T+⟨ub×n, v⟩∂T = (∇×u0, v)T+⟨(ub−u0)×n, v⟩∂T = (∇×u0, v)T,

which gives on each . Using (20), (9) and the integration by parts, we have

 0 =

Letting and in the above equation yield on each . Next, by letting and be the jump of on each interior face , we conclude that is continuous across each interior face in the normal direction.

Note that . Thus, there exists a potential function such that on . It follows from and the fact that is continuous that is strongly satisfied in . The boundary condition of (3) implies that on . Therefore, must be a constant on . The uniqueness of the solution of the Laplace equation implies that is the only solution of if is simply connected. Then we must have . Since , we have .

Since , we then have for any . It follows from the definition of and that

 0 = b(v, ph)=(v0,∇wph)h = −∑T∈Th(∇⋅v0, p0)T+∑T∈Th⟨v0⋅n, pb⟩∂T = ∑T∈Th(v0,∇p0)T,

where we have used the fact that on . Letting in (4) gives on each , i.e. is a constant on . Using the facts and on , we obtain .

## 5 Error Equations

For each element , denote by and the projections onto and respectively. Let be the projection onto . Then we can define two projections onto the finite element space and such that on each element ,

 Qhv={Q0v,Qbv=Qb(v1)t1+Qb(v2)t2},Qhq={Q0q,Qbq}.

In addition, denote by the local projection onto . The projection operators , and have some useful properties as stated in the following Lemma.

{lemma}

Let and be the projection operators onto the finite element spaces and respectively. Then, we have

 (22) ∇w×(Qhu)=Qh(∇×u)∀u∈H(curl;Ω)

and

 (23) ∇w(Qhq)=Q0(∇q)∀q∈H1(Ω).
{proof}

Using (12), the integration by parts, and the definition of and , we have

 (∇w×(Qhu),w)T = (Q0u,∇×w)T+⟨(Qbu)×n,w⟩∂T = (u,∇×w)T+⟨u×n,w⟩∂T = (∇×u,w)T=(Qh(∇×u),w)T

for any . This implies that (22) holds true.

As to (23), we use the definition of and the discrete gradient operator to obtain

 (∇w(Qhp),v)T = −(Q0p,∇⋅v)T+⟨Qbp,v⋅n⟩∂T = −(p,∇⋅v)T+⟨p,v⋅n⟩∂T = (∇p,v)T=(Q0(∇p),v)T

for all , which verifies the desired relation (23).

Define two error functions as follows

 (24) eh = {e0,eb}={Q0u−u0,Qbu−ub}, (25) εh = {ε0,εb}={Q0p−p0,Qbp−pb}.

The rest of this section is to derive some equations that the above error functions must satisfy. For simplicity of analysis, we assume that the coefficient in (5) is a piecewise constant function with respect to the finite element partition .

{lemma}

Let be the WG finite element solution arising from (19) and (20), and be the error between the WG finite element solution and the projection of the exact solution as defined in (24)-(25). Then, the following equations are satisfied

 (26) a(eh, v)−b(v, ϵh) = φu(v)∀v∈Vh,0, (27) b(eh, q)+s2(ϵh, q) = ϕu,p(q)∀q∈Wh,

where

 (28) φu(v) = s1(Qhu, v)−l1(u, v), (29) ϕu,p(q) = s2(Qhp, q)+l2(u,q),

and

 (30) l1(u, v) = ∑T∈Th⟨(I−Qh)∇×u, ν(vb−v0)×n⟩∂T (31) l2(u, q) = ∑T∈Th⟨q0−qb, (u−Q0u)⋅n⟩∂T.
{proof}

Using (22), (12), and the integration by parts we have

 (ν∇w×(Qhu),∇w×v)T = (νQh(∇×u),∇w×v)T = (νv0, ∇×Qh(∇×u))T+⟨νvb×n, Qh(∇×u)⟩∂T = (ν∇×v0,Qh(∇×u))T+⟨ν(vb−v0)×n, Qh(∇×u)⟩∂T = (ν∇×u,∇×v0)T+⟨Qh(∇×u),ν(vb−v0)×n⟩∂T.

It follows from (23) that

 (33) (∇w(Qhp),v0)T=(Q0∇p,v0)T=(∇p,v0)T.

Next, using the definition of and , we obtain

 (Q0u, ∇wq)T = −(q0,∇⋅(Q0u))T+⟨qb,Q0u⋅n⟩∂T = (∇q0, u)T−⟨q0−qb, Q0u⋅n⟩∂T.

Testing (1) by with gives

 (35) (∇×(ν∇×u),v0)−(∇p, v0)=(f,v0).

It follows from the integration by parts that

 (∇×(ν∇×u),v0)=∑T∈Th(ν∇×u, ∇×v0)T+∑T∈Th⟨ν(vb−v0)×n, ∇×u⟩∂T,

where we use the fact that . Using (5) and the equation above, we have

 (∇×(ν∇×u),v0) = (ν∇w×(Qhu), ∇w×v)h +∑T∈Th⟨(I−Qh)∇×u, ν(vb−v0)×n⟩∂T.

Substituting (33) and (5) into (35) yields

 (ν∇w×(Qhu), ∇w×v)h−(∇wQhp,v0)h=(f, v0)−l1(v, u).

Adding to both sides of the equation above gives

 (37) a(Qhu, v)−b(v, Qhp)=(f,v0)+φu(v).

To derive a second equation, we test equation (2) by with and then use the integration by parts to obtain

 (38) −∑T∈Th(u,∇q0)T+∑T∈Th⟨u⋅n, q0−qb⟩∂T=(g,q0),

where we have used the fact . Combining (5) with (38) gives

 ∑T∈Th(Q0u,∇wq)T=−(g,q0)+l2(u,q).

Adding to both sides of the equation above gives

 (39) b(Qhu,q)+s2(Qhp, q)=−(g,q0)+ϕu,p(q).

Finally, the differences of (37) and (19), (39) and (20) yield the error equations (26) and (27), respectively.

## 6 Preparation for Error Estimates

For , define as follows

 (40) |||v|||2=a(v,v)=∑T∈Thν∥∇w×v∥2T+∑T∈Thh−1∥(v0−vb)×n∥2∂T.

It is clear that defines merely a semi-norm for the linear space . A norm can be derived from the semi-norm by adding two more terms given as follows

 (41)

where is the jump of the function at each edge/face in the normal direction. The proof of Lemma 4 can be employed to verify that is indeed a norm in . For convenience, we also use the following notation:

 (42) |v|1,h:=⎛⎝∑T∈Thh−1∥(v0−vb)×n∥2∂T⎞⎠1/2.

The linear space can be equipped with the following norm

 |||q|||0=|q|0,h+h∥∇q∥0,h,

where

 |q|20,h=∑T∈Thh∥q0−qb∥2∂T

and

 (43) ∥∇q∥0,h=⎛⎝∑T∈Th∥∇q0∥2T⎞⎠12

for any .

The following Lemma provides some approximation estimates for the projections , , and .

{lemma}

Let be a WG shape regular partition of , , , and . Then, for , we have

 (44) ∑T∈Thh2sT∥w−Q0w∥2s,T≤Ch2(t+1)∥w∥