Notes on solving Maxwell equations, part I, finite elements method using vector elements

Notes on solving Maxwell equations, part I, finite elements

method using vector elements

Citation for published version (APA):

Rook, R. (2011). Notes on solving Maxwell equations, part I, finite elements method using vector elements. (CASA-report; Vol. 1110). Technische Universiteit Eindhoven.

Document status and date: Published: 01/01/2011

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

EINDHOVEN UNIVERSITY OF TECHNOLOGY

Department of Mathematics and Computer Science

CASA-Report 11-10

February 2011

Notes on solving Maxwell equations Part 1:

Finite elements method using vector elements

by

R. Rook

Centre for Analysis, Scientific computing and Applications

Department of Mathematics and Computer Science

Eindhoven University of Technology

P.O. Box 513

5600 MB Eindhoven, The Netherlands

ISSN: 0926-4507

Notes on solving Maxwell equations Part 1:

Finite Elements Method using vector elements

R.Rook

February 3, 2011

1 Introduction

The content of this document is a brief description of the steps taken in the process of designing and writing a simulation program for solving the Maxwell equations using the nite element method (FEM). It covers the theory from the literature by various authors. It deals with the conversion to certain conventions, like the orientation of axis and the time convention for time-harmonic solutions and the non-dimensionalization process. The mathematical foundation is far from complete, but the reader is referred to the literature for the details. The application is the scattering by an object of electromagnetic waves. Part 2 will cover the far-eld calculation given the near-eld solution obtained by the FEM.

2 Problem

A typical problem is an incident plane wave refracting and diracting by a 2D-periodic object (e.g. a resist line in air on top of a layered substrate, see gure 1). The Maxwell equations with the proper boundary conditions are the governing equations.

The materials are assumed to be isotropic and homogeneous within the materials (constant permittivity and permeability).

2.1 Maxwell Equations

The time convention used in this report is exp(−iωt) (as used by the authors

Chew, Monk), where angular frequency ω [s−1_{], ω > 0. This particular choice}

n2 n0 K θ x z y M L n₁ n₃ Figure 1: 2D-grating on Ω

corresponding direction of the axis, e.g. exp(ikx−iωt), with k > 0, travel in the positive x-direction. Note that converting a solution, say E, to the exp(iωt) time convention (as used by the author Zhu & Cangellaris) is simply the conjugated solution ¯E.

2.1.1 Dimensional

On a domain Ω ∈ R3_{, the time-harmonic electric eld E [Vm}−1_{], the}

time-harmonic magnetic eld H [Am−1_{], electric displacement D [Cm}−2_{] and}

