Optimal penalty parameters for symmetric discontinuous Galerkin discretisations of the time-harmonic Maxwell equations

Galerkin dis retisations of the time-harmoni Maxwell equations D. Sármány

∗,

1

,

3 , F. Izsák 1

,

2

and J.J.W. vander Vegt 1

January 25,2010

Abstra t

Weprovideoptimalparameterestimatesandapriorierrorboundsforsymmetri

dis ontin-uous Galerkin (DG) dis retisations of the se ond-order indenite time-harmoni Maxwell

equations. More spe i ally, we onsider two variations of symmetri DG methods: the

interiorpenalty DG (IP-DG) method and onethat makes useof the lo al lifting operator

in theuxformulation. Asanovelty,ourparameterestimatesanderrorbounds arei)valid

in thepre-asymptoti regime;ii) solely depend onthegeometry andthepolynomialorder;

and iii) are free of unspe ied onstants. Su h estimates are parti ularly important in

three-dimensional (3D) simulations be ause in pra ti e many 3D omputations o ur in

the pre-asymptoti regime. Therefore, it is vital that our numeri al experiments that

a ompanythetheoreti alresultsarealsoin3D.Theyare arriedoutontetrahedralmeshes

withhigh-order(

p

= 1, 2, 3, 4

)hierar hi

H

(curl)

- onformingpolynomialbasisfun tions.

Keywords: optimal parameter estimates; symmetri dis ontinuous Galerkin

meth-ods; Maxwellequations;

H

(curl)

- onformingve torelements.

1 Introdu tion

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

the presen e of material dis ontinuities and the fa t that the url operator has a large kernel.

Moreover, the unknown elds in the Maxwell equations have spe ial geometri hara teristi s.

These are most pronoun ed in the three-dimensional version of the equations, and manifest

