High-order accurate discontinuous Galerkin method for the indefinite time-harmonic Maxwell equations

for the indenite time-harmoni Maxwell equations

D. Sármány

∗,

1

, F. Izsák 1

,

2

and J.J.W. van der Vegt 1

De ember 2, 2008

Abstra t

We introdu e a high-order a urate dis ontinuous Galerkin (DG) method for the

indenite frequen y-domain Maxwell equations in three spatial dimensions. The

novelty of the method lies in the way the numeri al ux is omputed. Instead

of using the more popular lo al dis ontinuous Galerkin (LDG) or interior-penalty

dis ontinuous Galerkin (IP-DG) numeri al uxes, we opt for a formulation whi h

makesuseofthelo alliftingoperator. Thisallows usto hooseapenaltyparameter

that is independent of the mesh size and the polynomial order. Moreover, we use

a hierar hi onstru tion of

H(curl)

- onforming basis fun tions, the rst-order version of whi h orrespond to the se ond family of Nédéle elements. We also

provide a priori error boundsfor our formulation, and arry out three-dimensional

numeri al experimentsto validatethetheoreti alresults.

Keywords: numeri al ux with lo al lifting operator; IP-DG method; Maxwell

equations;

H(curl)

- onforming ve tor elements.

1 Introdu tion

The di ulties of solving the Maxwell equations usually lie in the omplexity of the

ge-ometry, the presen e of material dis ontinuities and the fa t that the url operator has

alarge kernel. Moreover, the unknown elds in the Maxwell equationshave spe ial

geo-metri hara teristi s. Theseare mostpronoun edinthethree-dimensionalversionof the