mag-netic induction B [T] for isotropic linear media are described by (Faraday's law) ∇ × E = iωB − M (Ampère's law) ∇ × H = J − iωD, (Ohm's law) J = Jc+ Jv,

(Gauss' electric eld law)

∇ · D = 1

(Gauss' magnetic eld law)

∇ · B = 1

iω∇ · M,

with the complex electric permittivities  [Fm−1_{] and permeability µ [Hm}−1_].

The electric current density Jv[Am−2_{] represents the impressed current sources}

in the domain (imposed electric eld), Jc [Am−2_{] accounts for induced current}

ow in the medium due to the rpecence pf conductor and/or dielectric loss and

M [Ts−1_{] is the non-physical magnetic current density source.}

The constitutive relations between the electromagnetic eld quantities and theit associated ux densities for isotropic materials are given by

D = E, B = µH, Jc = γE,

with γ [Sm−1_{] is the (scalar) conductivity.}

Then, in each isotropic material, the vector wave Maxwell equations can be derived as

∇ × (∇ × E) − k2_{E = iωµJ − ∇ × M,}

∇ × (∇ × H) − k2H = iωM + ∇ × J,

with k2 _{= ω}2_{µ and shows the symmetry between the electric and magnetic}

elds (the duality principle; Chew, page 9).

If the eective permittivity is ˜ =  − iγ/ω, a material with index i has the

(eective) relative refraction index n2

i = µ˜/µ00 = µr,ir,i, where µr,i [-] and

r,i [-] are the relative permittivity and permeability, respectively. In air at

low eld strengths (same as in vacuum γ = 0,  = 0, µ = µ0), the reference

wave number k0 [m−1_{] is dened as k}2

0 = ω2µ00 = ω2/c2, with c the speed of light in vacuum. Combining these with the relation for the wave length in air

λ0= 2πn0/k0 [m], we can derive an expression for the frequency squared wave

number for the i-th material k2

i = ω

2_{µ˜ = 4π}2_n2 i/λ

2

0 or ki= k0ni/n0.

From these relations (and zero magnetic current density M), we derive the time-harmonic vector wave Maxwell equation (Zhu & Cangellaris, page 15; Monk, page 7)

2.1.2 Non-dimensional

Spatial variables x, y and z are non-dimensionalized with respect to wave length

λ0 as the used wavelengths are comparable to the dimensions of the scatterer.

The domain is reduced to a rectangular domain Ω = [0, K]×[0, L]×[0, M], with nite K, L and M [m] in each direction.

The unknowns are the non-dimensional electric vector eld E(x, y, z) and the

magnetic vector eld H(x, y, z), where E0is a reference value, e.g. the amplitude

of the incident eld. Take H0=q0

µ0E0as reference value (intrinsic impedance

in free space) for the magnetic eld and Jv0 = H0/λ0 for the current density

source and M0= E0/λ0for the magnetic current density.

After non-dimensionalization, the eld relations are

H = 1 2πin0µr,i(∇ × E + M ), (1) and E = 1 2πin0r,i(Jv− ∇ × H), (2) where we used ω˜q_µ 0 0λ0= ωr,i0 q_µ 0 0λ0= ωr,i √ 0µ0λ0= r,ik0λ0= 2πn0r,i.

The electric eld E will satisfy the following (non-dimensional) Maxwell

equa-tion (from now on the refracequa-tion index for air n0= 1)

∇ × (∇ × E) − 4π2_n2

iE = 2πiµr,iJv.

Instead of the full eld E we can rewrite the problem for the scattered eld Escat= E − Einc,

with Einc _{a solution in air (i = 0). So, the scattered eld satises (no current}

sources, Jv= 0) ∇ × (∇ × Escat_{) − 4π}2_n2 iE scat_{= 4π}2_(n2 i − n 2 0)E inc_,

which is referred to as the scattered eld formulation, where materials (other than air) will act as sources.

Another description is the contrast formulation for the grating problem Econtr= E − Emulti,

y x φ p d ψ θ z

Figure 2: Conical incident wave

where Einc _{is replaced by the exact solution of the multiple layer stack E}multi

(thus without the grating). In this case Econtr _satises

∇ × (∇ × Econtr) − 4π2n2iE contr = 4π2(n2i − n 2 res)E multi ,

where nres= n0 in and nres= ni outside the scatterer, which makes the

right-hand side is non-zero in the resistor only.

2.1.3 Scalar wave equations and polarization

In an isotropic medium where the electromagnetic properties are varying in one direction, say the z-direction, the Maxwell equations can be reduced to two scalar equations that are decoupled. They characterize two types of waves:

the transverse-electric (TEz or TE) and the transverse-magnetic (TMz or TM)

waves. Assuming that the eld is linearly polarized, the electric eld in the TE waves are directed in the x, y-plane and can be rotated around the z-axis to point only in the y-direction E = Ey. A better way of characterizing this

wave is to use the z-component of the magnetic eld Hz. If Hz 6= 0, the eld

describe TE waves. Similarly, by duality, the Ez component characterizes TM

waves (Chew, page 45). 2.1.4 Incident wave

The general incident (polarized) wave in air is

Einc= p exp(2πin0d · x), p · d = 0. For a conical incident wave, the vector p is the unit vector

p = 



cos(ψ) cos(θ) cos(φ) − sin(ψ) sin(φ) cos(ψ) cos(θ) sin(φ) + sin(ψ) cos(φ)

cos(ψ) sin(θ)



where ψ is the polarization angle (π/2 for TE and 0 for TM), θ is the polar angle and φ is the azimuth of the incident wave. It is perpendicular to the wave direction d d =   sin(θ) cos(φ) sin(θ) sin(φ) − cos(θ)  .

2.2 Boundary conditions

In general, two types of boundary conditions are specied: Robin type Silver-Müller condition (Monk, page 11) and (essential) Dirichlet boundary conditions or PEC boundaries (Monk, page 9; Zhu & Cangellaris, page 6).

2.2.1 Silver-Müller

On the outside boundary we may pose the Silver-Müller condition as a Robin type condition

n × (∇ × E) − σET = −r, (3)

with

r = −n × (∇ × Einc) + σE_Tinc,

and where ET = n × (E × n), the tangential projection of E on the

bound-ary surface. Note that in the scattered eld or contrast eld formulation, the right-hand side vanishes r = 0 and the Silver-Müller condition reduces to a homogeneous condition.

2.2.2 Dirichlet

When the solution is known, we may pose Dirichlet or essential boundary con-ditions

E = Einc(x).

Note that the tangential part of the incident wave only contributes to the prob-lem.

2.2.3 Perfectly conducting

At perfect electrically conduction (PEC) boundaries n × E = 0 holds, with n the normal direction of the boundary.

2.2.4 Periodic boundaries

The quasi-periodic boundary conditions are dened as E(K, y, z) = E(0, y, z) exp(2πin0Kd1),

E(x, L, z) = E(x, 0, z) exp(2πin0Ld2),

E(x, y, M ) = E(x, y, 0) exp(2πin0M d3) in the x-, y- and z-direction, respectively.

3 Finite Element Method

As the geometry involves arbitrary shapes, the problem is solved using FEM with vector basis functions. The FEM allows us to create a mesh consisting of triangular elements that perfectly overlays the geometry. The vector functions solve a number of problems surfacing when using standard node-based scalar base functions: the problem of spurious modes. The rst is the enforcement of the tangential electric and magnetic eld continuity at the material interfaces. Also the continuity of the normal components of the electric and magnetic ux densities cannot be guaranteed. The second is the inability to model properly the null space of the curl operator (the null space consist of all vector elds

F = −∇Ψ, where Ψ is an arbitrary scalar eld) (Zhu & Cangellaris, page 15).

3.1 Variational formulation

In the Galerkin method, this solution satises the variational formulation Z Ω (∇ × E) · (∇ × φ) − 4π2n2_iEφdV + I n × (∇ × E)φTdΓ = Z Ω F φdV or when the Robin boundary condition (3) is included

Z Ω (∇ × E) · (∇ × φ) − 4π2n2_iEφdV + I σETφTdΓ = Z Ω F φdV + I rφTdΓ, (4)

for all φ ∈ H(curl; Ω).

The Dirichlet boundary conditions have the following variational form I

ETφTdΓ =

I

Note that within the variational formulation (4), if you set σ and r to a large

value, say 1010_{, you obtain a the same result (within a certain accuracy). This is}

also known as the penalty method, but this will destroy a well-conditioned mass matrix. Rather eliminate the Dirichlet boundary condition from the system of equations.

Assume that the perfectly conducting boundary conditions are located on Γ and the Silver-Müller conditions on Σ. In order for all the integrals to be well-dened, we should use the space X dened by

X =u ∈ H(curl; Ω) | n × u = 0 on Γ ∧ uT ∈ (L2_(Σ))3_{on Σ}

with the Sobolev space

H(curl; Ω) = v ∈ (L2(Ω))3| ∇ × v ∈ (L2_(Ω))3 _.

However, note that X is not ideal to dene the nite elements space on (Monk, page 99) as it contains the null-space for the curl operator.

3.2 Ane maps

A nite element is dened by the triple (K, PK, ΣK), where K is a geometric

domain, PK is the space of functions on K and ΣK a set of functionals on PK

also called the degrees of freedom (DOFs) of the nite element. An important step is to map a these elements to a reference element ( ˆK, P_Kˆ, Σ_Kˆ).

The solution E is approximated and will be the sum of the shape functions φi:

E =XLi(E)φi,

where Li are the DOFs in the system.

Suppose the domain is subdivided into elements K. We can dene the shape functions on these elements. These shape functions are mapped with an

iso-parametric mapping FK : ˆx → x from reference element ˆK to the physical

element K. For tetrahedra we take the ane mapping of the 4 vertices. The functions and their derivatives in the variational formulation are then trans-formed in the following way

(∇ × φ) ◦ FK = 1

det(dFK)dFK

ˆ ∇ × ˆφ

φ ◦ FK= (dFK)−Tφˆ

The (ane) mapping FK is constant on K and we have ˆE =P_{Li( ˆ}ˆ _{E) ˆ}φi, thus

ˆ

Li = Li. This means that the sets of DOFs on K are invariant under

trans-formation (see Monk, Lemma 5.34, page 131). The vector elements described below are formulated on these reference elements.

v0 e0 v1 e1 z y x ^ ^ e5 ^ s1 s0 s3 e3 v3 v2 e2 e4 s2

(a) Local numbering vertices and DOFs

(b) Basis function

3.3 Curl-conforming elements

Curl-conforming means that the nite element space that is a subspace of

H(curl; Ω). The rst family of curl-conforming elements (edge elements) has

shape functions that span the polynomial space Rk = (Pk−1)3 ⊕ Sk, with

Sk = {p ∈ ( ˜Pk)3 | x · p = 0}, Pk the space of polynomial functions with

or-der k ≥ 1 and ˜Pk the set of homogeneous polymials with exact order k (Monk,

page 108). A second family of curl-conforming elements can be obtained, where

the space of functions is larger PK = (Pk)3, e.g. hp-elements (Monk, page 202;

τ

Figure 3: Denition edge tangential 3.3.1 Whitney element

The Whitney element, or linear edge element is the lowest order curl-conforming element and resolves space R1.

The degree of freedom Lj is dened on the j-th edge in the mesh, and is the

integral of the inner product of the solution and the tangential vector τej dening

that edge: Lj(E) = 1 length(ej) Z ej τe_j· E ds.

In Whitney elements, the shape functions are constructed from the barycentric coordinates

λ0= 1 − ˆx − ˆy − ˆz, λ1= ˆx, λ2= ˆy, λ3= ˆz, and the shape functions are derived by

ˆ

φe{i,j}= λi∇λjˆ − λj∇λiˆ ,

where e {i, j} is the edge with vertices i and j (i, j = 0, ..., 3), on which the shape function is dened.

This results in the following set of k = 6 shape functions: ˆ φ0= ˆφe{0,1} = [1 − ˆy − ˆz, ˆx, ˆx] T ˆ φ1= ˆφe{1,2} = [−ˆy, ˆx, 0] T ˆ φ2= ˆφe{0,2} = [ˆy, 1 − ˆx − ˆz, ˆy] T ˆ φ3= ˆφe{0,3} = [ˆz, ˆz, 1 − ˆx − ˆy] T ˆ φ4= ˆφe{1,3} = [−ˆz, 0, ˆx] T ˆ φ5= ˆφe{2,3} = [0, −ˆz, ˆy] T

v0

e0

v1

∆

λ

₁

∆

λ

₀ 0

λ = 0

λ = 0

1

Figure 4: Barycentric coordinates and their curls

ˆ ∇ × ˆφ0= [0, −2, 2] T , ˆ∇ × ˆφ1= [0, 0, 2] T ˆ ∇ × ˆφ2= [2, 0, −2] T , ˆ∇ × ˆφ3= [−2, 2, 0] T ˆ ∇ × ˆφ4= [0, −2, 0]T, ˆ∇ × ˆφ5= [2, 0, 0]T

These functions are dened in the reference space ˆK.

3.3.2 Quadratic edge elements

R2is spanned by 20 polynomial vectors. The rst 12 spans the space (P1)3and

are the Whitney shape functions with 6 additional linear shape functions which related to the DOFs

L0_j(E) = 1

length(ej) Z

ej

E · τe_j(3 − 6s) ds.

The other 8 shape functions are related to the faces and their DOFs are

Lj(E) = 1 area(fj) Z fj E · (−8τ1j+ 4τ2j) dA and L0_j(E) = 1 area(fj) Z fj E · (4τ1j+ 4τ2j) dA

τ

τ

1

2

Figure 5: Denition face tangentials

where τ1jand τ2j span the 2 sides of the triangular face. Note that area(fj) =

1

2kτ1j× τ2jk.

The shape functions are the 6 Whitney shape functions from the previous section and 6 additional linear edge functions

ˆ φ0_e{i,j}= λi∇λˆ j+ λj∇λˆ i, or ˆ φ6= ˆφ0e{0,1}= [1 − 2ˆx − ˆy − ˆz, −ˆx, −ˆx] T ˆ φ7= ˆφ0e{1,2}= [ˆy, ˆx, 0] T ˆ φ8= ˆφ0_e{0,2}= [−ˆy, 1 − ˆx − 2ˆy − ˆz, −ˆy]T ˆ φ9= ˆφ0_e{0,3}= [−ˆz, −ˆz, 1 − ˆx − ˆy − 2ˆz]T ˆ φ10= ˆφ0_e{1,3}= [ˆz, 0, ˆx]T ˆ φ11= ˆφ0_e{2,3}= [0, ˆz, ˆy]T These result in 2 functions per edge.

The other 8 functions spans S2 and have the form

ˆ

φs{i,j,k}= λiλj∇λkˆ − λiλk∇λjˆ and

ˆ

φ0_s{i,j,k}= λiλj∇λkˆ − λjλk∇λiˆ These are related to the faces (also 2 functions per face):

ˆ φ12= ˆφs{0,1,2}= [−(1 − ˆx − ˆy − ˆz)ˆy, (1 − ˆx − ˆy − ˆz)ˆx, 0] T ˆ φ13= ˆφ0_s{0,1,2}= [ˆxˆy, (1 − ˆx − ˆz)ˆx, ˆxˆy]T ˆ φ14= ˆφs{0,1,3}= [−(1 − ˆx − ˆy − ˆz)ˆz, 0, (1 − ˆx − ˆy − ˆz)ˆx] T ˆ φ15= ˆφ0_s{0,1,3}= [ˆxˆz, ˆxˆz, (1 − ˆx − ˆy)ˆx]T ˆ φ16= ˆφs{1,2,3}= [0, −ˆxˆz, ˆxˆy] T ˆ φ17= ˆφ0s{1,2,3}= [−ˆy ˆz, 0, ˆxˆy] T ˆ φ18= ˆφs{0,2,3}= [0, −(1 − ˆx − ˆy − ˆz)ˆz, (1 − ˆx − ˆy − ˆz)ˆy] T ˆ φ19= ˆφ0s{0,2,3}= [ˆy ˆz, ˆy ˆz, (1 − ˆx − ˆy)ˆy] T

In general, the shape functions φj do not form a set of functions with the

δ-property, with respect to the linear forms Li(DOFs). However, we can construct

one. Suppose we have our basis {φj}and we are looking for another basis {θj}

such that θj = N X k=1 akjφk. So, Li(θj) = Li N X k=1 akjφk ! = N X k=1 akjLi(φk) = δij

which reduces to LA = I, with L = {Li(φk)}N

i,k=1. For each element p in the polynomial space we have

p =X j Lj(p)θj=X j Lj(p)X k akjφk=X k X j akjLj(p)φk=X k ck(p)φk with ck(p) = X j akjLj(p),

the DOFs corresponding to the basis {φk}.

The resulting matrix A for the second order Nedelec nite element is A =  I 0 −B I  , B = {Li(φj)}j=0,11_i=12,19 and explicitly, B =             −4 4 0 0 0 0 −11_{/3 −1}1_{/3 2}2_{/3 0 0 0} 4 0 4 0 0 0 −11_{/3 2}2_{/3 −1}1_{/3 0 0 0} −4 0 0 0 4 0 −11_{/3 0 0 2}2_{/3 −1}1_{/3 0} 4 0 0 4 0 0 −11_{/3 0 0 −1}1_{/3 2}2_{/3 0} 0 −4 0 0 0 4 0 −11_{/3 0 0 2}2_{/3 −1}1_/3 0 4 0 0 4 0 0 −11_{/3 0 0 −1}1_{/3 2}2_/3 0 0 −4 0 0 4 0 0 −11_{/3 2}2_{/3 0 −1}1_/3 0 0 4 4 0 0 0 0 −11_{/3 −1}1_{/3 0 2}2_/3            

For example, shape function θ0 = φ0+ 4φ12− 4φ13+ 4φ14− 4φ15 owns the

δ-property Li(θ0) = δi,0 and DOF c12(p) = 4L0(p) − 4L1(p) + 11_{/3L6(p) +}

11_{/3L7(p) − 2}2_{/3L8(p) + L12(p)}corresponds to shape function φ12.

Another way of constructing a new basis is to integrate the DOFs using exact in-tegration, a reference basis for the polynomial space and then use the δ-property to obtain the linear combination of the reference basis (see P. Monk, page 140). The δ-property is used to eliminate the DOFs corresponding to the essential boundary conditions from the system.

3.4 PML

The PML takes care of damping incoming (reected) waves. The coordinates perpendicular to the boundary is stretched by a complex number

˜

x = x + i

Z x

0

σ (ξ) dξ, d˜x = (1 + iσ)dx = γdx,

where the real functions σ is dened as σ(ξ) = ( σ∗ξ−x0 x1−x0 2 x0< ξ < x1 0 ξ < x0 ,

for some real constant σ∗_{, an upper bound x0and a lower bound x1. Generally,}

we have mapping G

G : x → ˜x.

3.4.1 Maxwell equations The electric eld transforms as

˜

E ◦ G = dG−TE.

The idea is that E is continued into the upper half of the complex plane in each direction, so that

˜

∇ × ( ˜∇ × ˜E) − 4π2n2_iE = ˜˜ F . If we change back to the physical coordinates we derive that

( ˜∇ × ˜u) ◦ G = A∇ × B ˜u,

A =   1 γ2γ3 0 0 0 _γ1 1γ3 0 0 0 _γ1 1γ2  , B =   γ1 0 0 0 γ2 0 0 0 γ3  .

From the new variable ¯E = B ˜E and right-hand side ¯F = B ˜F we obtain

∇ × (µ−1∇ × ¯E) − 4π2n2_iµ ¯E = µ ¯F , with µ = (BA)−1=   γ2γ3 γ1 0 0 0 γ1γ3 γ2 0 0 0 γ1γ2 γ3  .

Note that ¯E matches E at the interface of the PML because γ3(z) = 1 at the

interface.

3.4.2 Variational formulation

The variational formulation in the PML region will be (Monk, page 378) Z Ω (µ−1∇ × ¯E) · (∇ × φ) − 4π2n2_iµ ¯EφdV + I σ ¯ETφTdΓ = I ¯ rφTdΓ + Z Ω µ ¯F φdV. 3.4.3 Boundary conditions

The corresponding boundary condition in the PML is (Monk, page 381) n × (µ−1∇ × ¯E) = σ ¯ET − ¯r,

where ¯r should be derived from the known (stretched) incident wave ¯

3.5 Constraints

Direchlet and periodic boundary conditions can be looked upon as constraints to the DOFs. Suppose the mass matrix is constructed straightforward following the variational formulation, resulting in the linear system

Ax = b,

where x are the DOFs corresponding to the functional space {φi}. The

con-straint equation for the boundary conditions is an ane function x = Cy + d,

where C is a non-square complex-valued matrix and d a vector. DOFs y

cor-respond to functional subspace {Φj} ⊆ {φi} (where the boundary DOFs are

eliminated).

3.5.1 Boundary conditions Recall the expansion formula

E = N X

i=1

φixi, xi= Li(E),

and in terms of the functions in subspace (M ≤ N) you also may write E =

M X

j=1 Φjyj

and the DOFs are related by

xi = Li(E) = Li   M X j=1 Φjyj  = M X j=1 Li(Φj)yj= M X j=1 Cijyj.

Dirichlet boundary condition at location N (spanned by φN) in constraint form

is xN = LN(E) = dand LN(Φj) = 0, which gives CN j = 0, dN = d. A similar

constraint can be found for the periodic boundary conditions. For example,

construct the function Φ1 = φ1 + φN, a linear combination of the original

base functions at the boundaries, and apply the expansion above. This gives C11= L1(Φ1) = L1(φ1) = 1and CN 1= LN(φN) = 1.

3.5.2 Reduced system

The above system is solved by eliminating the constraint function for x and solving the following system for y

CHACy = CHb − CHAd.

3.6 Errors

The error dened in the H(curl; Ω) space is dened as E − Eexact _H( curl;Ω)= E − Eexact 2 (L2_(Ω))3+ ∇ × (E − Eexact) 2 (L2_(Ω))3 =

(x − xexact)HA(x − xexact),

where x is the vector of DOFs and A is a mass matrix of the system with natural boundary conditions n × (∇ × E) = 0.

The relative error is dened as

kE − Eexact_k

H(curl;Ω) kEexact_k

H(curl;Ω) .

Error estimates can be obtained by restriction on the coecients in the problem (Monk, pages 169 and 277).

3.7 Preconditioners

Preconditioners are used to design faster iterative solvers such as algebraic MG methods. The motivation to use Nested Multigrid Potential Preconditioner (NMGAV) preconditioner is the slow convergence of the FEM system. The rst reason is the presence of non-zero gradient components or spurious modes and the resolution of the low-frequency modes (Zhu & Cangellaris, page 160). Two types of linear subspaces are used in this algorithm and need transfer operators to project them to the original space and vice versa.

Basically, dene a linear subspace V of the solution space W and construct a transfer operator G

where W and V is a row vector of the base functions that spans the consecutive spaces.

If an approximated solution xW for system A(WT_{, W )xW = fW is obtained, we}

solve the following preconditioned system

A(VT, V )xv= GT(fW − A(WT, W )xW).

and correct xW as

xw7→ xw+ Gxv.

For the Maxwell equations the two types of subspaces are gradient space ∇Wi _⊂

Wi _{and the coarse grid spaces W}i

kh⊂ W(k−1)hi ⊂ · · · ⊂ W i h. 3.7.1 Potential formulation

In the case of the gradient subspace ∇W , we formulate the potential formulation of the problem. It follows from the Gauss' electric eld law constraining the solution (Zhu & Cangellaris, page 167) (non-dimensional)

∇ · (rE) = 1

2πi∇ · Jv,

with Jv is the electric current density with boundary condition Jv· n = 0.

We construct the weak formulation by multiplying it with the test function

φ ∈ ∇W: Z Ω ∇ · (rE)φ dV = 1 2πi Z Ω ∇ · Jvφ dV

and E ∈ ∇W . Partial integration gives Z Ω ∇ · (rEφ) − rE · ∇φ dV = 1 2πi Z Ω ∇ · (Jvφ) − Jv· ∇φ dV. On the rst part on both sides we can apply divergence theorem

Z Ω ∇ · (rEφ) = I ∂Ω rE · nφ dS

Z Ω ∇ · (Jvφ) = I ∂Ω Jv· nφ dS = 0.

Using the Ampère and Faraday relations, again the boundary conditions Jv·n =

0, φ = 0 and a vector identity, the right-hand side integral can be rewritten as

(Zhu & Cangellaris, page 167, non-dimensional) I ∂Ω rE · nφ dS = 1 2πi I ∂Ω (Jv− ∇ × H) · nφ dS = −1 2πi I ∂Ω n × H · ∇φ dS = −1 2πiµr I ∂Ω n × ∇ × E · ∇φ dS

This multiplied by 2πiµr, we get the divergence free constraint −4π2_n2 i Z Ω E · ∇φ dV + I ∂Ω n × ∇ × E · ∇φ dS = 2πiµr Z Ω Jv· ∇φ dV.

Including the Robin boundary condition we nally have

−4π2n2i Z Ω E · ∇φ dV + I ∂Ω σET · ∇φ dS − I ∂Ω r · ∇φ dS = 2πiµr Z Ω Jv· ∇φ dV.

The PML can easily be included

−4π2_n2 i Z Ω µ ¯E · ∇φ dV + I ∂Ω σ ¯ET· ∇φ dS − I ∂Ω ¯ r · ∇φ dS = Z µ ¯F · ∇φ dV.

From this, the potential formulation is in fact the same equation as the eld formulation (Zhu & Cangellaris, page 181), i.e.

A(∇WT, ∇W ) = GTA(WT, W )G

f∇W(∇WT) = GTfW

3.7.2 Multigrid

The nite element is divided into subelements and the corresponding base func-tions form a linear subspace.

The algorithm (Zhu & Cangellaris, page 170) for NMGAV

xh_W ← N M GAV (fh

E, n = 1),

where the superscript h denotes the nest grid n = 1 and the coarsest grid is N h.

1. if n = N, solve the curl formulation exactly. 2. else

(a) pre-smooth v times on the potential formulation

i. pre-smooth curl formulation using forward GS v times

ii. calculate new right-hand side f∇W ← GT_(fnh

W − Ax

nh

W)

iii. pre-smooth on potential formulation using forward GS v times iv. correct electric eld xnh

W ← x

nh

W + G x∇W

(b) calculate new right-hand side f(n+1)h

W ← Q T_(fnh W − AxnhW) (c) x(n+1)h W ← N M GV A(f (n+1)h W , n + 1) (d) xnh W ← xnhW + Q x (n+1)h W

(e) Post-smooth v times on the potential formulation

i. calculate new right-hand side f∇W ← GT_(fnh

W − Ax

nh

W), x∇W ←

0

ii. post-smooth on potential formulation using backward GS v times

iii. calculate new right-hand side fnh

W = f

nh

W − AG x∇W

iv. post-smooth on curl formulation using backward GS v times

v. correct electric eld xnh

W ← xnhW + G x∇W

Discarding steps 2(a)iv, 2(b), 2(c), 2(d) and 2(e)i in the above algorithm renders the single grid potential preconditioner (SGAV, Zhu & Cangellaris page 168).

4 Implementation

The previous section give the general algorithm. The following subsections cover in more detail the choices made in the implementation stage.

4.1 Boundary conditions

4.1.1 Silver-Müller

For the general incident wave, we assume that zero modes are dominant and approximate the scattered eld at the top boundary to

Escat≈ q exp(−2πin0d · x).

Which means we have to look for a suitable σ such that n × (∇ × Escat_{) =}

σn × (Escat× n). The curl of Escat_{, ∇ × E}scat_{= −2πin0d × E}scat_{, and assume}

almost perpendicular incident waves d ≈ −n, we approximate σ = −2πin0. For example for TE polarization in 2D, the plane wave incident eld

E₂inc(x, y, z) = B0exp(−2πin0z),

E1inc= E3inc= 0

is a solution of the Maxwell equation. At z = M, the electric eld E satises ∂E2

∂z − 2πin0E2= −2πi2n0B0exp(−2πin0z),

E1= E3= 0.

At z = 0, we have

−∂E2

∂z − 2πin3E2= 0,

E1= E3= 0.

For the scattered eld formulation with Einc _{subtracted on the whole domain}

we get at the top

∂Escat 2 ∂z − 2πin0E scat 2 = 0, E1= E3= 0

and at the bottom

−∂E scat 2 ∂z − 2πin0E scat 2 = −2πi(n0− n3)E inc_, E1= E3= 0.

4.1.2 PEC

A PEC boundary can be implemented as a homogeneous Dirichlet boundary condition for the DOFs (in the function space these are the tangential parts of the eld E).

4.2 PML

The PML is basically a material with anisotropic properties, which must be cho-sen carefully. As the problem is described by the Maxwell equations it requires a boundary conditions.

4.2.1 Parameters

The scattered eld inside the PML is multiplied by a function exp  −2πn0 Z z z0 σ3(ξ)dξ  .

The rule of thumb is that Rz

z0σ3(ξ) dξshould be constant for varying PML layer

thickness. The constant σ∗_{cannot be taken too large, because the incident wave}

amplitude increases exponentially towards the boundary of the domain, which requires a ner mesh.

4.2.2 Boundary conditions

A PML is supposed to damp out all waves in a certain direction and if the it performs well, PEC conditions would suce. But the radiation condition is also a good choice. We construct the equation which describes the scattered wave

˜

Escatapproximately at the boundary as

n × ( ˜∇ × ˜Escat) + 2πin0E˜_Tscat= 0.

The boundary condition becomes at the top n = [0, 0, 1]T_:

   ∂ ˜E3 ∂ ˜x − ∂ ˜E1 ∂ ˜z ∂ ˜E3 ∂ ˜y − ∂ ˜E2 ∂ ˜z 0   = −2πin0   − ˜E1 ˜ E2 0  − ˜r

For the planar wave with θ = 0, φ = 0 and ψ = π/2, at ˜z = ˜M, we have:

˜ E1= ˜E3= 0, ∂ ˜E2 ∂ ˜z − 2πin0 ˜ E2= −2πi2n0B0exp(−2πin0z),˜

9 11 3 8 18 10 0 20 24 21 23 25 22 12 13 17 19 ₁₅ 14 16 5 1 4 2 7 ₆ 26

Figure 6: (Re)numbering nodes for a 3x3x3 mesh with 3 periodic boundary conditions.

which means σ = −2πin0 and ˜r = −2πi2n0B0exp(−2πin0z)˜ in this case.

And at ˜z = 0, we simply have ˜ E1= ˜E3= 0 −∂ ˜E2 ∂ ˜z − 2πin3 ˜ E2= 0, so σ = −2πin3 and ˜r = 0.

4.3 Orientation

Important is the direction of the tangential vectors dening edges. The direc-tions should be the same for each (tetrahedral) element sharing the same face. The local vertex renumbering in these corresponding elements dene the direc-tions uniquely. Also the numbering of vertices on the matching faces at the corresponding periodic boundaries should count for the correct orientation of the edge and face DOFs. This is achieved by rst renumbering the corners (at least six periodic boundaries), edges (at least 4 periodic boundaries) and faces (at least 2 periodic boundaries) as depicted in Figure 6.

4.4 Constraints

Grating retain the full solution vector xT _{= [y z]}T_{, where y = Qx are the}

unconstraint DOFs and z = P x are the constraint DOFs. The system including the constraint equation is solved

 CH_{AC 0} −P C I   y z  =  CH_{b − C}H_Ad P d  .

1 2 4 3 0 0 1 2

Figure 7: Linear subspace

Let's eliminate the quasi-periodic boundary condition x0 = αx1 where x1 is a

unconstrained DOF. If the base function Φ0is dened as the linear combination

of two original base functions Φ0 = αφ0+ φ1, and copying the other base

functions Φj = φj+1, j ≥ 1 we can write E = y0Φ0+ y1Φ1+ · · · = y0αφ0+

y0φ1+ y1φ2· · ·. Then the DOFs y and z are directly related to the DOFs x as

y0= x1

y1= x2

· · ·

z0= x0= y0α,

The (non-square) matrix Q has a the identity at the unconstrained rows and zeros elsewhere. For more general boundary conditions, the constraints are written as

xi=X

j6=i

αijxj+ βi,

where αij, βi are complex numbers, there is also a direct correspondence (i.e.

Cij = αij and di= βi).

An illustrative example with one Dirichlet and one periodic boundary condition A standard nodal nite element space W contains 5 base functions, numbered from 0 to 4. The new space V ⊂ W is a linear subspace, in which the base functions are linear combinations of the original ones. Note that base function

1 in W is located on a Dirichlet boundary and is therefore eliminated. The

transfer matrix C in V = W C is C =       0 0 1 0 0 0 1 0 0 0 1 0 0 0 1      

and P and Q are dened as P =  1 0 0 0 0 0 1 0 0 0  , Q =   0 0 1 0 0 0 0 0 1 0 0 0 0 0 1  .

In case of vector base functions, the periodic boundary conditions also involve

the diculty of orientation. The pair of DOFs, say x0 and x1, need to be

identied and corrected for orientation. Since the mesh is unstructured, we have to consider the orientation of both element and edge. By smart numbering all the orientations are equal at both sides: rst number the eight corner vertices, the vertices on the 10 edges and the other vertices on the 6 sides of the domain, respectively.

4.5 Finite elements

Below a number of vector elements from the rst family are discussed. 4.5.1 Whitney

A DOF is dened on an edge, which are shared by the neighboring elements. On a element, the middle nodes are shared this way. So, the vertices have no DOFs and middle points (i = 0, ..., 5) are the location of the edge element DOFs. 4.5.2 Second order Nedelec

20 shape functions dene the solution. So, in addition to the Whitney functions, 6 extra edge functions (i = 0, ..., 11) and 8 face functions (i = 12, ..., 19) are implemented. Note that i = 0 and i = 1 are a Whitney and a linear edge functions dened on edge e0, and so on. Face nodes are located in the centre of the faces. The delta property for this element is forced by the explicitly do the linear combination as described above.

4.5.3 Third order Nedelec

45 shape functions are implemented. Adding 6 edge functions (i = 0, ..., 17), 16 face functions (i = 18, ..., 41) and 3 volume functions (i = 42, ..., 44) to the second order Nedelec shape functions. The volume DOFs are part of the element itself, so no node (or location) is associated to these DOFs.

4.6 Constraints

The constrained DOFs are eliminated from the mass matrix and right-hand side, but kept in the system as linear combination of other DOFs. What is the impact of applying transfer operators on these systems.

Suppose we have two constraint equations x = Cy+d, ˆx = ˆC ˆy + ˆdand a transfer

operator G that maps to these equations x = Gˆx. Thus,

Cy + d = G( ˆC ˆy + ˆd). The constrained preconditioned system can be derived

A(Cy + d) = b AG( ˆC ˆy + ˆd) = b AG ˆC ˆy = b − AG ˆd (G ˆC)TAG ˆC ˆy = (G ˆC)T(b − AG ˆd) ˆ CT(GTAG) ˆC ˆy = ˆCTGTb − ˆCT(GTAG) ˆd ˆ CTA ˆˆC ˆy = ˆCTˆb − ˆCTA ˆˆd

Which means that constraining is eectively a multiplicative combination of two

preconditioners ˆC,G and the constraint part of the system remains decoupled.

Can we construct a transfer operator for the constrained variables y and ˆy?

Note that d = G ˆdholds, so

Cy = G ˆC ˆy.

If multiplied with the transpose of C, the right-hand side operator is the con-strained transfer operator

i

0

i

1 0

j

1

j

right left

Figure 8: Reposition of elements by constraining. Left and right multiplication. j₀

1

j

i₀ i₁

Figure 9: Periodic base functions.

CHCy = CHG ˆC ˆy.

In the case of Dirichlet boundary conditions CH_{C = I, the identity matrix. For}

quasi-periodic boundaries this product is a diagonal matrix and in the general

case a full matrix. Note that in all cases CH_{C is invertible. So the operator ˆG}

transferring y to ˆy is

ˆ

G = CHC
−1CHG ˆC.

In the assembly stage the problem of a non-unit diagonal element is circum-vented by overwriting the entries in the operator, this is automatically achieved by the in-lined constraining as done in Grating.

Left and right multiplying the transfer matrix with the constraints will reposi-tion the contribureposi-tions. The quasi-periodic boundary condireposi-tions as dened above will ensure that the factor is always 1 (i.e. the complex constraint is only a shift in the angle).

4.7 Preconditioner

As an example, the NMGAV algorithm uses the spaces ∇W = {∇λ0, ∇λ1, ∇λ2, ∇λ3},

where λiare the barycentric coordinates and W = {φ0, φ1, φ2, φ3, φ4, φ5}the

ap-proximation space of lowest order Whitney elements (e.g. φ0= λ0∇λ1−λ1∇λ0,

7 3 9 8 6 4 0 1 5 2

Figure 10: Division tetrahedron

G =         −1 1 0 0 0 −1 1 0 −1 0 1 0 −1 0 0 1 0 −1 0 1 0 0 −1 1         .

The elements on the coarse grid are divided into 8 sub-elements

with numbering 0-4-6-7, 4-1-5-8, 6-5-2-9, 7-8-9-3, 4-5-6-8, 7-4-6-8, 6-8-5-9 and 7-8-6-9, respectively.

The numbering of the sub-elements is such that the sorted numbering on the outer sides of the sub-elements is preserved. So, for sub-element 0-4-6-7, sides

s0, s1 and s3 has the local 0-1-2, 0-1-3 and 0-2-3 numbering respectively. This

is needed for identifying the new boundary faces.

5 Results

The results listed below are obtained from simple congurations with analytical solutions.

5.1 Two-layer problem

The two-layer problem with is a prescribed incident plane wave that reects back a part of this wave. The exact solution (Pomplun, page 10) is used to determine the convergence behavior.

The solution can be written as a number of plane waves E0 = EA₀+ EB₀ and

x z B0 A0 A₁ n₀ n l 2 l 1 1 0

Figure 11: Exact solution two-layer case

EA₀(x) = A0exp(2πin0dA₀· x),

EB₀(x) = B0exp(2πin0dB₀· x),

EA1(x) = A1exp(2πin1dA1· x),

where

A0· dA₀ = B0· dB₀ = A1· dA₁ = 0

The outgoing wave is reected back at the same angle

(dB0)x= (dA0)x, (dB0)y= (dA0)y, (dB0)z= −(dA0)z,

and from the refraction preservation we have

n1(dA₁)x= n0(dA₀)x, n1(dA₁)y= n0(dA₀)y,

(dA₁)z= −q1 − (dA₁)2

x− (dA1)

2 y.

We take the minus sign here, because we are interested in (the physical) waves going toward z → −∞.

At the interface z = l1we have n × (E0− E1) = 0, thus

And n × (∇ × E0− ∇ × E1) = 0. Using the identity ∇ × E = 2πinid × E, we have

n0(dA0× A0)x,yexp(2πin0(dA0)zl1) + n0(dB0× B0)x,yexp(2πin0(dB0)zl1)

= n1(dA1× A1)x,yexp(2πin1(dA1)zl1).

For the six unknowns A1 and B0 we have six linear equations, which can be

solved:        

(dA₁)x (dA₁)y (dA₁)z 0 0 0

0 0 0 (dB₀)x (dB₀)y (dB₀)z

α 0 0 −β 0 0

0 α 0 0 −β 0

0 −n1(dA₁)zα n1(dA₁)yα 0 n0(dB₀)zβ −n0(dB₀)yβ

n1(dA₁)zα 0 −n1(dA₁)xα −n0(dB₀)zβ 0 n0(dB₀)xβ                 (A1)x (A1)y (A1)z (B0)x (B0)y (B0)z         = [0, 0, γ(A0)x, γ(A0)y, γn0(dA0× A0)x, γn0(dA0× A0)y] T

with α = exp(2πin1(dA1)zl1), β = exp(2πin0(dB0)zl1)and γ = exp(2πin0(dA0)zl1).

This linear system can be solved numerically (LU-decomposition) for a given

incident wave. For example n0 = 1, n1 = 2, l1 = 1₂ with perpendicular

in-cident wave dA0 = [0, 0, −1] T _{and A0 = [0, 1, 0]}T_{, we get B0 = [0, −}1 3, 0] T _, A1= [0, −2₃, 0]T_{, dB} 0= [0, 0, 1] T _{and dA} 1 = [0, 0, −1] T_.

Note that the solution for this problem can be obtained without solving the linear system above. We can reduce the problem a one-dimensional problem for TE and TM type waves. If we rotate the coordinate system around the

z-axis, the y-components of the electric and magnetic eld decouples, Ey =

ey(z) exp(±ikx0x0) and Hy = hy(z) exp(±ikx0x0). In our problem the

rota-tion is x = x0_{cos(φ) and y = x}0_{sin(φ). Furthermore, we have e}inc

y (z) =

pyexp(−2πin0cos(θ)z)for the incoming wave.

Then use the Fresnel reection and transmission coecients (part 2; Chew, page 48) and impose constraints on the interfaces to obtain the amplitudes of

the ey(z)in the dierent layers.

For TE waves the traverse components Ez and Ex are then calculated from

102 103 104 105 10−2 10−1 100 number of DOFs relative error

Figure 12: ◦ error H(curl)-norm and × L2-norm two-layer problem. Uniform mesh, Whitney elements.

5.2 Cavity problem

Figure 13 shows the result the simulation of the cavity problem. It has the solution

Eexact= p exp(2πid · x),

where p · d = 0 and d is a unit vector. It satises the Maxwell equations with

F = 0and ni= 1.

The solution to the second problem is an element in R1=_{u(x) = a + b × x | a, b ∈ C}2 _:

Eexact=   y − z 4 − (x − z) x − y  ,

with ∇ × ∇ × Eexact_{= 0and should be resolved exactly, when using Whitney}

elements (Monk, page 139).

Eexact=   z2_{− xy} x2_{− zy} y2_{− xz}  

is a function in R2 = (P1)3_{⊕ S2 in which elements of S2 are polynomials like}

102 ₁₀3 ₁₀4 ₁₀5 ₁₀6 10−3 10−2 10−1 number of DOFs relative error H curl 10−2 ₁₀−1 10−3 10−2 10−1 mesh size relative error H curl

Figure 13: Error H(curl)-norm cavity problem. ◦ Dirichlet; × x-periodic; 

x, y-periodic; ∗ Dirichlet, θ = 20◦. Domain K = 0.8, L = 0.2 and M = 0.2.

Uniform grid, Whitney elements.

∇ × ∇ × Eexact _{= [−3, −3, −3]}T_{.where a perfectly conducting region is placed}

in air.

5.3 Cylinder problem

The exact scattered solution Es _{for the innite cylinder problem is (Balanis,}

page 603) E_zs(r, φ) = −E0 ∞ X n=0 nin Jn(ka) Hn(1)(ka) H_n(1)(kr) cos(nφ),

where 0 = 1, n = 2, n > 0. Hankel function of the rst kind Hn(1) and Jn

Bessel function of the rst kind. k is the wave number and a is the radius of the cylinder.

5.4 Performance

Using the MKL PARDISO solver, we obtain a boost of factor 80 compared to the non-optimized PETSc direct LU solver. Note that PARDISO performs reordering for less ll-in (METIS), super node pivoting approach and maximum weighted matching algorithm. The convergence performance of the NMGAV algorithm is depicted in Figure 15 simulating the two-layer problem.

6 Appendix

105 10−0.6 10−0.5 10−0.4 error H curl number of DOFs

Figure 14: Error in H(curl) of cylinder problem.

0 5 10 15 20 25 30 35 40 45 10−8 10−6 10−4 10−2 100 102 number of iterations l2

−norm residual error

without

with

Figure 15: Convergence solver for two-layer problem with and without NMGAV-preconditioner.

6.1 Mapping

F (ˆx, ˆy, ˆz) =   x(ˆx, ˆy, ˆz) y(ˆx, ˆy, ˆz) z(ˆx, ˆy, ˆz)   dF =    ∂x ∂ ˆx ∂x ∂ ˆy ∂x ∂ ˆz ∂y ∂ ˆx ∂y ∂ ˆy ∂y ∂ ˆz ∂z ∂ ˆx ∂z ∂ ˆy ∂z ∂ ˆz    (dF )−T =    ∂ ˆx ∂x ∂ ˆy ∂x ∂ ˆz ∂x ∂ ˆx ∂y ∂ ˆy ∂y ∂ ˆz ∂y ∂ ˆx ∂z ∂ ˆy ∂z ∂ ˆz ∂z   = 1 det(dF )      ∂z ∂ ˆz ∂y ∂ ˆy− ∂z ∂ ˆy ∂y ∂ ˆz − _∂y ∂ ˆx ∂z ∂ ˆz− ∂y ∂ ˆz ∂z ∂ ˆx  _∂y ∂ ˆx ∂z ∂ ˆy − ∂y ∂ ˆy ∂z ∂ ˆx −∂x ∂ ˆy ∂z ∂ ˆz− ∂x ∂ ˆz ∂z ∂ ˆy  ∂x ∂ ˆx ∂z ∂ ˆz− ∂x ∂ ˆz ∂z ∂ ˆx −  ∂x ∂ ˆx ∂z ∂ ˆy − ∂x ∂ ˆy ∂z ∂ ˆx  ∂x ∂ ˆy ∂y ∂ ˆz− ∂y ∂ ˆy ∂x ∂ ˆz −  ∂x ∂ ˆx ∂y ∂ ˆz− ∂x ∂ ˆz ∂y ∂ ˆx  ∂x ∂ ˆx ∂y ∂ ˆy − ∂y ∂ ˆx ∂x ∂ ˆy      det(dF ) = ∂x ∂ ˆx  ∂y ∂ ˆy ∂z ∂ ˆz − ∂z ∂ ˆy ∂y ∂ ˆz  +∂y ∂ ˆx  ∂z ∂ ˆy ∂x ∂ ˆz − ∂x ∂ ˆy ∂z ∂ ˆz  +∂z ∂ ˆx  ∂x ∂ ˆy ∂y ∂ ˆz − ∂y ∂ ˆy ∂x ∂ ˆz  (φ ◦ F )(ˆx) = φ(F (ˆx)) n ◦ F = (dF ) −T_nˆ |(dF )−T_n|ˆ , τ ◦ F = dF ˆτ |dF ˆτ | Z V φ(x) dV = Z ˆ V φ(F (ˆx))| det(dF )| d ˆV Z e φds = b Z a φ(s(t)) ks0(t)k dt S = r(T ), Z S φdS = Z T φ[r(u, v)] ∂r ∂u × ∂r ∂v du dv

6.2 Cross product and curl operator

a × b =   a2b3− a3b2 −(a1b3− a3b1) a1b2− a2b1   a · (b × c) = (a × b) · c (a + b) × c = (a × c) + (b × c) a × b = −b × a, (a × b) × a = a × (b × a) a × (b × c) + b × (c × a) + c × (a × b) = 0 ∇ × ∇g = 0 ∇ · (∇ × E) = 0 ∇ · (A × B) = B · (∇ × A) − A · (∇ × B) A × (B × C) = (A · C)B − (A · B)C ∇ × (g F ) = (∇g) × F + g(∇ × F ) ∇ · (g F ) = g∇ · F + (∇g) · F A · (∇ × ∇ × B) − (∇ × ∇ × A) · B = ∇ · (B × ∇ × A − A × ∇ × B) ∇ × E =   ∂ ∂yE3− ∂ ∂zE2 −( ∂ ∂xE3− ∂ ∂zE1) ∂ ∂xE2−∂y∂ E1  

a

b

a x b

Figure 16: Cross product

References

[1] Peter Monk, Finite Element Methods for Maxwell's Equations, Oxford University Press Inc., New York, ISBN 0 19 850888 3

[2] Yu Zhu and Andreas Cangellaris, Multigrid Finite Element Methods for Electromagnetic Field Modelling, The IEEE Press, Wiley-Interscience, ISBN 13 978 0 471 74110 7, ISBN 10 0 471 74110 8

[3] Weng Cho Chew, Waves and Fields In Inhomogeneous Media, Van Nos-trand Reinhold, New York, ISBN 0 442 23816 9

[4] Constantine A. Balanis, Advanced Engineering electromagnetics, Wiley 1989, ISBN 0 471 62194 3

[5] Krzysztof A. Michalski, Electromagnetic Scattering and radiation by sur-faces of Arbitrary Shape in Layered Media, Part 1 & 2, IEEE Trans. An-tennas and Propagation, Vol. 38 Issue 3, 1990

[6] Pavel Solin, Karel Segeth, Ivo Dolezel, Higher-Order Finite Element Meth-ods, Chapman & Hall/CRC, London, ISBN 1 58488 438-X

[7] Jan Pomplun, Rigorous FEM-Simulation of Maxwell's Equations for EUV-Lithography" (Diplomarbeit), Konrad Zuse Institut, 2006

[8] Ronald Rook, Notes on solving Maxwell equations Part 2: Green's function for stratied media, CASA report 10-??

PREVIOUS PUBLICATIONS IN THIS SERIES:

Number Author(s)

Title

Month

11-06

11-07

11-08

11-09

11-10

M.E. Hochstenbach

L. Reichel

W. Hoitinga

E.H. van Brummelen

M.V. Ugryumova

W.H.A. Schilders

M.V. Ugryumova

W.H.A. Schilders

R. Rook

Fractional Tikhonov

regularization for linear

discrete ill-posed problems

A discontinuous Galerkin

finite-element method for

a 1D prototype of the

Boltzmann equation

Evaluation and comparison

of FEM and BEM

for extraction of

homogeneous substrates

Efficient simulation of

power MOS transistors

Notes on solving Maxwell

equations Part 1:

Finite elements method

using vector elements

Jan. ‘11

Jan. ‘11

Febr. ‘11

Febr. ‘11