themselvesinthedeRhamdiagram; seee.g.[6 ,17,21 ℄. However,manyofthepopularnumeri al

dis retisation te hniques do not satisfy the de Rham diagram at the dis rete level, and often

ontaminate the numeri al solutionbyprodu ingspurious modes. Onenotable ex eption isthe

H(curl)

- onformingnite-elementmethod,whi hmakesuseofspe ialve tor-valuedpolynomials tomimi thegeometri propertiesoftheele tromagneti eldsatthedis retelevel. Basedonthe

on eptintrodu edbyWhitneyinthe ontextofalgebrai topology [31℄,theywereproposedfor

theMaxwell systembyNédéle and Bossavit[5 , 22 ,23 ℄. Ahierar hi onstru tionof high-order

basisfun tionsthatsatisfythesamepropertiesaregiven in[1 ℄fortetrahedralmeshesandin[27℄

formoregeneralthree-dimensional meshes. Thefa tthatthesefun tionspreserve thegeometri

1

DepartmentofAppliedMathemati s,University ofTwente,P.O.Box217,7500 AEEns hede,Netherlands.

E-mail: [d.sarmany,f.izsak,j.j.w.vand erveg t℄m ath.u twen te.nl .

2

Eötvös Loránd University, Department of Applied Analysis and Computational Mathemati s, H-1117,

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

3

S hoolofComputing,UniversityofLeeds,LS29JT,Leeds,UnitedKingdom. E-mail: d.sarmanyleeds.a .

uk.

∗

numeri al dis retisation intheframework ofdierential geometry [7 ,17℄.

However,su helementssuerfroma oupleofpra ti alhurdles. Inparti ular,althoughthey

are apable of handling omplex geometri al featuresand materialdis ontinuities,

implementa-tion is in reasingly di ult when high-order basis fun tions are used. Furthermore, extending

the approa h to non- onforming mesheswhere the lo al polynomial order an vary between

elementsand hangingnodes an be presentposes onsiderable di ulties.

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

an handlenon- onforming meshesrelatively easily and theimplementation ofhigh-order basis

fun tionsis also omparatively straightforward. Resear h in theeld of DG methods has been

verya tive inthe pasttenyears orso; seethere ent books [13℄and [16 ℄and referen es therein.

In the ontext of the Maxwell equations, a nodal approa h was developed in [14℄, and further

studiedin[15 ℄. Thisapproa hhadoriginallybeenbasedonLax-Friedri hstypenumeri aluxes,

andwaslaterapplied tothe lo aldis ontinuousGalerkinmethod [29℄. Inthemeantime, various

DGdis retisationsofthelow-frequen y Maxwell equations[19, 20℄aswellasthehigh-frequen y

Maxwellequations[18 ,10,9℄havealsobeenextensivelystudied. Thequestionofspuriousmodes

inDGdis retisations hasbeenaddressedin[10 ,29 ,9℄for onformingmeshesand,morere ently,

in[11 ℄for two-dimensional non- onforming meshes.

Inthis work,we investigatethetime-harmoni Maxwell equations ina losslessmediumwith

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

E

= E(x)

thatsatises

∇ ×

_µ

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

withboundary

Γ = ∂Ω

and outward normal unit ve tor

n

. The right-hand side

J

is the external sour e and

k

is the (real-valued) wave numberwiththeassumptionthat

k

2

isnot aMaxwelleigenvalue. Throughoutthis hapter

the (relative) permittivity and the (relative) permeability orrespond to va uum (or dry air).

Thatis, we set

ε

r

= 1

and

µ

r

= 1

.

Out of the many dierent in arnationsof DG dis retisations for (1) we fo uson symmetri

ones, simply be ause they provide the possibility to use linear solvers  su h as MINRES  that

aree ient butonlyappli abletosymmetri matri es. Thesymmetri interiorpenaltyDG

(IP-DG)method isprobablythemost popularsu h method thanksto thesimplepenalisationterm

inthe uxformulation. However, thepenalisation term grows quite sharply as thepolynomial

order is in reased or the mesh is rened. As an alternative, one may opt for a numeri al ux

formulationthatmakesuseof alo alliftingoperator,su hastheonesintrodu ed in[4℄and[8℄.

Theseformulations, togetherwitha large number ofother ux hoi es, were analysedin [2 ℄for

theLapla eoperator,and we refer to thatwork andreferen es thereinfor further details.

Theasymptoti onvergen ebehaviouroftheIP-DGdis retisationfor(1)wasrstestablished

in[18 ℄. In[10 ℄,the asymptoti spe tralproperties oftheasso iated eigenvalueproblem

∇ ×

1

µ

r

∇ × E − k

2

_ε

r

E

= 0 in Ω,

n

_{× E = 0 on Γ,}

(2)

wereanalysed for the IP, in omplete IP, non-symmetri IP, and lo al DG (LDG) methods. An

aprioriestimate forea h of thesemethods results asadire t orollary ofthespe tralanalysis.

Wetakeaslightlydierentapproa h inthis hapter. Iftheproblemisthree-dimensionalitis

[29℄, where it was shown that for a given mesh the dis rete eigenvalues of the symmetri LDG

method tend to the

H(curl)

- onforming dis rete eigenvalues asthe penaltyparameter tendsto innity. The same result isnaturally valid for other symmetri DG dis retisations, su h asthe

ones onsidered here.

However,takingatoolargepenaltyterm omesata omputational ost. Itresultsinalarger

numberofiterationswhenaniterativesolverisusedfor thedis retelinearsystem orresponding

to (1) or (2). Furthermore, if that system is used as a semi-dis rete system in time-domain

omputations, a large penalty term results in a parti ularly stringent time-step restri tion for

expli it time-integration methods. It is therefore essential that an optimal estimate for the

penalty parameter be given that guarantees stability but does not signi antly ompromise

omputationale ien y.

An expli it expression of the IP parameter for the Poisson equations on simpli ial meshes

wasderived in [26℄ andmore re ently in[12 ℄. We extend these results to theMaxwell equation

(1)for IP-DGandalsoprovide anexpli itexpressionoftheDGmethodoriginally introdu edin

[8℄asa slightlymodiedversionof [4 ℄. Ourresults arebasedonthetra einverseinequality[30℄

andon anextension of an a urateestimate for the liftingoperators[25℄.

For ourDGdis retisationwe useahierar hi onstru tionof

H(curl)

- onforming basis fun -tions[1 , 27 ℄. Theysatisfytheglobal deRham diagraminthe ontinuous niteelement setting.

However, be auseof the dis ontinuous natureof the methods dis ussed here, we annot expe t

ourdis retisationtobeglobally

H(curl)

- onformingandtosatisfythedeRhamdiagram. Never-theless,we believethattheuseof

H(curl)

- onforming basisfun tionisbene ial,sin eitentails that the average a ross any fa e is also

H(curl)

- onforming. For higher-order polynomials, it also results in a sparser stiness matrix (i.e. dis rete url- url operator) than standard s alar

H

1

- onforming basisfun tions.

We implement the basis fun tions up to order ve. In prin iple, it is possible to in rease

theorderfurther, but implementation inthree dimensionsis hindered by anumber ofpra ti al

di ulties. First,high-order(i.e.

p > 9

)quadraturerulesfortetrahedraarestillsub-optimaland omputationally expensive, makingthe assembly a lengthy pro edure. Se ond, iterative solvers

forindenite linearsystemsareknownto onvergeslowly,a propertyexa erbated by theuseof

veryhigh-order

H(curl)

- onforming basisfun tions.

The hapterisorganisedasfollows. Wedenethetessellationandfun tionspa esinSe tion2

and derive the DG dis retisation for (1) inSe tion 3. We derive expli it lower bounds for the

penaltyparametersintheDGmethodsandaprioriupperboundsfortheDGmethodsthemselves

inSe tion4. Three-dimensionalnumeri al omputationsare arriedout inSe tion5toshowthe

validityofthe estimates. Finally,inSe tion 6,we on lude andprovide an outlook.

2 Tessellation and fun tion spa es

We onsidera tessellation

T

h

that partitionsthe polyhedral domain

Ω ⊂ R

3

into a setof

tetra-hedra

{K}

. Throughout the hapter we assume that the mesh is shape-regular and that ea h tetrahedronis straight-sided. Thenotations

F

h

,

F

i

h

and

F

b

h

standrespe tively for theset of all fa es

{F }

,theset ofallinternal fa es,and theset ofallboundaryfa es. For aboundeddomain

D ⊂ R

d

,

d = 2, 3

,we denoteby

H

s

_(D)

the standard Sobolev spa e offun tions withregularity

exponent

s ≥ 0

andnorm

k · k

s,D

. When

D = Ω

,we write

k · k

s

. Onthe omputationaldomain

Ω

,weintrodu e the spa e

H(curl; Ω) :=

n

u

_∈

L

2

(Ω)



3

: ∇ × u ∈

L

2

_(Ω)



3

o

,

withthenorm

kuk

2

curl

= kuk

2

0

+ k∇ × uk

2

0

. Let

H

0

(curl; Ω)

denotethesubspa e of

H(curl; Ω)

of fun tionswithzerotangential tra e. Wewill also usethenotation

(·, ·)

D

for thestandard inner produ tin

L

2

_(D)



3

,

(u, v)

_D

=

Z

D

u

_{· v dV,}

andtheoperator

∇

h

for theelementwiseappli 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)

- onforming polynomial spa eisdened as

Q

p

=

n

u

_{∈ [P}

p

(K)]

3

; u

T

|

s

i

∈ [P

p

(s

i

)]

2

; u · τ

j

|

e

j

∈ P

p

(e

j

)

o

,

(3)

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 the element;

u

T

is the tangential omponent of

u

;and

τ

j

is thedire tedtangential ve tor on edge

e

j

. Wedene the spa e

Σ

p

h

as

Σ

p

_h

:=

n

σ

_{∈ [L}

2

(Ω)]

3







σ

|

K

∈ Q

p

, ∀K ∈ T

h

o

.

Consider an interfa 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 normal ve tors. We dene the tangential jump and

theaverage ofthe 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 thevalues of the tra e of

u

at

∂K

L

and

∂K

R

, respe tively.

At theboundary

Γ

,we set

{{u}} = u

and

[[u]]

T

= n × u

. In aseweonly needtheaverageofthe tangential omponents, we usethe notation

{{u}}

T

.

For theanalysisinSe tion 4,wealso dene 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

denotesthe

L

2

_(F)

norm,andh

(x) = h

F

,whi histhediameteroffa e

F

ontaining

x

, i.e.

k

h

−

1

2

[[v]]

T

k

2

0,F

h

=

P

F∈F

h

h

F

[[v]]

T

k

2

0,F

. Similarly,

h

K

denotes the diameter of element

K

. Note that the shape-regular property of themesh implies that there is a positive onstant

C

d

independent of themesh size su h thatfor all fa es

F

and the asso iated elements

K

R

and

K

L

wehave

h

F

≤ C

d

min{h

K

L

, h

_K

R

}.

(4)

To derive the DG formulations (in thenext se tion)werst need to introdu e global lifting

operatorsfor

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

,

(5)

andtheglobal lifting operator

R:

L

2

_(F

h

)



3

→ Σ

p

h

as

(R(u), v)

Ω

=

Z

F

h

u

_{· {{v}} dA,}

_{∀v ∈ Σ}

p

_h

.

(6)

Foragivenfa e

F ∈ F

h

,wewillalsoneedthelo alliftingoperator

R

F

:

L

2

_{(F )}



3

→ Σ

p

h

,dened as

(R

F

(u), v)

_Ω

=

Z

F

u

_{· {{v}} dA,}

_{∀v ∈ Σ}

p

_h

.

(7)

Note that

R

F

(u)

vanishes outside the elements onne ted to the fa e

F

so that for a given element

K ∈ T

h

we have the relation

R(u) =

X

F

∈F

h

R

F

(u),

∀u ∈

L

2

(F

h

)



3

.

(8)

We also use the notation

H

r

_(Ω)

for the Sobolev spa e (with a possibly non-integer exponent)

andthenotation

H

r

_(T

h

) :=

u ∈ L

2

(Ω) : ∇ × u|

K

∈ H

r

(K), ∀K ∈ T

h

,

3 Dis ontinuous Galerkin dis retisation

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

hoi eofthenumeri aluxisnot yetspe ied. Thenwe onsidertwodierentdenitionsofthe

numeri al ux, ea hof whi hresults ina symmetri algebrai system.

3.1 Derivation of the bilinear form

Thederivation follows the same lines astheone in[28℄ for theLapla e operator. However, this

timeit is arriedoutfor the url- url operator. Wealso referto [2 ℄for auniedanalysison DG

methods for ellipti problems.

Werstintrodu etheauxiliaryvariable

q

∈

L

2

_(Ω)



3

sothat, insteadof(1),we an onsider

the rst-order system

∇ × q − k

2

E

= J

in Ω,

q

_{= ∇ × E in Ω,}

(9)

n

_{× E = g on Γ.}

From here we follow the standard DG approa h (given, for example, in [2℄ or [28 ℄ for ellipti

operators): a) multiply both equations in (9) with arbitrary test fun tions

φ

, π ∈ Σ

p

h

and integrateby parts; b)in theelement boundary integrals substitutethe numeri al uxes

q

∗

h

and

E

∗

_h

for their original ounterparts; ) and nally integrateagain the se ond equation in (9) by parts. Thenweseekthepair

(E

h

, q

h

) ∈ Σ

p

h

×Σ

p

h

su hthatforalltestfun tions

(φ, π) ∈ Σ

p

h

×Σ

p

h

:

(q

h

, ∇

h

× φ)

_Ω

− k

2

(E

h

, φ)

_Ω

+

X

K∈T

h

(n × q

∗

h

, φ)

∂K

= (J, φ)

Ω

,

(10)

(q

_h

, π)

_Ω

_{= (∇}

h

× E

h

, π)

_Ω

+

X

K∈T

h

(n × (E

∗

h

− E

h

) , π)

∂K

.

(11)

Beforewepro eed, we makeuse ofthefollowing result: for anygiven

u

, v ∈ Σ

p

h

,theidentity

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

(12)

(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, φ)

Ω

(13) 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.

(14)

We an use the lifting operators to expressand thus eliminatethe auxiliary variable

q

h

as a fun tionof

E

h

. From(14)and fromthedenitionofthelifting operators (5)and(6), itfollows that

q

_h

_{= ∇}

h

× E

h

− L({{E

∗

h

− E

h

}}) + R([[E

∗

h

− E

h

]]

T

).

(15) Herewehave alsousedtheboundary denitionof

[[·]]

T

. Substituting (15)into(13)andapplying (11)results inthe weakform

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}

(16)

+

Z

F

b

h

(n × (E

∗

h

− E

h

)) · (∇

h

× φ) dA −

Z

F

b

h

q

∗

_h

_{· (n × φ) dA = (J, φ)}

_Ω

.

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

numeri al uxes

E

∗

h

and

q

∗

h

that are most suitable for the problem. An overview of dierent uxesfor the Poissonequationis given in[2 ℄.

3.2 Numeri aluxes

At this point, we spe ifythe numeri al uxes

E

∗

h

and

q

∗

h

in (16). We investigate two dierent formulations, one of whi h results in the IP-DG formulation that was thoroughly analysed in

[18℄. The other is similarto the stabilised entral ux, ex eptthat inthestabilisation term we

use the lo al lifting operator (7). Note that in both ases the numeri al uxes are onsistent,

i.e.

∀E, q ∈ H(curl, Ω)

the relations

{{E}}

T

= n × E

,

{{q}} = n × q

h

,

[[E]]

T

= 0

and

[[q]]

T

= 0

hold. 3.2.1 Interior-penalty ux

First,we dene the numeri al uxes sothatthey orrespond to theIPux,

E

∗

_h

_{= {{E}

h

}} ,

q

∗

h

= {{∇

h

× E

h

}} −

a

F

[[E

h

]]

T

,

if

F ∈ F

i

h

,

n

_{× E}

∗

_h

= g,

q

∗

_h

_{= ∇}

h

× E

h

−

a

F

(n × E

h

) +

a

F

g,

if

F ∈ F

b

h

,

(17)

witha

F

beingthe penaltyparameter. We an nowtransform thefollowing fa eintegrals 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

a

F

[[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

a

F

(n × E

h

) · (n × φ) dA −

Z

F

b

h

a

F

g

· (n × φ) dA,

while the other fa e integrals are zero. If we plug these ba k to (16), dene the bilinear form

B

h

ip

: Σ

p

h

× Σ

p

h

→ R

as

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

a

F

[[E]]

T

· [[φ]]

T

dA

(18)

andthelinearform

J

ip

h

: Σ

p

h

→ R

as

J

h

ip

(φ) := (J , φ)

Ω

−

Z

F

b

h

g

_{· (∇}

h

× φ) dA +

Z

F

b

h

a

F

g

· (n × φ) dA,

(19) we have the IP-DG method for the time-harmoni Maxwell equations, formulated as follows.

Find

E

h

∈ Σ

p

h

su hthat forall

φ

∈ Σ

p

h

therelation

B

h

ip

(E

h

, φ) = J

ip

h

(φ)

(20)

issatised. Notethat in(18)we nolonger distinguishexpli itlybetween internal andboundary

fa es. Thisispermissiblethankstothedenitionsoftheaverageandthetangential jumpatthe

boundary.

3.2.2 Numeri alux of Brezzi formulation

Asanext step, we dene the numeri al uxes inthemanner ofBrezzi etal.[8 ℄:

E

∗

_h

_{= {{E}

h

}} ,

q

∗

h

= {{q

h

}} − α

R

([[E

h

]]

_T

),

if

F ∈ F

i

h

,

n

_{× E}

∗

_h

= g,

q

∗

_h

= q

_h

_{− α}

R

(n × E

h

) + α

R

(g),

if

F ∈ F

b

h

.

(21) where

α

R

(u) = η

F

{{R

F

(u

h

)}}

for

F ∈ F

h

and

η

F

∈ R

+

. Following the same line of argument

asbeforeand using (15),the bilinearform (16)nowtransformsas

B(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.

(22)

We annowusetherelation

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), R

F

([[φ]]

T

))

Ω

,

where

n

f

isthe numberof fa esof an element. Letus introdu e the bilinear form

B

br

h

: Σ

p

h

× Σ

p

h

→ R

and thelinearform

J

br

h

: Σ

p

h

→ R

as

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

+

X

F

∈F

h

(η

F

+ n

f

) (R

F

([[E]]

_T

), R

F

([[φ]]

_T

))

_Ω

,

(23) 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 × φ))

Ω

,

(24)

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

writtenasfollows. Find

E

h

∈ Σ

p

h

su hthatfor all

φ

∈ Σ

p

h

therelation

B

br

h

(E

h

, φ) = J

h

br

(φ)

(25)

issatised.

The dis rete ounterparts of the eigenvalue problem (2) for the IP and Brezzi type DG

methods naturally follow from (20) and (25), i.e. nd

k

2

_{∈ R}

+

0

su h that for some

E

h

∈ Σ

p

h

, respe tively,

B

ip

h

(E

h

, φ) = 0

and

B

br

h

(E

h

, φ) = 0

aresatised forall

φ

∈ Σ

p

h

.

4 Expli it parameter and error estimates

Boththe IPand theBrezzi type DG formulations, given respe tively by(20) and (25) , ontain

for these parameters. First, we present an a urate lower bound for thelifting operator

R

F

on tetrahedralelements,extendingtheproofin[25℄forhexahedra. Next,were allthestatementsin

[18℄,whi h arene essary for the onvergen e proof and keep tra kof all onstant terms. Using

these results we provide optimal penalty parameter for both the IP and the Brezzi type DG

method. We also point out that these onditions are su ient for a spurious-free onvergen e

fortheasso iated eigenvalue problems,dis ussed in[10℄.

In the onse utive estimates

K

L

and

K

R

denote theadja ent elements to thefa e

F ∈ F

h

andwe introdu e

M

F

= max



S(F )

V (K

L

₎

,

S(F )

V (K

R

₎



,

where

S

and

V

denote thesurfa eand volume, respe tively.

4.1 Bounds for the lifting operator

Lemma 1 For an arbitrary fa e

F

K

of

K ∈ T

h

any

v

∈ Σ

p

h

satises the inequality

2

3

p

2

F

2

_(p)

S(F

K

)

V (K)

k [[v]]

T

k

2

0,F

K

≤ kR

F

([[v]]

T

)k

2

K

,

(26) where

F

2

_{(p) = 8}

P

p

i=

p

2

1

2i+3

if

p

is even and

F

2

_{(p) =}

8p

2

(p+1)

2

P

p

i=

p−1

2

+1

1

2i+3

if

p

is odd. Proof: The proofis dividedinto three steps.

Step 1  Extension operator on the referen e tetrahedron. We rst onsider a referen e

tetra-hedron

ˆ

K

with verti es

(1, 1, 1), (−1, 1, 1), (1, −1, 1), (1, 1, −1)

and dene an extension opera-tor orresponding to the fa e

F

ˆ

opposite to

(1, 1, 1)

. Let

∆

s

denote a triangle with verti es

(s, 1, 1), (1, s, 1), (1, 1, s)

. Anarbitrary point

(ξ, η, ζ)

an be representedas

(ξ, η, ζ) = (1, s, 1) + u(0, 1 − s, s − 1) + v(s − 1, 1 − s, 0),

(27) where

0 ≤ u, v, u + v ≤ 1

and

−1 ≤ s ≤ 1

, hen e

F = ∆

ˆ

−1

. The Ja obian of the mapping

(ξ, η, ζ) → (u, v, s)

is





0

_{s − 1}

v

1 − s 1 − s 1 − u − v

s − 1

0

u





(28)

withthe determinant

(1 − s)

2

and underthis transformation the fa e

ˆ

F

is mapped to thefa e

˜

F

.

We now dene the extension of the polynomial

φ : ˜

˜

F → R

, whi h is given in terms of the lo al oordinates

(u, v)

. Note that thetransformation

(ξ, η, ζ) → (u, v, s)

islinear from

ˆ

F

to

˜

F

andtherefore

Z

˜

F

| ˜

φ|

2

₌

S( ˜

F )

S( ˆ

F )

Z

ˆ

F

| ˆ

φ|

2

₌

1

4

√

3

Z

ˆ

F

| ˆ

φ|

2

_.

(29)

Iftheorder

p

of thepolynomial

˜

φ

iseven,the extension

ˆ

E( ˜

φ)

isdened as

ˆ

E( ˜

φ)(u, v, s) =

2

p

p

X

j=

p

2

+1

P

_j

(0,2)

_{(−s) ˜}

φ(u, v),

(30) where

P

(0,2)

j

denotesthe

j

th-orderJa obipolynomialon

(−1, 1)

withtheweightfun tion

w(x) =

(1 + x)

2

and

P

(0,2)

j

(1) = 1

. It isalso knownthat

Z

1

−1

(1 + x)

2

P

_i

(0,2)

(x)P

_j

(0,2)

(x) dx =

2

3

_{· Γ(j + 3)Γ(j + 1)}

j! · (2j + 3)Γ(j + 3)

8 δ

ij

=

8 δ

ij

2j + 3

.

Theidentity in(30)givesthat

ˆ

E( ˜

_{φ)(u, v, −1) = ˜}

φ(u, v)

. Interms of

ξ, η, ζ

,we have, using (27) with

φ(u, v) = ˆ

˜

φ(ξ, η)

that

ˆ

E( ˆ

φ)(ξ, η, ζ) = ˆ

φ(ξ, η)

at

F ,

ˆ

hen e

E( ˆ

ˆ

φ)

isinfa tan extension of

φ

ˆ

. Using (28) , (30)and (29), wehave

Z

ˆ

K

| ˆ

E( ˆ

_{φ)(ξ, η, ζ)|}

2

=

Z

1

−1

Z

1

0

Z

1−v

0

| ˆ

E( ˜

_{φ)(u, v, s)|}

2

_{(1 − s)}

2

du dv ds

=

Z

1

−1

Z

1

0

Z

1−v

0













2

p

p

X

i=

p

2

+1

P

_i

(0,2)

_{(−s) ˜}

φ(u, v)













2

(1 − s)

2

du dv ds

(31)

=

4

p

2

Z

1

−1

p

X

i=

p

2

+1

p

X

j=

p

2

+1

P

_i

(0,2)

_(−s)P

_j

(0,2)

_{(−s)(1 − s)}

2

ds

Z

1

0

Z

1−v

0

| ˜

φ(u, v)|

2

_{du dv}

=

4

p

2

p

X

i=

p

2

+1

8

2i + 3

Z

˜

F

| ˜

φ|

2

₌

1

p

2

p

X

i=

p

2

+1

8

√

3

1

2i + 3

Z

ˆ

F

| ˆ

φ|

2

_.

Asaresult, we obtain the relation

k ˆ

E( ˆ

_φ)k

_{0, ˆ}

_K

=

√

4

1

3

1

p

F (p)k ˆ

φk

0, ˆ

F

,

(32) where

F

2

(p) = 8

p

X

i=

p

2

+1

1

2i + 3

if

p

iseven

.

(33)

Analogously,for odd

p

we dene the extension as

ˆ

E( ˜

φ)(u, v, s) =

2

p + 1

p

X

i=

p−1

2

+1

P

_i

(0,2)

_{(−s) ˜}

φ(u, v)

andthesame derivationasin(31) gives that

k ˆ

E( ˆ

_φ)k

2

_{0, ˆ}

_K

=

√

1

3(p + 1)

2

p

X

i=

p−1

2

+1

8

2i + 3

Z

ˆ

F

| ˆ

φ|

2

₌

_√

1

3

F

2

(p)

p

2

k ˆ

φk

2

0, ˆ

F

,

(34)

su hthatwe have

F

2

(p) =

8p

2

(p + 1)

2

p

X

i=

p−1

2

+1

1

2i + 3

if

p

is odd

.

(35)

For omputing the norm of the extension operator

E

ˆ

, both for odd and even

p

, we use the estimates

p

X

i=

p

2

+1

1

2i + 3

≤

Z

p

p

2

1

2t + 3

dt =

1

2

ln

 2p + 3

p + 3



≤

1

₂

ln 2

and

p

X

i=

p−1

2

+1

1

2i + 3

≤

Z

p

p−1

2

1

2t + 3

dt =

1

2

ln

 2p + 3

p + 2



≤

1

₂

ln 2

F

2

_{(p) ≤ 4 ln 2.}

(36) Theestimatein(36) issharpas

lim p → ∞

.

Step 2  Extension operator on a general tetrahedron. For an arbitrary tetrahedron

K

with a fa e

F

K

we dene the anetransformation

T

K

: ˆ

K → K

as

T

K

(ˆ

x) = J

K

x

ˆ

+ b,

where

b

∈ R

3

_{, J}

K

∈ R

3×3

and

T

K

( ˆ

F ) = F

K

.

Theextension

E

of a fun tion

φ : F

K

→ R

is given thenasfollows:

ˆ We dene thefun tion

φ : ˆ

ˆ

F → R

with

ˆ

φ(ˆ

x) := φ(T

K

x).

ˆ

ˆ We extend

ˆ

φ

to

E( ˆ

ˆ

φ)

using the methodinStep 1. ˆ The extension to

K

is given by

E(φ)(x) := ˆ

E( ˆ

φ)(T

_K

−1

x).

As

J

K

islinear,we anapplyasimple hangeofvariables

x

= T

K

(ˆ

x)

for omputingtheintegral ofany

g ∈ L

1

(K)

:

Z

K

g(x) = |

det

J

K

|

Z

ˆ

K

ˆ

g(ˆ

x) =

V (K)

V ( ˆ

K)

Z

ˆ

K

ˆ

g(ˆ

x).

(37)

Sin etherestri tion of

J

K

to the fa e

F

K

of

K

remains ane, we also have, asin(29), that

Z

F

K

g(x) =

S(F

K

)

S( ˆ

F )

Z

ˆ

F

ˆ

g(ˆ

x).

(38)

Using (37)withthe relations (32), (34)and (38)we obtain

kE(φ)k

2

0,K

=

V (K)

V ( ˆ

K)

k ˆ

E( ˆ

φ)k

2

0, ˆ

K

=

V (K)

V ( ˆ

K)

1

√

3

F

2

(p)

p

2

k ˆ

φk

2

0, ˆ

F

=

V (K)

V ( ˆ

K)

1

√

3

S( ˆ

F )

S(F

K

)

F

2

(p)

p

2

kφk

2

0,F

K

=

S( ˆ

F )

V ( ˆ

K)

V (K)

S(F

K

)

1

√

3

F

2

(p)

p

2

kφk

2

0,F

K

.

(39)

Onthereferen etetrahedron

ˆ

K

we extended

ˆ

φ

fromthe fa e

ˆ

F

with

S( ˆ

F ) = 2

√

3

and we have

V ( ˆ

K) =

4

₃

,therefore(39) redu es to

kE(φ)k

2

0,K

=

3

2

V (K)

S(F

K

)

F

2

_(p)

p

2

kφk

2

0,F

K

.

(40)

Step 3  The inequality for the jump term. Using theestimate in (40) , the denition of

R

F

in (7)withthe fa tthat

E([[v]]

T

)

is ontinuous on

∂K

weobtain

k [[v]]

T

k

2

0,F

K

=

Z

F

[[v]]

_T

_{· E([[v]]}

_T

) =

Z

K

R

F

([[v]]

_T

) · E([[v]]

_T

)

≤ kR

F

([[v]]

T

)k

0,K

 3

2

V (K)

S(F

K

)

F

2

(p)

p

2



1

2

k [[v]]

T

k

0,F

K

,

Remark: Sin e

K

isanarbitraryelementadja ent to

F

K

,we anrewritetheestimateinLemma 1as

2

3

M

F

p

2

F

2

_(p)

k [[v]]

T

k

2

0,F

≤ kR

F

[[v]]

_T

k

2

0,K

.

(41) Inthefollowinglemma, we will make useofthe inverse tra einequality onanarbitraryfa e

F

oftheelement

K

kwk

2

0,F

≤

(p + 1)(p + 3)

3

S(F )

V (K)

kwk

2

0,K

(42) in

Σ

p

h

,whi h isproved inTheorem 4in[30 ℄. Lemma 2 For every fa e

F ∈ F

h

and every

v

∈ Σ

p

h

we have the inequality

kR

F

([[v]]

T

)k

0

≤

r

M

F

(p + 1)(p + 3)

6

k [[v]]

T

k

0,F

.

(43)

Proof. The denition of the

[L

2

(Ω)]

3

norm and the tra e inequality in (42) give that for an

arbitrary

v

∈ Σ

p

h

kR

F

([[v]]

_T

)k

0

= sup

w

∈Σ

p

_h

R

Ω

R

F

([[v]]

T

) · w

kwk

0

= sup

w

∈Σ

p

_h

R

F

[[v]]

T

· {{w}}

kwk

0

≤ sup

w

∈Σ

p

_h

k [[v]]

T

k

0,F



R

F



_w

_|

∂KL

+

w

|

∂KR

2



2



1

2

kwk

0

≤ sup

w

∈Σ

p

_h

k [[v]]

T

k

0,F

1

2

(kwk

2

∂K

L

+ kwk

2

_∂K

R

)

1

2

kwk

0

≤ sup

w

∈Σ

p

_h

q

M

F

(p+1)(p+3)

₃

k [[v]]

T

k

0,F

kwk

0

·

 1

₂

 V (K

L

₎

S(F )

3

(p + 1)(p + 3)

kwk

2

∂K

L

+

V (K

R

)

S(F )

3

(p + 1)(p + 3)

kwk

2

∂K

R



1

2

≤ sup

w

∈Σ

p

_h

q

M

F

(p+1)(p+3)

₃

k [[v]]

_T

k

0,F



1

2



kwk

2

_0,K

L

+ kwk

2

_0,K

R



1

₂

kwk

0

≤ sup

w

∈Σ

p

_h

q

M

F

(p+1)(p+3)

₆

k [[v]]

T

k

0,F

kwk

0

kwk

0

=

r

M

F

(p + 1)(p + 3)

6

k [[v]]

T

k

0,F

,

asstated.



4.2 Gårding inequalities and ontinuity estimates

We beginbyproving the Gårdinginequality forthebilinear form oftheBrezzi typeDG

formu-lation(25) .

Lemma 3 There exist onstants

{η

F,0

}

F

∈F

h

, independent of the dis retisation parameter

h =

max

K∈T

h

diam

K

andthe wavenumber

k

, su hthat for all

v

∈ Σ

p

h

and all parameters

η

F

≥ η

F,0

we have the following inequality

β

2

_(k∇

h

× vk

2

0

+ k

h

−

1

2

[[v]]

T

k

2

0,F

h

) − k

2

kvk

2

0

.

Therefore,using (23)itis su ient to 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

).

(45) These ondterm on theleft handside an be estimatedwithanypositive

C

K

L

and

C

K

R

as,

2

Z

F

h

[[v]]

_T

_{· {{∇}

h

× v}} dA

=

X

F∈F

h

Z

F

2

p1 − β

2

h

−

1

2

F

C

−1

K

L

[[v]]

_T

· C

K

L

p1 − β

2

2

h

1

2

F

∇

h

× v

L

|

F

+

2

p1 − β

2

h

−

1

2

F

C

K

−1

R

[[v]]

_T

· C

K

R

p1 − β

2

2

h

1

2

F

∇

h

× v

R

|

F

dA

≤

1

1 − β

2

X

F

∈F

h

h

−1

_F

C

_K

−2

L

k [[v]]

T

k

2

0,F

+

1 − β

2

4

X

F

∈F

h

h

F

C

_K

2

L

k∇

h

× v

L

k

2

0,F

+

1

1 − β

2

X

F∈F

h

h

−1

_F

C

_K

−2

R

k [[v]]

_T

k

2

0,F

+

1 − β

2

4

X

F

∈F

h

h

F

C

_K

2

R

k∇

h

× v

R

k

2

0,F

.

(46)

Applying(42) to the urltermson theright-hand side of (46), we obtain

1 − β

2

4

h

F

C

2

K

L

k∇

h

× v

L

k

2

0,F

≤

1 − β

2

4

h

F

C

2

K

L

(p + 1)(p + 3)

3

S(F )

V (K

L

₎

k∇

h

× v

L

k

2

0,K

L

,

(47)

andinthesame way

1 − β

2

4

h

F

C

2

K

R

k∇

h

× v

R

k

2

0,F

≤

1 − β

2

4

h

F

C

2

K

R

(p + 1)(p + 3)

3

S(F )

V (K

R

₎

k∇

h

× v

R

k

2

0,K

R

.

(48)

For thejumpterms, using(26) , we obtain

C

_K

−2

L

h

−1

F

k [[v]]

T

k

2

0,F

≤ C

K

−2

L

h

−1

F

3

2

V (K

L

₎

S(F )

F

2

(p)

p

2

kR

F

([[v]]

T

)k

2

0

,

(49)

andinthesame way

C

_K

−2

R

h

−1

F

k [[v]]

T

k

2

0,F

≤ C

K

−2

R

h

−1

F

3

2

V (K

R

)

S(F )

F

2

(p)

p

2

kR

F

([[v]]

T

)k

2

0

.

(50) Choosing

C

_K

L

=

s

3 · V (K

L

₎

h

F

(p + 1)(p + 3) · S(F )

and

C

K

R

=

s

3 · V (K

R

₎

h

F

(p + 1)(p + 3) · S(F )

tetrahedra)givesthat

1 − β

2

4

X

F

∈F

h

h

F

C

_K

2

L

k∇

h

× v

L

k

2

0,F

+

1 − β

2

4

X

F

∈F

h

h

F

C

_K

2

R

k∇

h

× v

R

k

2

0,F

≤ (1 − β

2

)k∇

h

× vk

2

0

(51) andsimilarly,summation of (49)and (50)givesthat

1

1 − β

2

X

F

∈F

h

C

_K

−2

L

h

−1

F

k [[v]]

T

k

2

0,F

+

1

1 − β

2

X

F

∈F

h

C

_K

−2

R

h

−1

F

k [[v]]

T

k

2

0,F

≤

1

1 − β

2

F

2

(p)(p + 1)(p + 3)

p

2

kR([[v]]

T

)k

2

0

.

(52) Usingestimates (51)and (52)in(46) weobtain that

2

Z

F

h

[[v]]

_T

_{· {{∇}

h

× v}} dA

≤

1

1 − β

2

F

2

(p)(p + 1)(p + 3)

p

2

kR([[v]]

T

)k

2

0

+ (1 − β

2

)k∇

h

× vk

2

0

.

(53)

Therefore,usingalso (41)we an estimatethe left hand sideof (45)as

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

+

X

F

∈F

h



n

f

+ η

F

−

1

1 − β

2

F

2

_{(p)(p + 1)(p + 3)}

p

2



kR

F

([[v]]

T

)k

2

0

≥ β

2

k∇

h

× vk

2

0

+

X

F∈F

h

h

F



n

f

+ η

F

−

1

1 − β

2

F

2

(p)(p + 1)(p + 3)

p

2



·

p

2

F

2

_(p)

2

3

M

F

h

−1

F

k [[v]]

T

k

2

0

.

(54) Therefore,we have to hoose

η

F

su h that

h

F

M

F

·



n

f

+ η

F

−

1

1 − β

2

F

2

_{(p)(p + 1)(p + 3)}

p

2

 2

3

p

2

F

2

_(p)

≥ β

2

(55)

andwiththis (45)is satised.



Remarks:

1. Given that

n

f

= 4

for tetrahedra we anmake the onditionfor

η

F

expli it,

η

F,0

=

F

2

_(p)

p

2



3β

2

2h

F

M

F

+

(p + 1)(p + 3)

1 − β

2



− 4.

(56)

2. The oer ivity onstant

β

is, however, still undened. Using the a priori error analysis, whi h will be dis ussed inthe nextse tion, we an nd anoptimal valuefor

η

F,0

.

3. A straightforward estimationgivesthat

F

2

(p)(p+1)(p+3)

p

2

≥ 1

,whi htogether with(55)gives that

n

f

+ η

F

≥ 1

if

0 ≤ β

2

_{< 1.}

(57)

Observe thatfor anarbitrary

K

we havediam

K = h

F

≥ m

F

,where

F

is afa eof

K

and

m

F

isthe height orrespondingto

F

. Hen e,

S(F )h

F

≥ S(F )m

F

= 3V (K)

and therefore,

max

F

∈F

h

h

F

M

F

≥ max

F

∈F

h

h

F

max



S(F )

V (K

L

₎

,

S(F )

V (K

R

₎



≥ 3.

(58)

Using the method inLemma 3 we analso obtain aboundfor thepenalty parameterinthe

interiorpenalty(IP)method(18) su h thatthe Gårdinginequalityisvalid.

Lemma 4 There exist onstants a

F,0

, independent of the dis retisation parameter

h =

max

K∈T

h

diam

K

andthewave number

k

,su h that for all

v

∈ Σ

p

h

andall parameters a

F

≥

a

F,0

we have the following inequality

B

ip

h

(v, v) ≥ β

2

kvk

2

DG

− (k

2

+ β

2

)kvk

2

0

.

(59)

Proof. A ording to theproof ofLemma 3 itis su ient to prove that

k∇

h

× vk

2

0

− 2

Z

F

h

[[v]]

_T

_{· {{∇}

h

× v}} dA +

X

F

∈F

h

a

F

k [[v]]

T

k

2

0,F

≥ β

2

(k∇

h

× vk

2

0

+ k

h

−

1

2

[[v]]

T

k

2

0,F

h

).

(60) Withthesame hoi e of oe ients

C

K

L

and

C

K

R

asinLemma 3,weobtain theinequality

2

Z

F

h

[[v]]

_T

_{· {{∇}

h

× v}} dA

≤

1

1 − β

2

X

F

∈F

h

C

_K

−2

L

h

−1

_F

k [[v]]

_T

k

2

0,F

+

1

1 − β

2

X

F

∈F

h

C

_K

−2

R

h

−1

_F

k [[v]]

_T

k

2

0,F

≤

1

1 − β

2

(p + 1)(p + 3)

3



S(F )

V (K

L

₎

+

S(F )

V (K

R

₎



k [[v]]

T

k

2

F

h

,0

.

Substituting (46) into the right-hand side and using also (51) , the left hand side of (60) is

estimatedas

k∇

h

× vk

2

0

− 2

Z

F

h

[[v]]

_T

_{· {{∇}

h

× v}} dA +

X

F

∈F

h

a

F

k [[v]]

T

k

2

0

≥ β

2

k∇

h

× vk

2

0

(61)

+

X

F

∈F

h

h

F



a

F

−

1

1 − β

2

(p + 1)(p + 3)

3



S(F )

V (K

L

₎

+

S(F )

V (K

R

₎



h

−1

_F

_{k [[v]]}

_T

_k

2

_F

_h

_,0

.

We have to hoose thenthe parametera

F

onthe fa e

F

su hthat

h

F



a

F

−

1

1 − β

2

1

3

(p + 1)(p + 3)



S(F )

V (K

L

₎

+

S(F )

V (K

R

₎



≥ β

2

,

a

F,0

≥

β

2

h

F

+

1

1 − β

2

1

3

(p + 1)(p + 3)



S(F )

V (K

L

₎

+

S(F )

V (K

R

₎



.

(62)

Thisprovesthe lemma.



Intheerror analysisone hasto onsider(see [18 ℄)theextended( f. (23))bilinear form

B

br

: (H

0

(

url

, Ω) + Σ

p

h

) × (H

0

(

url

, Ω) + Σ

p

h

) → R,

whi h isgiven as

B

br

(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

))

Ω

andthelinearform

J

h

: H

0

(

url

, Ω) + Σ

p

h

→ R

,dened as

J

h

(v) = (J , v)

Ω

whenzeroboundary onditions are onsidered. Infollowing twolemmas we usethenotation

M = max

F

∈F

h

r

h

F

M

F

(p + 1)(p + 3)

6

.

Using (58)for

p ≥ 1

wehave that

M ≥ 2

.

Using the inverse tra einequality (42)we also havethat

k∇

h

× u

L

k

2

0,F

≤

(p + 1)(p + 3)

3

S(F )

V (K

L

₎

k∇

h

× uk

2

0,K

L

≤ h

−1

F

_F

max

_∈F

h

M

F

h

F

(p + 1)(p + 3)

3

k∇

h

× uk

2

0,K

L

≤ 2h

−1

F

M

2

k∇

h

× uk

2

0,K

L

(63)

anda similarestimateholds for theneighboringelement

K

R

.

Lemma 5 The bilinearform

B

br

is ontinuous on

(H

0

(

url

, Ω) + Σ

p

h

) × (H

0

(

url

, Ω) + Σ

p

h

)

with respe t to the DGnorm,i.e. the followinginequalityholds 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

br

(u, v) ≤ Ckuk

DG

kvk

DG

,

(64) where

C = max

F

∈F

h



k

2

,

5

4

M

2

_(n

f

+ η

F

)



.

Proof. Using the triangleinequality,Lemma2,theresultoftheeigenvalue problemdis ussedin

theAppendix, theestimate

M ≥ 2

and (57)we obtainthat

B

br

(u, v) ≤ |(∇

h

× u, ∇

h

× v)

Ω

| + k

2

|(u, v)

Ω

| +

X

F

∈F

h