equations, and manifest themselves in the de Rham diagram; see e.g. [8, 21, 29℄.

How-ever, manyofthe popularnumeri aldis retisationte hniquesdonotsatisfy thede Rham

diagram at the dis rete level, and often ontaminate the numeri al solution by

produ -ing spurious modes. One notableex eption are url- onforming nite-element methods,

whi h are spe ial ve tor-valued polynomials that mimi the geometri properties of the

ele tromagneti elds atthe dis retelevel. Based onthe on ept introdu edby Whitney

1

DepartmentofAppliedMathemati s,UniversityofTwente,P.O.Box217,7500AEEns hede,

Nether-lands. E-mail: [d.sarmany,f.izsak,j.j.w.vandervegt℄math.utwente.nl.

2

EötvösLorándUniversity,DepartmentofAppliedAnalysisandComputationalMathemati s,H-1117,

Pázmánysétány1/C,Budapest,Hungary. E-mail: izsakf s.elte.hu.

∗

Nédéle and Bossavit [7, 30, 31℄. A hierar hi onstru tion of high-order basis fun tions

thatsatisfythesamepropertiesaregivenin[1℄fortetrahedralmeshesandin[34℄formore

general three-dimensional meshes. The fa t that these fun tions preserve the geometri

properties of the Maxwell equations has motivated many to study the Maxwell system

and itsnumeri al dis retisationin the framework of dierentialgeometry [9,21℄.

However, su h elements suer from a ouple of pra ti al hurdles. In parti ular,

al-thoughthey are apableof handling omplex geometri al features and material

dis onti-nuities,implementationbe omesin reasinglydi ultwhenhigh-orderbasisfun tionsare

used. Furthermore, extending the approa h to non- onforming mesheswhere the lo al

polynomial order an vary between elements and hanging nodes an be presentposes

onsiderabledi ulties.

One attra tive alternative is the dis ontinuous Galerkin (DG) nite element method.

It an handle non- onforming meshes relatively easily and the implementation of

high-order basis fun tions is also omparatively straightforward. Resear h in the eld of DG

methodshasbeenverya tiveinthepasttenyearsorso;seethere entbooks[14℄ and[20℄

and referen es therein. In the ontext of the Maxwell equations, a nodal approa h was

developed in [18℄, and further studied in [19℄. This approa h had originally been based

onLax-Friedri hstype numeri aluxes, and waslater appliedto the lo aldis ontinuous

Galerkin method [37℄. In the meantime, various DG dis retisations of the low-frequen y

Maxwell equations [23, 24℄ as well as the high-frequen y Maxwell equations [22, 12, 11℄

havealsobeen extensivelystudied. The questionof spuriousmodes inDGdis retisations

has been addressed in [12, 37, 11℄ for onforming meshes and, more re ently, in [13℄ for

two-dimensionalnon- onforming meshes.

Inthiswork,weinvestigatethetime-harmoni Maxwellequationsinalosslessmedium

with inhomogeneous boundary onditions, i.e. nd the (s aled) ele tri eld

E

= E(x)

that satises

∇ ×

1

µ

r

∇ × E − k

2

ε

r

E

= J

in Ω,

n

× E = g

on Γ,

(1)

where

Ω

is an open bounded Lips hitz polyhedron on

R

3

with boundary

Γ = ∂Ω

and outward normal unit ve tor

n

. The right-hand side

J

is the external sour e and

k

is the (real-valued) wave number with the assumption that

k

2

isnot aMaxwell eigenvalue.

Throughout this arti le the (relative) permittivity and the (relative) permeability

orre-spond tova uum(ordry air). That is, weset

ε

r

= 1

and

µ

r

= 1

.

The most important new feature of the high-order DG method dis ussed here is the

ux formulations we apply. In three-dimensional (or, indeed, in any dimensional)

om-putations of the Maxwell equations the most widely used numeri al uxes are the

Lax-Friedri hs ux, the LDG ux and the IP ux. See [20℄ for an overview. While we also

study the omputational performan e of the IP-DG method, the fo us of this arti le is

on a numeri al ux whi h makes use of the lo al lifting operator. This formulation was

originallyintrodu edin [10℄ and further analysedtogether with a large number of other

ux hoi esin [4℄. It has yielded promising results in the dis retisation of the Lapla e

pro eedsalong thelines of [22℄, and isthereforerestri ted tothe ase of smoothmaterial

oe ients. However, we believe that the analysis in [12℄, whi h overs dis ontinuous

materials, an be extended to the DG method presented here. Our theoreti al results

demonstrate the main advantage of the formulation. Namely, that it allows us to use

a stabilising parameter that is independent of both the mesh size and the polynomial

order. This is espe ially important in three-dimensional omputations, where there are

still relatively few experiments available to help us tune a mesh-dependent parameter,

su h as that inthe IP-DG method.

For our DG dis retisation we use a hierar hi onstru tion of

H(curl)

- onforming basisfun tions [1,34℄, whi h satisfythe global de Rhamdiagraminthe ontinuous nite

element setting. However, be ause of the dis ontinuous nature of the method dis ussed

here, we annot expe t our dis retisation to be globally url- onforming and to satisfy

theglobaldeRhamdiagram. Nevertheless, webelievethattheuse of

H(curl)

- onforming basisfun tionisbene ial,sin eitentailsthatthe averagea rossanyfa eisalso

H(curl)

- onforming. Furthermore, the lo al lifting operator is approximated by the same lo al

polynomialbasis as the unknown eld.

Weimplementthebasisfun tionsuptoorderve. Inprin iple,itispossibletoin rease

the order further, but implementation in three dimensions is hindered by a number of

pra ti aldi ulties. First,high-order(i.e.

p > 9

)quadraturerules fortetrahedraare still sub-optimal and omputationally expensive, making the assembly a lengthy pro edure.

Se ond, iterative solvers for indenite linear systems are known to onverge slowly. This

property isexa erbated bythe use ofveryhigh-order

H(curl)

- onformingbasisfun tions, asnding suitablepre onditioners then be omes more of hallenge.

The outlineofthis arti leis asfollows. Wedenethe tessellation and fun tionspa es

inSe tion 2, derive the DG dis retisation for (1) inSe tion 3, and provide a priori error

bounds in Se tion 4. The issue of pre onditioning is very briey addressed in Se tion 5.

We verify and ompare the numeri al methods on both onvex and on ave domains in

Se tion6. Finally, inSe tion 7,we on lude and provide anoutlook.

2 Tessellation and fun tion spa es

We onsider a tessellation

T

h

that partitionsthe polyhedral domain

Ω ⊂ R

3

into aset of

tetrahedra

{K}

. Throughout the arti le we assume that the mesh is shape-regular and that ea htetrahedron is straight-sided. The notations

F

h

,

F

i

h

and

F

b

h

stand respe tively forthe set of allfa es

{F }

, the set of allinternalfa es, andthe set ofall boundaryfa es. Foraboundeddomain

D ⊂ R

d

,

d = 2, 3

,wedenoteby

H

s

_(D)

thestandard Sobolevspa e

of fun tions with regularity exponent

s ≥ 0

and norm

k · k

s,D

. When

D = Ω

, we write

k · k

s

. On the omputationaldomain

Ω

, we introdu e the spa e

H(curl; Ω) :=

n

u

∈

L

2

_(Ω)



3

: ∇ × u ∈

L

2

_(Ω)



3

o

,

with the norm

kuk

2

curl

= kuk

2

0

+ k∇ × uk

2

0

. Let

H

0

(curl; Ω)

denote the subspa e of

H(curl; Ω)

offun tionswithzero tangentialtra e. Wewillalsouse thenotation

(·, ·)

the standard inner produ t in

[L

2

_(D)]

3

,

(u, v)

_D

=

Z

D

u

· v dV,

and the operator

∇

h

for the elementwise appli ation of

∇ = (∂/∂x, ∂/∂y, ∂/∂z)

T

.

We now introdu e the nite element spa e asso iated with the tessellation

T

h

. Let

P

p

(K)

be the spa e of polynomials of degree at most

p ≥ 1

on

K ∈ T

h

. Over ea h element

K

the

H(curl)

- onformingpolynomialspa e is dened as

Q

p

=

u ∈ [P

p

(K)]

3

; u

T

|

s

i

∈ [P

p

(s

i

)]

2

; u · τ

j

|

e

j

∈ P

p

(e

j

) ,

(2) where

s

i

,

i = 1, 2, 3, 4

are the fa es of the element;

e

j

,

j = 1, 2, 3, 4, 5, 6

are the edges of theelement;

u

T

isthetangential omponentof

u

;and

τ

j

isthedire tedtangentialve tor onedge

e

j

. Wedene the spa e

Σ

p

h

as

Σ

p

_h

:=

n

σ

∈ [L

2

(Ω)]

3







σ|

K

∈ Q

p

_{, ∀K ∈ T}

h

o

.

Consider aninterfa e

F ∈ F

h

between element

K

L

and element

K

R

, and let

n

L

and

n

R

represent their respe tive outward pointing normalve tors. We dene the tangential jumpand the average of the quantity

u

a rossinterfa e

F

as

[[u]]

_T

= n

L

× u

L

+ n

R

× u

R

and

{{u}} = u

L

+ u

R

/2,

respe tively. Here

u

L

and

u

R

are the values of the tra e of

u

at

∂K

L

and

∂K

R

,

respe -tively. Atthe boundary

Γ

, we set

{{u}} = u

and

[[u]]

T

= n × u

. In ase weonly need the averageof the tangential omponents, we use the notation

{{u}}

T

. Forthe analysis inSe tion 4,we alsodene the DG norm

kvk

DG

= (kvk

2

0

+ k∇

h

× vk

2

0

+ k

h

−

1

2

[[v]]

T

k

2

0,F

h

)

1

2

,

where

k · k

0,F

h

denotes the the

L

2

_{(F )}

norm, andh

(x) = h

F

, whi h isthe diameteroffa e

F

ontaining

x

. Similarly,

h

K

denotes the diameterof element

K

. Note that the shape-regularproperty of the mesh implies that there is a positive onstant

C

d

independent of the mesh size su h that forall fa es

F

and the asso iated elements

K

R

and

K

L

we have

h

F

≤ C

d

min{h

K

L

, h

_K

R

}.

(3)

To derivethe DG formulations (inthe next se tion)we rst need to introdu e global

liftingoperatorsfor

u

∈ Σ

p

h

. The globallifting operator

L : [L

2

_(F

i

h

)]

3

→ Σ

p

_h

isdened as

(L(u), v)

_Ω

=

Z

F

i

h

u

· [[v]]

_T

dA,

∀v ∈ Σ

p

_h

,

(4)

and the global liftingoperator

R(u) : [L

2

_(F

h

)]

3

→ Σ

p

h

as

(R(u), v)

_Ω

=

Z

F

h

u

· {{v}} dA,

∀v ∈ Σ

p

_h

.

(5)

Fora given fa e

F ∈ F

h

, we willalso need the lo alliftingoperator

R

F

(u) : [L

2

_{(F )]}

3

_→

Σ

p

h

, dened as

(R

F

(u), v)

_Ω

=

Z

F

u

· {{v}} dA,

∀v ∈ Σ

p

_h

.

(6)

Notethat

R

F

(u)

vanishesoutsidetheelements onne tedtothefa e

F

sothatforagiven element

K ∈ T

h

we have the relation

R(u) =

X

F

∈F

h

R

F

(u),

∀u ∈

L

2

(F

h

)



3

.

(7)

3 Dis ontinuous Galerkin dis retisation

We now derive the DG formulation for (1). We rst provide a general bilinear form

wherethe hoi eof thenumeri aluxisnotyetspe ied. Thenwe onsider twodierent

denitions of the numeri alux, ea h of whi hresults in asymmetri algebrai system.

3.1 Derivation of the bilinear form

Thederivationfollowsthesamelinesastheone in[35℄forthe Lapla eoperator. However,

this time it is arried out for the url- url operator. We also refer to [4℄ for a unied

analysis onDG methods for ellipti problems.

We rst introdu e the auxiliaryvariable

q

∈ [L

2

_(Ω)]

3

so that, instead of (1), we an

onsider the rst-order system

∇ × q − k

2

E

= J

in Ω,

q

= ∇ × E

in Ω,

n

× E = g

on Γ.

(8)

From here we follow the standard DG approa h (given, for example, in [4℄ or [35℄ for

ellipti operators): a)integratebothequationsin(8)byparts;b)intheelementboundary

integralssubstitutethe numeri aluxes

q

∗

h

and

E

∗

h

fortheiroriginal ounterparts; ) and nallyintegrateagainthese ondequationin(8)byparts. Thenweseekthepair

(E

h

, q

h

)

su h that for alltest fun tions

(φ, π) ∈ Σ

p

h

× Σ

p

h

:

(q

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

+

X

K∈T

h

(n × q

∗

h

, φ)

∂K

= (J , φ)

Ω

,

(9)

(q

h

, π)

Ω

= (∇

h

× E

h

, π)

_Ω

+

X

K∈T

h

(n × (E

∗

h

− E

h

) , π)

_∂K

.

(10)

Before we pro eed, we make use of the following result: for any given

u, v ∈ Σ

p

h

, the identity

X

K∈T

h

(n × u, v)

_∂K

=

−

Z

F

i

h

{{u}} · [[v]]

_T

dA +

Z

F

i

h

{{v}} · [[u]]

_T

dA +

Z

F

b

h

(n × u) · v dA

(11)

(q

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

i

h

{{q

∗

h

}} · [[φ]]

T

dA

+

Z

F

i

h

{{φ}} · [[q

∗

h

]]

T

dA +

Z

F

b

h

(n × q

∗

h

) · φ dA = (J , φ)

Ω

(12) and

(q

h

, π)

Ω

= (∇

h

× E

h

, π)

_Ω

−

Z

F

i

h

{{E

∗

h

− E

h

}} · [[π]]

T

dA

+

Z

F

i

h

{{π}} · [[E

∗

h

− E

h

]]

_T

dA +

Z

F

b

h

(n × (E

∗

h

− E

h

)) · π dA.

(13)

We an use the liftingoperators to expressand thus eliminatethe auxiliaryvariable

q

h

as a fun tion of

E

h

. From (13) and from the denition of the lifting operators (4) and (5), itfollows that

q

_h

= ∇

h

× E

h

− L({{E

∗

h

− E

h

}}) + R([[E

∗

h

− E

h

]]

T

).

(14) Here we have also used the boundary denition of

[[·]]

T

. Substituting (14) into (12) and applying(10) result in the weak form

B(E

h

, φ) := (∇

h

× E

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

i

h

{{E

∗

_h

− E

h

}} · [[∇

h

× φ]]

T

dA +

Z

F

i

h

[[E

∗

_h

− E

h

]]

T

· {{∇

h

× φ}} dA

−

Z

F

i

h

{{q

∗

h

}} · [[φ]]

T

dA +

Z

F

i

h

[[q

∗

h

]]

T

· {{φ}} dA

+

Z

F

b

h

(n × (E

∗

h

− E

h

)) · (∇

h

× φ) dA −

Z

F

b

h

q

∗

_h

· (n × φ) dA = (J , φ)

_Ω

.

(15)

This is the general primal formulation where one still has freedom to make the hoi es

aboutthenumeri aluxes

E

∗

h

and

q

∗

h

thataremostsuitablefortheproblem. Anoverview of dierent uxes for the Poissonequation is given in[4℄.

3.2 Numeri al uxes

At this point, we spe ify the numeri al uxes

E

∗

h

and

q

∗

h

in (15). We investigate two dierentformulations,one ofwhi hresultsinthe IP-DGformulationthatwasthoroughly

analysed in [22℄. The other is similar to the stabilised entral ux, ex ept that in the

stabilisation term we use the lo al lifting operator (6). Note that in both ases the

numeri al uxes are onsistent, i.e.

∀E, q ∈ H(curl, Ω)

the relations

{{E}}

T

= n × E

,

{{q}}

_T

= n × q

h

,

[[E]]

T

= 0

and

[[q]]

T

= 0

hold. The onsisten y of the DG formulation with the numeri alux of Brezziet al.[10℄ is dis ussed in Appendix A.

First,we dene the numeri al uxes sothat they orrespond to the IP ux,

E

∗

_h

= {{E

h

}} ,

q

∗

h

= {{∇

h

× E

h

}} − τ [[E

h

]]

_T

,

if

F ∈ F

i

h

,

(16)

n

× E

∗

h

= g,

q

∗

h

= ∇

h

× E

h

− τ (n × E

h

) + τ g,

if

F ∈ F

b

h

,

with

τ

being the penalty parameter. We an now transform the following fa e integrals as

Z

F

i

h

[[E

∗

_h

− E

h

]]

T

· {{∇

h

× φ}} dA = −

Z

F

i

h

[[E

h

]]

T

· {{∇

h

× φ}} dA,

Z

F

b

h

(n × (E

∗

_h

− E

h

)) · (∇

h

× φ) dA =

Z

F

b

h

(g − n × E

h

) · (∇

h

× φ) dA,

Z

F

i

h

{{q

∗

h

}} · [[φ]]

T

dA =

Z

F

i

h

{{∇

h

× E

h

}} · [[φ]]

T

dA −

Z

F

i

h

τ [[E

h

]]

T

· [[φ]]

T

dA,

Z

F

b

h

(n × q

∗

h

) · φ dA = −

Z

F

b

h

(∇

h

× E

h

) · (n × φ) dA

+

Z

F

b

h

τ (n × E

h

) · (n × φ) dA −

Z

F

b

h

τ g · (n × φ) dA,

whilethe otherfa e integrals are zero. Ifwe plug these ba k to(15), we havethe IP-DG

method forthe time-harmoni Maxwellequations,

B

h

ip

(E

h

, φ) := (∇

h

× E

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

h

[[E

h

]]

_T

· {{∇

h

× φ}} dA

−

Z

F

h

{{∇

h

× E

h

}} · [[φ]]

T

dA +

Z

F

h

τ [[E]]

_T

· [[φ]]

_T

dA

= (J , φ)

_Ω

−

Z

F

b

h

g

· (∇

h

× φ) dA +

Z

F

b

h

τ g · (n × φ) dA.

(17)

Note that in the left-hand side we no longer distinguish expli itly between internal and

boundary fa es. This is permissible thanks to the denitions of the average and the

tangentialjumpat the boundary.

3.2.2 Numeri al ux of Brezzi et al.

As anext step, we denethe numeri aluxes in the mannerof Brezzi etal. [10℄:

E

∗

_h

= {{E

h

}} ,

q

∗

h

= {{q

h

}} − α

R

([[E

h

]]

T

),

if

F ∈ F

i

h

,

(18)

n

× E

∗

_h

= g,

q

∗

_h

= q

h

− α

R

(n × E

h

) + α

R

(g),

if

F ∈ F

b

h

.

where

α

R

(u) = η

F

{{R

F

(u

h

)}}

for

F ∈ F

h

and

η

F

∈ R

+

. Following the same line of

argument as beforeand using (14), the bilinearform(15) now transformsas

B

h

br

(E

h

, φ) := (∇

h

× E

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

h

[[E

h

]]

_T

· {{∇

h

× φ}} dA −

Z

F

h

{{∇

h

× E

h

}} · [[φ]]

_T

dA

−

Z

F

h

{{R([[E

∗

_h

− E

h

]]

T

)}} · [[φ]]

T

dA +

X

F

∈F

h

Z

F

η

F

{{R

F

([[E

h

]]

T

)}} · [[φ]]

T

dA

+

Z

F

b

h

g

· (∇

h

× φ) dA −

X

F

∈F

b

h

Z

F

η

F

R

F

(g) · (n × φ) dA.

(19)

We an now use the relation

Z

F

h

{{R([[E

∗

h

− E

h

]]

_T

)}} · [[φ]]

_T

dA = (R([[E

∗

h

− E

h

]]

_T

), R([[φ]]

_T

))

_Ω

≈ n

f

X

F

∈F

h

(R

F

([[E

∗

h

− E

h

]]

_T

), R

F

([[φ]]

_T

))

_Ω

= −n

f

X

F∈F

i

h

(R

F

([[E

h

]]

_T

), R

F

([[φ]]

_T

))

_Ω

+ n

f

X

F

∈F

b

h

(R

F

([[g − E

h

]]

_T

), R

F

([[φ]]

_T

))

_Ω

= −n

f

X

F

∈F

h

(R

F

([[E

h

]]

_T

), R

F

([[φ]]

_T

))

_Ω

+ n

f

X

F

∈F

b

h

(R

F

([[g]]

_T

), R

F

([[φ]]

_T

))

_Ω

,

where

n

f

isthe numberof fa esof anelement. Letusintrodu ethebilinearform

B

br

h

: Σ

p

h

×Σ

p

h

→ R

andthelinearform

J

br

h

: Σ

p

h

→ R

as

B

br

h

(E

h

, φ) = (∇

h

× E

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

h

[[E

h

]]

T

· {{∇

h

× φ}} dA

−

Z

F

h

{{∇

h

× E

h

}} · [[φ]]

T

dA +

X

F

∈F

h

(η

F

+ n

f

) (R

F

([[E]]

T

), R

F

([[φ]]

T

))

Ω

,

(20) and

J

h

br

(φ) = (J , φ)

Ω

−

Z

F

b

h

g

· (∇

h

× φ) dA +

X

F

∈F

b

h

(η

F

+ n

f

) (R

F

(g), R

F

(n × φ))

_Ω

,

(21)

respe tively, then the dis rete formulation for the time-harmoni Maxwell equations an

bewritten asfollows. Find

E

h

∈ Σ

p

h

su h that for all

φ

∈ Σ

p

h

the relation

B

h

br

(E

h

, φ) = J

h

br

(φ)

(22)

Following the lines of [22℄ we prove that the dis ontinuous Galerkin dis retisation in

Se tion3 resultsinanoptimal onvergen e rate ofthe numeri alsolutionwith respe t to

the lassi alDG norm. Forsimpli ity,werestri t the analysisto homogeneousboundary

onditions.

A key ingredient inthe error estimation formulais the following Gårdinginequality.

Lemma 1 There exists a onstant

η

, independent of the dis retisation parameter

h =

max

K∈T

h

diam

K

and the wave number

k

, su h that for all

v

∈ Σ

p

h

we have the following inequality

B

h

br

(v, v) ≥ β

2

kvk

2

DG

− (k

2

+ β

2

)kvk

2

0

,

(23) with

β

2

_{= min{}

1

2

,

1

˜

C

inv

}

and

C

˜

inv determinedby inequality (30).

Proof: The righthand side of (23) an berewritten as

β

2

(k∇

h

× vk

2

₀

+ k

h

−

1

2

[[v]]

T

k

2

0,F

h

) − k

2

_kvk

2

0

.

Therefore, using (20) we haveto prove that

k∇

h

× vk

2

0

− 2

Z

F

h

[[v]]

_T

· {{∇

h

× v}} dA +

X

F

∈F

h

(n

f

+ η

F

)kR

F

([[v]]

_T

)k

2

0

≥ β

2

_(k∇

h

× vk

2

0

+ k

h

−

1

2

[[v]]

T

k

2

0,F

h

).

(24)

The se ondterm on the left handside an be estimated with any positive

C

1

as

4

Z

F

h

[[v]]

_T

· {{∇

h

× v}} dA =

Z

F

h

2 · 2

h

−

1

2

C

−1

1

[[v]]

T

·

h

1

2

C

1

{{∇

h

× v}} dA

≤

Z

F

h

4

h

−1

_C

−2

1

| [[v]]

T

|

2

₊

h

C

2

1

| {{∇

h

× v}} |

2

dA

≤ 4C

₁

−2

k

h

−

1

2

[[v]]

T

k

2

0,F

h

+ C

2

1

X

F

∈F

h

h

F

k {{∇

h

× v}} k

2

0,F

.

(25) Next,we use in

Σ

p

h

the inverse inequality

kwk

2

_0,∂K

≤ C

inv

(h

K

)

−1

_kwk

2

0,K

,

(26)

whi h is a straightforward onsequen e of the lo alinverse inequality (Lemma1.38) and

thetra e theorem (B.52)in[14℄. The onstant

C

inv

doesnot depend onthe meshsize, as

we have a shape-regular mesh. Applying (26) inthe se ond term on the right hand side

of (25) for an arbitrary fa e

F

asso iated with elements

K

L

and

K

R

and using (3) we obtain

h

F

C

1

2

k {{∇

h

× v}} k

2

0,F

≤

1

2

h

F

C

2

1

(k∇

h

× v

L

k

2

0,F

+ k∇

h

× v

R

k

2

0,F

)

≤

1

2

C

d

min{h

K

L

, h

K

R

}C

2

1

(k∇

h

× v

L

k

2

0,F

+ k∇

h

× v

R

k

2

0,F

)

≤

1

2

C

2

1

C

inv

C

d

(k∇

h

× vk

2

0,K

L

+ k∇

h

× vk

2

_0,K

R

).

(27)

X

F∈F

h

h

F

C

₁

2

k {{∇

h

× v}} k

2

_0,F

≤ k∇

h

× vk

2

₀

,

(28)

with the hoi e

C

2

1

=

2

C

inv

C

d

n

f

.

(29)

The estimatepre eeding (3.52) in[14℄ givesthat for all

v

∈ Σ

p

h

¯

C

inv

X

F

∈F

h

kR

F

([[v]]

T

)k

2

0

≤ k

h

−

1

2

[[v]]

T

k

2

0,F

h

≤ ˜

C

inv

X

F

∈F

h

kR

F

([[v]]

T

)k

2

0

(30)

with the onstants

C

¯

inv

and

C

˜

inv

independent of

h

. Therefore, using (29) weobtain

C

₁

−2

k

h

−

1

2

[[v]]

T

k

2

0,F

h

=

n

f

C

d

C

inv

2

k

h

−

1

2

[[v]]

T

k

2

0,F

h

≤

n

f

C

d

C

inv

2

C

˜

inv

X

F

∈F

h

kR

F

([[v]]

_T

)k

2

₀

.

(31)

Introdu ing the sum of (28) and (31) in(25) then gives

2

Z

F

h

[[v]]

_T

· {{∇

h

× v}} dA ≤

1

2

k∇

h

× vk

2

0

+ n

f

C

d

C

inv

˜

C

inv

X

F

∈F

h

kR

F

([[v]]

_T

)k

2

₀

,

(32)

and,using (30) we obtain

k∇

h

× vk

2

0

− 2

Z

F

h

[[v]]

_T

· {{∇

h

× v}} dA +

X

F

∈F

h

(n

f

+ η

F

)kR

F

([[v]]

T

)k

2

0

≥

1

2

k∇

h

× vk

2

0

+

X

F

∈F

h

(n

f

+ η

F

− n

f

C

d

C

inv

˜

C

inv

)kR

F

([[v]]

T

)k

2

0

≥ min

F

∈F

h

{

1

2

, n

f

+ η

F

− n

f

C

d

C

inv

˜

C

inv

}(k∇

h

× vk

2

0

+

X

F

∈F

h

kR

F

([[v]]

T

)k

2

0

)

≥ min

F

∈F

h

{

1

2

, n

f

+ η

F

− n

f

C

d

C

inv

˜

C

inv

}(k∇

h

× vk

2

0

+

1

˜

C

inv

k

h

−

1

2

[[v]]

T

k

2

0,F

h

).

Inequality (24) is,therefore satisedif

min

F

∈F

h

(n

f

+ η

F

− n

f

C

d

C

inv

˜

C

inv

) ≥ min

 1

2

,

1

˜

C

inv



,

or,equivalently,

η ≥ n

f

(C

d

C

inv

˜

C

inv

− 1) + min

 1

2

,

1

˜

C

inv



(33)

with

η = min

F

∈F

h

η

F

.



Remark: The onstant

C

inv

in (26) an beestimated using Theorem 4 in[38℄.

In the analysis we onsider now the extended ( f.with (20)) bilinear form

B

h

: (H

0

(

url

, Ω) + Σ

p

h

) × (H

0

(

url

, Ω) + Σ

p

h

) → R,

whi h isgiven as

B

h

(u, v) = (∇

h

× u, ∇

h

× v)

Ω

− k

2

(u, v)

Ω

−

X

F

∈F

h

(R

F

([[u]]

_T

), ∇

h

× v)

Ω

+ (R

F

([[v]]

_T

), ∇

h

× u)

Ω

+

X

F

∈F

h

(n

f

+ η

F

)(R

F

([[u]]

_T

), R

F

([[v]]

_T

))

Ω

and the linear form

J

h

: H

0

(

url

, Ω) + Σ

p

h

→ R

, dened as

J

h

(v) = (J , v)

Ω

.

Lemma 2 Thebilinearform

B

h

is ontinuouson

(H

0

(

url

, Ω) + Σ

p

h

) × (H

0

(

url

, Ω) + Σ

p

h

)

with respe t to the DG norm, i.e. the following inequality holds for all

u

= u

0

+ u

h

and

v

= v

0

+ v

h

with

u

0

, v

0

∈ H

0

(

url

, Ω)

and

u

h

, v

h

∈ Σ

p

h

:

B

h

(u, v) ≤ Ckuk

DG

kvk

DG

.

(34)

Proof: Usingthe result of Lemma 4.8in[22℄ itis su ientto prove that

X

F∈F

h

|(R

F

([[u]]

_T

), R

F

([[v]]

_T

))

Ω

| ≤ Ckuk

DG

kvk

DG

with some mesh-independent onstant

C

. Using the de omposition of

u

and

v

, whi h implies that

[[u

0

]]

T

= [[v

0

]]

T

= 0

, the S hwarz inequality and the estimate in (30) we obtain

X

F

∈F

h

|(R

F

([[u]]

T

), R

F

([[v]]

T

))

Ω

| =

X

F

∈F

h

|(R

F

([[u

h

]]

T

), R

F

([[v

h

]]

T

))

Ω

|

≤

s

X

F

∈F

h

kR

F

([[u

h

]]

_T

)k

2

0

s

X

F

∈F

h

kR

F

([[v

h

]]

_T

)k

2

0

≤ ¯

C

inv

−1

k

h

−

1

2

[[u

h

]]

_T

k

0,F

h

k

h

−

1

2

[[v

h

]]

_T

k

0,F

h

= C

inv

−1

k

h

−

1

2

[[u]]

T

k

0,F

h

k

h

−

1

2

[[v]]

T

k

0,F

h

≤ C

inv

−1

kuk

DG

kvk

DG

.



Next,we willestimatetheresidualterm orrespondingtothe bilinearform

B

h

, whi h isdened for anarbitrary

v

∈ Σ

p

h

, as

r

h

(v) = B

h

(E, v) − J

h

(v),

(35)

where

E

isthe exa t solutionof (1). We an easilyobtain anequivalent form

Lemma 3 Using the notation

Π

h

: [L

2

(Ω)]

3

_{→ Σ}

p

h

for the

L

2

proje tion to the nite element spa e we have for all

v

∈ Σ

p

h

the relation

r

h

(v) =

Z

F

h

[[v]]

_T

· {{∇ × E − Π

h

(∇ × E)}} dA.

(37) If

∇ × E ∈ [H

s

_(Ω)]

3

also holds for some

s

with

s >

1

2

then the following upper bound is valid for all

v

∈ Σ

p

h

:

|r

h

(v)| ≤ Ckvk

DG

s

n

f

X

K∈T

h

h

2 min{p,s}

_K

k∇ × Ek

2

s,K

,

(38)

where

C

does not depend on the mesh size.

Proof: First, we note that for

v

∈ Σ

p

h

∩ H

0

(

url

, Ω)

we have

[[v]]

T

= 0

and therefore, the bilinearform

B

h

(E, v)

simpliesto

B

h

(E, v) = (∇ × E, ∇ × v) − k

2

(E, v).

(39)

Using(35) and the weak formulationof the time-harmoni Maxwellequations we obtain

that

r

h

(v) = B

h

(E, v) − J

h

(v) = (∇ × E, ∇ × v)

Ω

− k

2

(E, v)

Ω

− (J, v) = 0

(40) isvalidforall

v

∈ Σ

p

h

∩ H

0

(

url

, Ω)

. This implies the orthogonalityrelation

B

h

(E − E

h

, v) = 0 ∀v ∈ Σ

p

h

∩ H

0

(

url

, Ω).

The Green formula, the fa t that the tangential omponents of

∇

h

× E

and

{{∇

h

× E}}

oin ide, and the denition of the

L

2

proje tion result in the following relation for an arbitrary

v

∈ Σ

p

h

r

h

(v) = B

h

(E, v) − J

h

(v)

= (∇ × E, ∇

h

× v)

Ω

− k

2

(E, v)

Ω

−

X

F

∈F

h

(R

F

([[v]]

T

), ∇ × E)

Ω

− (J , v)

Ω

= (∇ × ∇ × E, v) +

Z

F

h

∇ × E · [[v]]

_T

dA − k

2

(E, v)

−

X

F

∈F

h

(R

F

([[v]]

_T

), ∇ × E)

Ω

− (J, v)

Ω

=

Z

F

h

{{∇ × E}} · [[v]]

_T

dA −

X

F

∈F

h

(R

F

([[v]]

_T

), Π

h

(∇ × E))

Ω

=

Z

F

h

[[v]]

_T

· {{∇ × E − Π

h

(∇ × E)}}

d

A,

(41) whi h proves (37). If

v

∈ [H

s

_(K)]

3

with

s >

1

2

we use the interpolation estimate, see Theorem 1.7 in[2℄,

wherethe onstant

C

K

isindependentofthemeshsize

h

andthepolynomialorder

p ∈ Σ

p

h

. Using(41), the S hwarz inequality and elementwise summation of (42) givethat

r

h

(v) =

Z

F

h

[[v]]

_T

· {{∇ × E − Π

h

(∇ × E)}}

d

A

=

s

Z

F

h

h

−1

_{| [[v]]}

T

|

2

d

A

s

Z

F

h

h

| {{∇ × E − Π

h

(∇ × E)}} |

2

d

A

≤ kvk

DG

s

n

f

X

K∈T

h

C

K

h

2 min{p,s}

K

k∇ × Ek

2

s,K

.

Using the shape-regular property of the mesh we obtain the upper bound as stated in

(38).



Remark: Using the mesh parameter

h

, the upper bound (38) an be rewritten as

|r

h

(v)| ≤ Ch

min{p,s}

kvk

DG

k∇

h

× Ek

s

.

(43)

Forthe errorestimation we split the omputationalerrorintothree parts.

Lemma 4 Let

η

satisfy the assumption in Lemma 1. Then the omputational error of the solution

E

h

∈ Σ

p

h

of (1) an be estimated as

kE − E

h

k

DG

≤ C

inf

v

∈Σ

p

_h

kE − vk

DG

+ sup

0

₆₌

v

_∈Σ

p

h

r

h

(v)

kvk

DG

+ sup

0

₆₌

v

_∈Σ

p

h

(v, E − E

h

)

Ω

kvk

DG

!

,

(44) where the onstant

C

does not depend on the mesh size.

Proof: Using the triangle inequality and the Gårding inequality stated in Lemma 1 we

obtainthat for anarbitrary

v

∈ Σ

p

h

kE − E

h

k

DG

≤ kE − vk

DG

+ kE

h

− vk

DG

≤ kE − vk

DG

+

1

β

2

sup

06=

w

∈Σ

p

_h

B

h

(E

h

− v, w)

kwk

DG

+

k

2

_{+ β}

2

β

2

sup

06=

w

∈Σ

p

_h

(E

h

− v, w)

kwk

DG

≤ kE − vk

DG

+

1

β

2

sup

06=

w

∈Σ

p

_h

B

h

(E

h

− E, w)

kwk

DG

+

1

β

2

sup

06=

w

∈Σ

p

_h

B

h

(E − v, w)

kwk

DG

+

k

2

_{+ β}

2

β

2

sup

06=

w

∈Σ

p

_h

(E

h

− E, w)

kwk

DG

+

k

2

_{+ β}

2

β

2

sup

06=

w

∈Σ

p

_h

(E − v, w)

kwk

DG

≤ kE − vk

DG

+

1

β

2

sup

06=

w

∈Σ

p

_h

r

h

(w)

kwk

DG

+

C

β

2

kE − vk

DG

+

k

2

_{+ β}

2

β

2

sup

06=

w

∈Σ

p

_h

(E

h

− E, w)

kwk

DG

+

k

2

_{+ β}

2

β

2

kE − vk

DG

= (1 +

C

β

2

+

k

2

_{+ β}

2

β

2

)kE − vk

DG

+

1

β

2

sup

06=

w

∈Σ

p

_h

r

h

(w)

kwk

DG

+

k

2

_{+ β}

2

β

2

sup

06=

w

∈Σ

p

_h

(E

h

− E, w)

kwk

DG

.

Takingthe inmum over all

v

∈ Σ

p

h

we obtain the statementin the lemma.



Usingsmoothnessassumptionstotheexa tsolutionof (1),we annowprovethemain

resultof this se tion.

Theorem 1 Assume that

η

satises the ondition in Lemma 1 and for some parameter

s >

1

2

the exa t solution of (1) satises

E

∈ H

s

(Ω)

and

∇ × E ∈ H

s

(Ω).

Thenfor a meshsize

h

su iently small, we have the followingoptimal error bound

kE − E

h

k

DG

≤ Ch

min{p,s}

_(kEk

s,0

+ k∇ × Ek

s,0

) .

(45)

Proof: We have to estimate the last two terms on the right hand side of (44). For the

estimation of the third term we refer to Proposition 5.2 in [22℄. Note that in the proof

therein the bilinear form

˜a

is investigated on

(H

0

(

url

, Ω) + Σ

p

h

) × (H

0

(

url

, Ω) + Σ

p

h

)

, whi h oin ides with

B

h

onthe same domain. Therefore, the same estimation holds, i.e.

sup

0

₆₌

v

_∈Σ

p

h

(v, E − E

h

)

Ω

kvk

DG

≤ Ch

s

kE − E

h

k

DG

.

(46)

For the estimation of the term

infv

∈Σ

p

_h

kE − vk

DG

we take the Nédéle interpolant

v

=

Π

N,h

E

,whi hsatises the relation

kE − Π

N,h

Ek

DG

= kE − Π

N,h

Ek

url

≤ Ch

min{p,s}

_(kEk

s

+ k∇ × Ek

s

).

(47)

SeeTheorem 5.41and Theorem8.15in[29℄forthe rstand se ondfamilyofNédéle

ele-ments. Using(46),(47)andtheresultin(43)weobtaintheestimateinthe theorem.



Remarks:

1. This result guarantees the same onvergen e rate as Theorem 3.2 in [22℄. The

profound dieren e between the IP method (dis ussed in [22℄) and the approa h

analysedhereisthatwedonothavetousemeshdependent onstantsinthebilinear

form. More pre isely,

η

doesnotdepend on

h

,whi hresultsintheuniformstability of the solution.

2. The upper bound is not sharp if the ratio

max

_K∈Th

h

K

min

_K∈Th

h

K

is large. In those ases one

shouldusethe result ofLemma 3andthe elementwiseinterpolationestimate

orre-spondingto (47).

3. A ording to (46) the mesh size should satisfy the inequality

Ch

s

_{< 1}

.

The dis retisations (17) and (20) result in linear algebrai systems. In this se tion we

simplifynotationand write the linear system as

Ax = b,

(48)

whi h annowdenoteeitherofthelinearsystemsthat resultfromthedis retisations(17)

and (20), respe tively. The system (48) an be solved by a dire tmethod oraniterative

method, typi ally one of the numerous Krylov subspa e methods. Dire t methods are

reliable but not well suited for large-s ale problems be ause one also needs to store the

matrix ll-induring the solution pro ess. On the other hand, e ient iterative methods

require a good pre onditioner. Sin e matrix

A

is symmetri but not denite, nding an e ientpre onditioner isa non-trivialtask. See for example[5, 6,15,33, 36℄.

Theapproa hweadopthere isthatofshifted pre onditionersforse ond-order deriv

a-tiveoperators. Let

S

bethe matrixthat orresponds tothedis retisationof the url- url operator, and let

M

be the dis retisation of the se ond term in (20) (and also in (17)). One an then rewrite(48) as

(S − M)x = b.

(49)

We use pre onditioners that take the form

P = S

B

+ γM

B

,

(50)

where

S

B

onsistsof

N

K

-by-

N

K

diagonalblo ks of

S

,with

N

K

being thelo alnumberof degrees of freedomin element

K

. The positive real number

γ

is the shift parameter and we determine itsvalue experimentally. The matrix

M

B

is either taken tobe the identity matrix

I

or hosen to be

M

itself. Finally, we ompute the Cholesky fa torisation of

P

and use its resulting matri esto apply two-sided pre onditioning to the linear system (48).

In prin iple, it is possible to hoose

S

B

to be

S

, and this is often the ase when

S

represents the dis rete ounterpart of the Lapla ian. However, in the present situation

the Cholesky fa torisationturns out to be omputationally very expensive if

S

B

= S

. The approa h we have just presented here is similar to the one proposed in [36℄ for

the dis reteHelmholtz operator. But the analysis in that work is not dire tly appli able

to the DG dis retisations dis ussed here. Instead, we rely on omputational eviden e

to verify the approa h. On the whole, we have found that, if

γ

is suitably hosen, the pre onditioners speed up the onvergen e of the MINRES method in terms of iteration

steps as well as CPU time. The optimal value of

γ

depends slightly on the mesh-size

h

and more onsiderably on the polynomial order of the approximation

p

. We will briey illustratethis behaviour inSe tion6.3.

6 Numeri al experiments

We now provide three-dimensional numeri al examples to verify and demonstrate the

onstitutethe rst-orderrst-familyofNédéle elements[30℄. Thersttwelveofthe basis

fun tions used here are not the same as those that form the rst-order se ond-family of

Nédéle elements. However, they span exa tlythe samespa e and have exa tlythe same

properties as those dened in [31℄.

In the IP-DG dis retisation,the interior-penalty parameterin(17) is oftendened as

τ = C

IP

p

2

h

F

,

(51)

where

C

IP

= 10

,

p

isthepolynomialorder,and

h

F

isthediameterofthefa e

F ∈ F

h

. The valueofthe onstant

C

IP

isbasedonnumeri alexperiments. Ithasbeendemonstratedfor arange ofproblems that this hoi e of

C

IP

issu iently largeto guarantee stability. See [25℄andreferen estherein. Wehavefound,however, thatinthe aseofthree-dimensional

omputationsfortheMaxwellsystem,thispenaltyisnotlargeenough. Soinsteadof(51),

wedene the above parameter as

τ = C

IP

(p + 1)

2

h

F

,

with

C

IP

= 10

. See alsothe re ent arti le [13℄, where the same penalty term is used for omputingMaxwell eigenvalues ontwo-dimensionaldomains.

The parameter

η

F

(20), resultingfromthe stabilisation term

α

R

in(18), is set to

η

F

= 20,

butthe resultsdonot signi antlydependonthis hoi e,aslongas

η

F

satises ondition (33).

All numeri al omputations have been arried out in the framework of

hp

GEM [32℄, asoftware environment forDG dis retisationssuitableforavarietyof physi alproblems.

6.1 Example 1: smooth solution

As a rst example, we onsider the Maxwell equations (1) with

k

2

_{= 1}

in the domain

Ω = (0, 1)

3

and assumethe boundarytobeaperfe t ele tri ondu tor (PEC),i.e.

g

= 0

in(1). The sour e term is given as

J(x, y, z) = 2π

2

− 1





sin(πy) sin(πz)

sin(πz) sin(πx)

sin(πx) sin(πy)





,

(52)

sowe havethe exa t solution

E(x, y, z) =





sin(πy) sin(πz)

sin(πz) sin(πx)

sin(πx) sin(πy)





.

(53)

The omputations are performed on two dierent sequen es of meshes. The rst are

highlystru tured meshes and onstru ted as follows. The domain

Ω = (0, 1)

3

into

n × n × n

number of ongruent sub ubes, with integer

n = 2

m

and nonnegative

integer

m

. We then divide ea h of these sub ubes into ve tetrahedra, four of whi hare ongruentand havevolumeone-sixthofthe original ube. The fthhas volumeone-third

oftheoriginal ube. Althoughthemeshisnotuniform,thishasproved tobeasimpleand

onvenient way ofmeasuring onvergen e, asea htimewerene themesh,the maximum

ofthefa ediameter

h

F

willbeexa tlyhalfofthatofthepreviousmesh. The onvergen e resultsonstru turedmeshesare showninTableI forthe IP-DGmethodandinTable III

forthe method of Brezzietal.

Wehavealsorunthesameexampleonasequen eofunstru turedmeshes. Themeshes

were generated by CentaurSoft (http://www. entaursoft. om), a pa kage suitable for

generatingavarietyofhybridmesheswith omplexgeometries. Inthissequen eofmeshes,

webeginwitha oarsemeshof

54

tetrahedra. Thenwedivideea htetrahedronintoeight smaller tetrahedra to get the next (ner) mesh. The onvergen e results on stru tured

meshes are depi ted in Table II for the IP-DG dis retisation and in Table IV for the

method ofBrezzi etal.

BasedonouranalysisinSe tion4andontheanalysison[22℄,theoptimal onvergen e

rate for this example is

O(h

p+1

₎

inthe

L

2

_(Ω)

-norm and

O(h

p

₎

inthe DG norm. We an

see that,for both methodsonstru tured meshes,optimal onvergen e rate isa hieved in

the

L

2

_(Ω)

-norm,andhigher-than-optimal onvergen eratesareobservedintheDGnorm.

Onunstru turedmeshes,weonlyhaveanestimated onvergen eratewith

h ∼ N

−

1

3

el

. The onvergen e rates are slightly suboptimal, in part be ause we have toestimate the rates

of onvergen e, and inpart be ause we are stillin the pre-asymptoti regime.

6.2 Example 2: Fi hera ube

As a se ond example, we investigate the DG dis retisations (17) and (22) of the

Maxwell equations (1) with non-smooth solution. To this end, we onsider the domain

Ω = (−1, 1)

3

\ [−1, 0]

3

, and sele t

J

and the non-homogeneous boundary onditions so that the analyti alsolution isgiven by

E

= ∇φ(r),

with

φ(r) = e

r

_sin(r),

(54)

where

r =

px

2

_{+ y}

2

_{+ z}

2

. Theanalyti alsolution ontains asingularityatthe re-entrant

ornerlo atedattheoriginof

Ω

. Asaresult,

E

liesintheSobolevspa e

[H

1−ǫ

_(Ω)]

3

,

ǫ > 0

. Again, based on our analysis in Se tion 4 and on the analysis of [22℄, the theoreti ally

predi ted asymptoti al onvergen e rate for (54)is

O(h

min(1,p+1)

₎

inthe

L

2

_(Ω)

-norm and

O(h

min(1,p)

₎

inthe DG norm.

We ompute the numeri al solutions for a sequen e of globally rened unstru tured

meshes. Theglobalrenementis,of ourse,farfrombeing optimalfor singularsolutions,

but it is nonetheless suitable for verifying the theoreti al results. The onvergen e rates

are again estimated with

h ∼ N

−

1

3

el

.

We show the errors of the IP-DG approximations for (54) in Table V. The errors of

thedis retisation withthe lo alliftingoperatorare depi ted inTable VI. We an expe t

rst-order onvergen e in both norms, whi h is what we approximatelyobserve. Pre ise

onvergen eratesarenota hievedbe ausethemeshesarenotsubsequentlyrened(i.e.we

by the Sobolev oe ientand not by the order of the approximating polynomials.

6.3 Performan e of the iterative solver

To illustrate the ee t of the pre onditioning te hnique briey des ribed in Se tion 5,

we onsider two ases of Example 1 fromSe tion 6.1. One we allmesh

1

4

and it onsists of 27648 elements with polynomial order

p = 1

. The other, mesh

2

3

, onsists of 3456 elementswithpolynomialorder

p = 2

. In TableVIIweshowthe relativeresidualand the omputationalworkafter 10000iterationsfor dierent values of

γ

. We an see that with the orre tvalueof

γ

, the onvergen e rate anbeimproved signi antly. Letusassume, for example, that we want to a hieve a toleran e of tol

= 10

−6

for the given problem

withmesh

1

4

. With

γ = 10

6

ittakes15324iterationsand4990storea hthat toleran e. By

ontrast, the relative residual for the system without pre onditioning is stillonly

0.0085

after 40000 iterations and 6097s. Finally, we note that there is a limit to how far the

parameter

γ

an be in reased without ompromising the solution of the dis rete weak formulation(s). That limit is around

10

10



10

11

, and may thus mean that for high-order

approximation weneed to resort toa less-than-optimal value for

γ

.

7 Con luding remarks and outlook

We have introdu ed and analysed a dis ontinuous Galerkin method for the indenite

time-harmoni Maxwell equations in three dimensions. The novelty of this approa h

is twofold. First, we make use of a lo al lifting operator to ompute the (stabilising)

numeri alux. Thisapproa his mu hinthe spiritof [10℄,and ithas been appliedto the

Lapla e operator in a number of arti les sin e. It also allows us to hoose the onstant

in the ux formulation independent of the mesh size and the polynomial order. Se ond,

we use

H(curl)

- onforming ve tor-valued fun tions to build the lo al polynomial basis, whi h is a very natural hoi e for the Maxwell equations and is widely used in

H(curl)

- onformingnite elementdis retisations [29℄.

We have presented a ouple of numeri alexperiments whi h have demonstrated that

the method onverges at an optimal rate on both stru tured and unstru tured meshes.

Wehavealso arriedout the sameexperimentsusingthe IP-DG uxwiththe same basis

fun tions. Optimal onvergen e rate is, too,observed onboth typesof meshes. However,

for the IP-DG method the penalty parameter depends on the mesh size as well as on

the polynomial order. Furthermore, we found that we need a slightly larger penalisation

than the empiri al value given in [22, 25℄, whi h is mostly based on two-dimensional

experiments. See also[13℄.

It is also lear that a number of questions have remained unanswered. The most

obviousone on ernsthespe tralpropertiesofthemethod. In[12℄theauthorsprovideda

generalframeworkforstudyingthespe tral orre tnessofDGmethods(seealso[11,13℄).

The method introdu ed in this arti le ts into that framework, and the study of its

spe tralproperties is urrently underway.

Another important step would be to provide a posteriori error indi ator for the

A knowledgments

ˆ Thisresear hwassupportedbytheDut hgovernmentthroughthenationalprogram

BSIK:knowledgeandresear h apa ity,intheICTproje t BRICKS(http://www.

bsik-bri ks.nl), theme MSV1.

ˆ We gratefully appre iate the ontributionof Mike Bot hev to Se tion5.

A Appendix: onsisten y

In order to show onsisten y for the bilinear forms resulting from the ux formulation

(18), webegin with the generalprimal formulation

B

h

(E

h

, φ) = (∇

h

× E

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

i

h

{{E

∗

h

− E

h

}} · [[∇

h

× φ]]

T

dA +

Z

F

i

h

[[E

∗

h

− E

h

]]

T

· {{∇

h

× φ}} dA

−

Z

F

i

h

{{q

∗

h

}} · [[φ]]

T

dA +

Z

F

i

h

[[q

∗

h

]]

T

· {{φ}} dA

+

Z

F

b

h

(n × (E

∗

h

− E

h

)) · (∇

h

× φ) dA +

Z

F

b

h

(n × q

∗

h

) · φ dA.

Usingthe identity

(∇

h

× E

h

, ∇

h

× φ)

_Ω

= (∇

h

× ∇

h

× E

h

, φ)

_Ω

−

Z

F

i

h

{{φ}} · [[∇

h

× E

h

]]

T

dA +

Z

F

i

h

{{∇

h

× E

h

}} · [[φ]]

T

dA

+

Z

F

b

h

(n × φ) · (∇

h

× E

h

) dA,

wehave the equivalentformulation

B

h

(E

h

, φ) = (∇

h

× ∇

h

× E

h

, φ)

_Ω

− k

2

(E

h

, φ)

_Ω

−

Z

F

i

h

{{E

∗

h

− E

h

}} · [[∇

h

× φ]]

_T

dA +

Z

F

i

h

[[E

∗

h

− E

h

]]

_T

· {{∇

h

× φ}} dA

−

Z

F

i

h

{{q

∗

h

− ∇

h

× E

h

}} · [[φ]]

_T

dA +

Z

F

i

h

[[q

∗

h

− ∇

h

× E

h

]]

_T

· {{φ}} dA

+

Z

F

b

h

(n × (E

∗

h

− E

h

)) · (∇

h

× φ) dA +

Z

F

b

h

(n × (q

∗

h

− ∇

h

× E

h

)) · φ dA.