c.5559

,ie:
:

tot-5,: .3?
a." sand.
I! .‘i.!\

 

f!

2. 1%,.

RH?
amt»

”hm—3h

.43....

.n

_
. .

E. :

fl “in?0

x

. 1‘! .
t... .éflmmmwmmwu

1mg:
.1? V
3.? .1

$4.13.
. .5 n...

, w 3%“?!
my.» , ,w

.,w

. ..
$1 an...
.A ..
"3.... .
. .Ezflmx

953.:
sky Stir.

Ru

.1... . ‘1 ..

v. v1»».‘tc
:9. IV.»
‘1!- $9.?

.3 .

’ ‘ 9

.. ...u S .
u “my Mekfllmmialgdrifl
2 ‘ 2.91.8.» Z. . .

£22234.

kw?

 

q" P816

2 .LIBRARY
260:5 Midligan State
mversity

This is to certify that the
dissertation entitled

GENERALIZED FINITE ELEMENT
METHOD FOR ELECTROMAGNETIC
ANALYSIS

presented by

Chuan Lu

has been accepted towards fulfillment
of the requirements for the

Ph. D. degree in Electrical and Computer

 

 

Engineerinl

 

/ Major Professor’s Signature

02-04-2008

 

Date

MSU is an affirmative-action, equal-opportunity employer

 

 

_.—n-c-I-I-I-I-I-o-o-u-o-o-o-I-I-o-o-I-I_I-I-I-l-o—o...- -

 

PLACE IN RETURN Box to remove this checkout from your record.
TO AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

5/08 K:IProj/Acc&Pres/ClRCIDateDue.indd

GENERALIZED FINITE ELEMENT METHOD
FOR ELECTROMAGNETIC ANALYSIS

By

Chuan Lu

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY
Electrical and Computer Engineering

2008

ABSTRACT

GENERALIZED FIN ITE ELEMENT NIETHOD
FOR ELECTROMAGNETIC ANALYSIS
By
Chuan Lu

The generalized finite element method (GFEM), first introduced by Babuska, is
a partition of unity-based solver for scalar partial differential equations (PDEs). To
date, they have been applied extensively to the solution of elliptic and parabolic
PDEs. This technique is a generalization of a host of well known methods for
solving PDEs, specially the classical finite element method(FEM), element free
galerkin(EFG), hp clouds, etc. The main goal of this dissertation is for developing
a similar methodology for vector electromagnetic problems. Developing a solution
to these problems necessitates addressing the following problems: (i) The vector
nature of the problem and the different continuity requirements on each component
imply that basis functions developed should share similar characteristics; (ii) The
basis functions have to be able to represent divergence free electromagnetic fields
(in a source free region). (iii) Development of appropriate boundary conditions to
truncate the computational domain is necessary. (iv) High condition number of the
resulting system also plagues GFEM solver, as it does other high order solvers. So-
lution to these problems, and the developments of the GFEM solver is presented
here for both time and frequency domains. In any case, the h- and p- convergence

Of the method is presented.

To My Dear Parents

iii

ACKNOWLEDGMENTS

I’d like to take this opportunity to express my gratitude to all my friends and col-
leagues who helped me to complete this dissertation and who have played important
roles in my life at Michigan State University.

At first, I should thank my parents for your continuous and unconditional support
for my life and study since childhood. Especially in the last few months you take
care of my daily life and help me concentrate on my work.

I like to thank my advisor, Dr. Shanker Balasubramaniam, for leading me to
this interesting research area. You taught me every skills useful for the project; you
helped and encouraged me with open-minded discussion and creative ideas; you also
helped me improve my presentation and writing skills. In four years in this group,
I have grown up from a student with little electromagnetic background and poor
English to an experienced researcher with skills necessary for a successful career. It
would have been impossible without your help and continuously pushing me to be
better.

Thanks to Dr. Leo Kempel for your valuable suggestion in group discussion and
your contribution to HPC system, which is very helpful for my research. Thanks to
Dr. Nyquist for your reviewing my research in group meeting and comprehensive
exam. Thanks to Dr. Guowei Wei for giving me suggestion for my career. Thanks
to Dr. Edward Rothwell for consolidating my background and theory in advanced
electromagnetics research, your book was crucial for portion of my research. Special
thanks to Dr. Satish and Dr. Lalita for taking care Of my research and study in my

first year in US.

iv

A thank you also goes to my teammates, Jun Gao, Jun Yuan, He Huang, Gregory
Kobidze, Tuncer Ozgur, Naveen Nair and Vikram Reddy, Jorge M. Villa, Carlos A.
Jaramillo and Pedro Barba. I am always happy to discuss with you and learn a
great deal from it. You are one reason I like being in the office. It’s also exciting to

have dinner, play table tennis and practice with you together.

TABLE OF CONTENTS

LIST OF FIGURES ................................ viii

LIST OF TABLES ................................. xi

KEY TO SYMBOLS AND ABBREVIATIONS ................. xii
CHAPTER 1

Introduction ..................................... 1
CHAPTER 2

Scalar Generalized Finite Element Method .................... 8

2.1 Introduction ................................ 8

2.2 Basic framework of GFEM solver .................... 10

2.2.1 Formula of basis function ..................... 10

2.2.2 Evaluation of inner products .................. 17

2.3 Imposing boundary conditions ...................... 20

2.3.1 Neumann Boundary Condition: ................. 21

2.3.2 Dirichlet Boundary Condition: .................. 21

2.3.3 Global Radiation Boundary Condition: ............. 24

2.3.4 Local Radiation Boundary Condition: .............. 28

2.4 Numerical experiments .......................... 29
CHAPTER 3

Vector Generalized Finite Element Method .................... 49

3.1 Introduction ................................ 49

3.2 GFEM solver for homogeneous domain ................. 51

3.2.1 Statement of the Problem .................... 51

3.2.2 Basis functions .......................... 52

3.2.3 Error bound ............................ 57

3.2.4 Imposition of boundary conditions and the resulting bilinear

form ................................ 64

3.2.5 Numerical experiments ...................... 66

3.3 GFEM solver for piecewise homogeneous domain with flat interface . 78

3.3.1 Basis function ........................... 78

3.3.2 Proof of the completeness of the vector basis function ..... 83

3.3.3 Numerical experiments ...................... 84

3.4 GFEM solver for piecewise homogeneous domain with curved interface 87

3.4.1 Basis function ........................... 87

3.4.2 Numerical Experiments ...................... 95

3.5 Overcoming linear dependence and condition number issues ...... 105

vi

CHAPTER 4

Time domain Generalized Finite Element Method ................ 110
4.1 Introduction ................................ 110
4.2 Formulation for Time domain GFEM solver .............. 111

4.2.1 Discretization of time domain wave function .......... 111

4.2.2 Causal system for time domain wave equation ......... 115

4.2.3 Error bound of Time domain GFEM solver ........... 116

4.3 Analysis Of stability ............................ 119

4.4 Iterative Solver .............................. 122

4.5 Numerical Experiments .......................... 125
CHAPTER 5

Conclusion and areas for future work ....................... 140

vii

Figure 1.1
Figure 2.1
Figure 2.2
Figure 2.3
Figure 2.4

Figure 2.5

Figure 2.6

Figure 2.7
Figure 2.8

Figure 2.9

Figure 2.10

Figure 2.11

Figure 2.12

Figure 2.13
Figure 2.14
Figure 2.15
Figure 2.16

Figure 2.17

Figure 2.18

Figure 2.19

Figure 2.20

LIST OF FIGURES

Meshes for FEM ...........................
Construction of the computational domain. ............
One dimensional partition of unity function. ...........
Construction of one dimensional Basis function. .........

Subdivision of the domain of integration based on piecewise
weight function Wi- .........................

Subdivision of Dglj based on the number of patches at different
location. ...............................

Integration in the sub-patch E?)- that overlaps with boundary of
computation domain I‘. .......................

Geometry transform .........................

(a) Definition of geometry for imposing the boundary integral; (b)
Definition of the geometry for application of the PML. ......

Error in the L2 norm of the numerical and analytical solutions of
the PDE with the Neumann boundary condition. .........

h, p convergence of the numerical scheme applied to the solution
of a PDE with Neumann boundary conditions. ..........

Error in the L2 norm of the numerical and analytical solutions of
the PDE with the Dirichlet boundary condition. .........

h, p convergence of the numerical scheme applied to the solution
of a PDE with Dirichlet boundary condition. ...........

h, p convergence of GFEM-BI (CFIE) ...............
Error in E—field over a range of [ca with GFEM-BI (EFIE) . . . .
Error in E—field over a range of [ca with GFEM-BI (CFIE) . . . .

Comparison between numerical and analytical data obtained for
scattering from a perfectly conducting cylinder of radius 2.7 A:
(a) the electric field at Fe; (b) Echo-width of the cylinder.

Comparison between numerical and analytical data Obtained for
scattering from the coated perfectly conducting cylinder if radius
2.88 A: (a) the electric field at Fe; (b) Echo-width Of the cylinder.

Comparison between numerical and analytical data obtained for
scattering from a coated perfectly conducting cylinder Of radius

3.84 A: (a) the electric field at I‘e; (b) Echo-width. ........

Comparison Of data obtained for scattering from an L—shaped Ob-
ject using GFEM-BI and GFEM-PML. ..............

Geometric description Of the wedge ................

viii

20

25

31

39

40

Figure 2.21
Figure 2.22
Figure 2.23
Figure 2.24
Figure 2.25
Figure 3.1
Figure 3.2
Figure 3.3

Figure 3.4

Figure 3.5
Figure 3.6

Figure 3.7
Figure 3.8

Figure 3.9

Figure 3.10
Figure 3.11
Figure 3.12
Figure 3.13

Figure 3.14
Figure 3.15
Figure 3.16
Figure 3.17

Figure 3.18
Figure 3.19
Figure 3.20
Figure 3.21
Figure 3.22
Figure 3.23

Analytical result and graph of error for scalar wedge problem . . 44
Error plot of scalar wedge problem ................. 45
Error of FEM and GF EM solver .................. 46
Error plot of FEM solver ....................... 47
Error plot Of GFEM solver ...................... 48
Patch cross the material interface .................. 52
Vector plots of homogeneous local approximation function v,- . . . 56

Cpnvergence graphs for vector Neumarm and Dirichlet BVPs by
uap .................................. 68

anvergence graphs for vector N eurnann and Dirichlet BVPs by
“up .................................. 69

Convergence graphs for 3D homogeneous vector Neumann BVPs . 71

Convergence graphs for 3D homogeneous vector Neumann BVPs

(Plane Wave) ............................. 71
Geometric description of the wedge ................ 73
Analytical result for scattering from a wedge and error plot for y

componnet .............................. 75
Analytical result for scattering from a wedge and error plot for y

component .............................. 76
h, p convergence for wedge problem ................. 77
Interface Of two media ........................ 79
Vector plots of inhomogeneous local approximation function V, . 82

Convergence graphs for 2D inhomogeneous vector Neumann and
Dirichlet BVPs ............................ 86

Convergence graphs for 3D inhomogeneous vector Dirichlet BVPs 87

Figure Of patch on the curved structure ............... 89
Local coordinate of patch on the curved structure ......... 91
Vector plots of inhomogeneous local approximation function V,

for curved interface .......................... 92
Geometric transform ......................... 94
Figure of patch on the curved structure ............... 96
h—, p— convergence of the local approximation .......... 96
Figure of rotated curved structure .................. 98
Error with different curved surface ................. 98
Geometric description of dielectric-coated PEC cylinder ...... 99

ix

Figure 3.24
Figure 3.25
Figure 3.26
Figure 3.27
Figure 3.28

Figure 3.29

Figure 3.30

Figure 3.31
Figure 3.32

Figure 4.1
Figure 4.2
Figure 4.3
Figure 4.4
Figure 4.5
Figure 4.6
Figure 4.7
Figure 4.8
Figure 4.9
Figure 4.10

Figure 4.11

Figure 4.12

Figure 4.13

Figure 4.14
Figure 4.15

Convergence graphs of GFEM solver with Dirichlet BC ...... 99
Analytical and numerical result of p component of the field . . . . 100
Analytical and numerical result Of a component of the field . . . . 101
Figure of patch on the curved structure ............... 102
Analytical and numerical result of large 2D inhomogeneous prob—

lem in the interface .......................... 103
Analytical and numerical result Of large 2D inhomogeneous prob-

lem out Of the interface ........................ 104
Condition number of system with differen order of basis func-

tions(p = 1,...,4) ........................... 108
Error of the solver with differen order of basis functions ...... 108
Condition number of system with different order of basis func-

tions(p = 1,...,10) ........................... 109
APSW function in time and frequency domain .......... 114
Stability 0f syStem With Nsamp = 4 ................ 123
Stability of system with Nsamp = 5 ................ 123
Stability of system with Nsamp = 4 ................ 124
Stability of system with Nsamp = 5 ................ 124
Analytical and numerical result with Dirichlet BC ......... 127
Analytical and numerical result with Neumann BC ........ 128
Analytical and numerical result with Impedance BC ....... 129
h— p- convergence of the spatial basis function ........... 130
Convergence Of different temporal scheme in 2D homogeneous

problem ................................ 132
Convergence Of different temporal scheme in 2D inhomogeneous

problem ................................ 133
Convergence of different temporal scheme in 3D homogeneous

problem ................................ 135
Convergence of different temporal scheme in 3D inhomogeneous

problem ................................ 136
Error with iterative solver ...................... 139
Number of iterative steps in each time step ............. 139

Table 2.1
Table 2.2
Table 2.3
Table 3.1

Table 3.2

Table 3.3
Table 4.1

LIST OF TABLES

Error of eigenvalues in rectangular waveguide ........... 33
Error of eigenvalues in circular waveguide ............. 34
Error of eigenvalues in coaxial waveguide ............. 34
Error in eigenvalues computed using tic],p ............. 72
Error in eigenvalues computed using 1131, ............. 72
Error in eigenvalues computed using 113,1, ............. 88
The coefficients of characteristic equation .............. 121

FDTD:

MOM:

FEM:

FMM:

PDE:

GFEM:

SPH:

DEM:

EFG:

MLS:

APSW:

PML:

SVD:

CFIE:

EFIE:

PEC:

EW:

TE:

TM:

MOT:

KEY TO SYMBOLS AND ABBREVIATIONS

Finite-Difference Time-Domain
Method Of Moments

Finite Element Method

Fast Multiple Method

Partial Differential Equation
Generalized Finite Element Method
Smooth Particle Hydrodynamics
Diffuse Element Method

Element Free Galerkin

Moving Least Square

Approximate Prolate Spheroidal Wave
Perfectly Matched Layer

Singular Value Decomposition
Combined Field Integral Equation
Electric Field Integral Equation
Perfect electrically conductor

Echo Width

Transverse Electric

'Ifansverse Magnetic

Matching On Time

CHAPTER 1

INTRODUCTION

Computational electromagnetics has grown over the last few decades, in terms of
both development of the underlying tools and application to the design of practical
systems. Typical methods of analysis are either based on integral equation or dif-
ferential equation. As is well known, method of moment(MOM) and other integral
equation—based method is expensive in terms of both memory and CPU-time. Aug-
menting this with fast multiple method(FMM) has considerably alleviated the cost.
In differential equation based methods, the most popular are the finite difference
time domain(FDTD) and finite element method(FEM).

The state of art of FEM tools [1] for electromagnetic analysis has grown by leaps
and bounds over the past few decades. Classical methods require an underlying
tesselation on which basis functions are defined. These basis functions are based
on a span of polynomials, have finite support, and obey conditions at inter-element
boundaries. For instance, Figure 1.1 shows the discretization of elliptical computa-
tion domain. Here, the edge of the meshes coincide with each other and they also

coincide with the boundary of the computation domain as well.

 

Figure 1.1. Meshes for FEM

In classical FEM, the basis function is interpolatory. FEM solvers using this

basis have found widespread applications in solving problem in mechanical injury,
stress strain, vibration analysis, heat transfer, etc [2, 3]. In the electromagnetics,
early application Of FEM did focus on directly using these space of functions [4, 5].
However, it was discovered that these basis functions do result in spurious modes as
some conditions in the E and H are not satisfied. Furthermore, it is not possible to
readily impose the boundary conditions.

Fortunately, vector basis function, e.g. edge element vector developed in 19803
[6, 7] is free of all the previous shortcomings. The zero order basis function is
divergence free, as its tangential component is continuous across the edge Of the
mesh. In homogeneous domain, the normal component Of the basis function is not
continuous across the edge of the mesh, but the continuity of the normal component
of numerical solution is enforced by the variation formula of the vector wave equation.

High order vector basis function are defined in two ways: hierarchical [8] and
interpolatory [9]. Both have continuous tangential component cross the edge of
the mesh, but they are not divergence free. Until now, high order basis functions
have been continually refined [10, 11]. Likewise, basis functions that can be used
in a h, p—convergence setting have been presented and applied to several engineering
problems. Indeed, in a series of excellent papers, it has been shown that it is possible
to have true h, p—convergence [12, 13].

But there are some limitation of classical FEM, that may prove to be a bottleneck.
First, it requires a simplical structure. The edge of the mesh has to coincide with each
other and the edge of the mesh has to coincide with the boundary of the computation
domain as well. Having good quality meshes is also necessary for accuracy of the
FEM solver. Both properties make mesh generation laborious and time consuming.
When millions of unknowns are needed, time required for mesh generation might
be over 90% of total CPU time. Developing h,p—convergence meshes or analyzing
time varying phenomena that requires re—meshing at every time instant can be even

more laborious. Another handicap of classical FEM is that the ansatz space used

to approximate the local behavior is a span of polynomials. It implies that solving
problems in non-Lipschitzian domain is difficult. This is due to the fact that the
field or its derivative are singular at the geometric singularity. This can not be
simulated well by a span of polynomials even in dense meshes. If analytic local
behavior is known, then it might be possible to use functions other than polynomials
to approximate local behavior.

Those limitations result from the inter-element continuity requirements. The
polynomial basis functions can enforce the requirements by its conforming property.
There are still a variety Of methods to realize it with non-polynomial basis. One
is proposed by [14, 15]. Around the geometric singularity, there are one layer of
meshes supporting the singular basis functions. Those basis functions are conformal
to the neighboring regular elements. In this approach, the size of such meshes has
to be almost one wavelength. Another category of totally different methods, hybrid
methods, can construct the local function space in each mesh separately based on
the physical analysis of the solution and enforce the inter element continuity require-
ments by additional constraints [16, 17]. In [16], lagrange multiplier is introduced
to enforce it. This results in a larger linear system cause of the addition unknowns
related to the multipliers. Stability consideration also links the number of unknowns
related to basis functions and the number of unknowns related to multipliers.

In this dissertation, I will introduce a new method that is free of the men-
tioned shortcoming. This method will be referred to as the generalized finite ele-
ment method(GFEM)[18]. Comparing this method with regular FEM and hybrid
methods, regular FEM enforces inter element continuity requirements by conformal
finite elements but lacks the capability of physics of the problem. Hybrid methods
have the freedom of choosing local approximation but additional constraints between
elements are necessary. In one of hybrid methods, the function space for multiplier
has to be defined carefully. GFEM shares both of the above advantages, e.g. any

function space can be chosen to fit for the prior knowledge with conformal elements

enforcing the inter elements continuity requirements intrinsically which is free of
above drawbacks.

GFEM belongs to the category of meshless methods [19]. The starting point of
Meshless method is smooth particle hydrodynamics (SPH) method by Lucy in 1977
[20], The objective of meshless method is to eliminate the mesh generation proce-
dure by constructing the approximation around the nodes. After that, the rate of
publication about this topic was very modest until about the 19903. In 1992, Nay-
roles [21] used the moving least square (MLS) approximations in a Galerkin method
and called his method diffuse element method (DEM). Two years later, Belytschko
[22] redefined and modified this method, and called it element free galerking(EFG).
In 1996, there was a leap forward in the development of meshless methods with the
introduction of GFEM by Babuska and Melenk [18, 23]. This method is based on
partition of unity ftmction and propose to a very powerful extension of the EFG
approximation. In the same year, Duarte and Oden [24] extend the concept of par-
tition of unity function to the MLS shape function and called it hp-cloud method.
With this improvement, the EFG method based on MLS can be regarded as special
instance of GFEM.

While meshless methods are well developed in some research area such as me-
chanical engineer, civil engineer, etc, they are largely used to solve scalar equation
and usually solve the elliptic or parabolic PDEs. More specifically, it is used for
specialized application such as crack propagation, etc [25]. The application of this
technique to solving problems in electromagnetics has not been extensive. Prin-
cipally (this is not a complete list), research has been conducted by [26—37]. The
work in this area by Shanker and his colleagues revolved around developing meshless
methods for solving the diffusion equation in both the frequency and time domain as
applied to non-destructive evaluation [27—29]. The basis functions used in this analy—
sis relied on the element free galerkin method (EFG) [19]. The method, particularly

its scalar implementation, has been gaining a foothold in the research landscape

insofar as application to magnetic field analysis. Recently, they have explored the
viability of suitably modifying the EFG method to enable the analysis of vector
fields [31] and have deve10ped meshless-PMLS to enable the analysis of open region
problems [38] within the context of the EFG method. However, while the EFG
method can be thought of as a subset of GFEM, it does not lend itself readily to
hp—adaptivity, whereas this is inherent in other GFEM methods.

The above solvers are developed for the scalar differential equation. As was men-
tioned earlier, extending this to the vector differential equation that are encountered
in EM results in spurious modes. These modes can be readily identified by studying
the eigenfunction of same Operator for defined geometry. Two kinds of modes can
be distinguished. One is divergence free and another has high divergence value. The
latter non physical modes plagued EM simulation and is known as spurious modes.
The main reason of spurious problem is that the scalar FEM can not impose the
divergence free criteria and satisfy the continuity requirements intrinsically. In clas-
sical FEM, several ways have been suggested for years to remove the spurious modes
in FEM such as penalizing the divergence [39, 40] or using divergence free criteria to
reduce 1 / 3 of dependent unknowns in three dimension[41—43]. The introduction of
Whitney basis function, that is divergence free and tangentially continuous, provided
the method to overcome this deficiency.

It follows from the above discussion that the vector GFEM solver should satisfy
the following criteria:

The basis function developed should be able to approximate the field which
is (i)divergence free (ii)satisfying all the continuity requirements on the material
interface.

Except for the spurious problem, another problem is to impose essential boundary
condition in GFEM solver when it’s trivial in classical FEM. It also results in the
difficulties to implement the boundary integral techniques with EFIE formula.

Another problem that plagues GFEM, as it does any higher order method is the

condition number Of the resulting system. Indeed, this is more of a problem here as
opposed to classical FEM as the space of basis functions is not interpolatory.
Tsukerman [35] tried to develop the vector basis function which can be used in
a GFEM solver, but he did not demonstrate the convergence of the GFEM solver
and the basis function can not handle inhomogeneous domain.
The objective of this dissertation is to develop vector GFEM method for EM
analysis. The new method should be able to simulate the vector fields with high

order accuracy. The principal goals of this research are as follows:

0 Develop a scheme for implementing GFEM for the Helmholtz equation.

0 When solving closed domain problem, an adaptation of Nitsche’s method is
developed for implementing the Dirichlet boundary condition. When solving
open domain problem, the hybrid GFEM-BI method for domain truncation is

developed.

0 Develop vector basis function to approximate the divergence free field in ho-

mogeneous domain and prove bounds on their convergence.

0 Develop basis functions for inhomogeneous domain that will enable the sat-
isfaction of boundary conditions across flat and curved material interfaces.
These conditions will be on the tangential components of the fields and their
normal derivatives, and the normal component Of the fields and their nor-
mal derivatives. We will mathematically prove the accuracy of the proposed
method and demonstrate the convergence of the method for two dimension

and three dimension.

0 Develop a preconditioning scheme tO solve the ill—conditioned system resulting

from the high order system, and derive error bound and cost Of this scheme.

0 Develop time domain GFEM solver and prove and demonstrate the conver-

gence of this method.

o Propose high order time stepping schemes based on approximate prolate
spheroidal wave (APSW) and analyze the stability of the high order time

stepping scheme. Demonstrate the convergence of the proposed scheme.

This dissertation is organized as followings: in the chapter 2, basic framework of
GFEM method is introduced. The techniques for imposing different boundary con-
ditions for both Open domain and closed domain are elucidated. In chapter 3, vector
GFEM method is developed for homogeneous domain and piecewise homogeneous
domain with flat and curved material interface. The convergence of the solver will
be proved and demonstrated in both two dimension and three dimension problems.
A preconditioning scheme is also introduced in the same chapter with error bound
and running cost provided. In chapter 4, high order time stepping schemes are pro-
posed and demonstrated in two dimension and three dimension. Stability of high
order time stepping schemes are analyzed. Finally, in chapter 5, a surmnary of the

work and area for future research are collected.

CHAPTER 2

SCALAR GENERALIZED FINITE ELEMENT METHOD

2.1 Introduction

Intuitively, GFEM solver works as follows: the domain being considered is parti-
tioned into a union of patches or a “partition of unity,” and on these patches, the
local approximation is constructed using a span of ftmctions [44]. Thus, the repre-
sentation of the ftmction is achieved via two functions; one that is defined on the
partitions of unity and the other on each of the patches. The basis functions describ-
ing the unknowns inherits the higher order nature of approximation from the local
basis functions and the smoothness Of the functions defined on the partition of unity.
As with classical FEM, using a span of different local approximations in different
regions is also possible. Thus, GFEM retains several features Of classical FEM, pro-
vides additional flexibility in terms of functions that are used and obfuscates the
need for a simplical partition of the domain.

The mathematical foundation of GFEM has been laid out in great detail [18, 23,
45], and it has been shown that h—, p- and hp-adaptivity is easily achieved. Likewise,
the efficacy of using a space of harmonic functions as local approximants have been
demonstrated [45]. Nevertheless, when some meshless solver has been developed for
solving the diffusion equation in both the frequency and time domain as applied to
non-destructive evaluation [27-29], it does not lend itself readily to hp—adaptivity.

Proper understanding of the sources of error and the means through which one
may control them is important. We have found that in most of the implementation,
reason that the convergence is not of the same order of the underlying scheme is
largely due to the improper imposition of boundary conditions. Proper imposition
of boundary conditions within a meshless scheme is challenging. Unlike classical

FEM, the space Of approximation functions are not interpolatory. This poses severe

challenges in imposing Dirichlet boundary conditions. More specifically, one cannot
use the Lagrange multiplier technique as the approximation spaces have to obey
the inf-sup condition, and it is not always possible to construct such spaces for
meshless methods. This deficiency is directly linked to the difficulty in truncating
the computational domain using boundary integrals. This is due to two facts; (i)
to solve the hybrid problem, one typically defines an auxiliary set of basis functions
and unknowns to represent the tangential components of the fields [1]. It implies
from the above discussion [46] that these basis functions, together with those used in
the interior, should satisfy the inf-sup (Babuska-Brezzi) condition, and (ii) there are
practical situations wherein one uses a first kind Fredholrn integral equation as the
boundary is not closed. In these cases, the BI enforces a Dirichlet type condition.

Thus, the principal content of this chapter is four fold:
1. A scheme for implementing GFEM for the Helmholtz equation is presented.

2. An adaptation of N itsche’s method for implementing the Dirichlet boundary

condition is developed.

3. Hybrid GFEM-BI technique for domain truncation for both open and closed

domains is presented.

4. Methodology wherein local boundary conditions can be integrated with
GFEM-more specifically, the perfectly matched layer (PML), is presented.
The development of this technique is a by-product of the need to have an ad-
ditional modality of validating the results obtained by the GFEM-BI scheme.

The rationale for embarking upon this specific problem is as follows: (a) This
approach presented in this chapter (basis functions/ means to impose boundary con-
ditions, etc) can readily used for solution of quasi-static electromagnetic phenomena
and scalar wave equations; (b) It permits us to work out several mathematical and

numerical hurdles—the principal being the application Of Dirichlet boundary condi-

9

tion and the accurate evaluation of integrals; (c) it is equally important to under-
stand the manner in which boundary integral based techniques can be hybridized
with this scheme. The advantage of hybridizing GFEM with BI is readily appar-
ent; it imposes an exact boundary condition for open domain problems, and the
computational cost can be amortized using recent advances in the integral equation
technolog, namely the fast multipole method or a host of FFT-based schemes.
This chapter proceeds along the following lines; in the next section, we formulate
the problem in detail. Here, we introduce the concepts of GFEM, discretization of
the domain, basis functions, methods for integration, and methods for implementing
various boundary condition. The last includes different types of boundary integral
techniques and a local absorbing boundary condition. Next, we demonstrate the
accuracy and convergence of the GFEM and GFEM-BI via a series of analytical
comparisons. We shall also demonstrate the accuracy of this scheme by comparing
the results obtained against those obtained by truncating the domain using a PML
as an absorbing boundary condition. Finally, we will demonstrate the advantages of

GFEM solver by scattering analysis with geometric singularity.

2.2 Basic framework of GFEM solver

2.2.1 Formula of basis function

Consider a multiply connected domain (2 whose interior boundaries are denoted by
60:: I‘ = U,- I‘i. It is assumed that this domain is embedded in a domain 96 and
its exterior boundary Fe is defined as Fe := (Te 0 5—). Interior to the domain {2, the

function u(r) satisfies

(v- [a(r)V] + w27(r)) um = fl!)
3,- {u(r)} = 9,;(r) for r E F,- (2-1)

Be {u(r)} = 960‘) for r E Iwe

10

In the above equations, it is assumed that r 6 Rd, 88 and B,- are differential op—
erators, and gz-(r) is the function that is imposed on Pi. Here, d = 2,3, and a(r)
and 7(r) are material parameters. The function of interest, u(r), is used to denote
the 2 component of either the electric or magnetic field. The rationale of defining
Fe explicitly is to impose appr0priate boundary conditions that enable the analysis
of scattering problems. The parameters a(r) and ,6 (r) can stand for either the per-
mittivity or permeability, depending on the variable that u(r) represents. Solution
to this problem using the standard finite element method requires an underlying
tesselation on which basis functions are defined. These basis functions are based
on a span of polynomials, have finite support, and Obey‘conditions at inter-element
boundaries. For instance, Whitney elements that are typically used in computational
electromagnetics satisfy either tangential or normal continuity across inter-element
boundaries. Higher order basis functions based on these elements have also been
presented and have been continually refined [10, 11]. Likewise, basis functions that
can be used in a h, p-convergence setting have been presented and applied to several
engineering problems. Indeed, in a series of excellent papers, it has been shown that
it is possible to have true hp—adaptivity [12, 13]. On the other hand, meshless meth-
ods attach a patch or volumina to each point whose union forms an open covering of
the domain. The local shape functions are constructed within each domain. Several
different flavors of these methods exists [23]. In this dissertation we will base our
development on the generalized finite element method (GFEM) [23].

Our presentation of the fundamentals of basis functions is a repetition of those
in [45, 47]. GFEM is based on a set of N nodes located at r,- in the vicinity of the
, domain 9 such that {1'13 6 Rd : L; 6 9,2' = 1, - - - ,N}. Associated with each of these
nodes is a patch or volumina denoted by $2,; of size h,- such that Q C Cg := U,- Q,-
and Q,- = {r 6 Rd: ||r — rill S hi} C Rd. Specifically, a patch 52,- is defined as
Q,- = ®g=1 522(k), ng) = {1"(k) 6 RI r106) — TU“) IS hgk)}. Figure 2.1 describes

such a construction. Typically, there are no restrictions on the shape of the domain.

11

TO a large extent, these are chosen depending on the underlying basis functions.

 

Patches

Figure 2.1. Construction of the computational domain.

Associated with each patch are basis functions that will be used for Galerkin
testing and source. The basis function is a product Of two functions, if),- (r) and
v,- (r): functions w,- (r) form a partition of unity subordinate to the cover (22;, and a
space of functions v,- (r) E span {222710)} that are local to the domain Qi-

The global approximation of the variable is then a space of functions denoted
by V,;(r) = ¢i(r)v,- (r). The basic theory of the GFEM method using this space of
functions was originally developed by Babuska and Melenk [45] and is summarized

below. Note that the definitions listed below are key to some of the proposed tasks.

Definition 2.1 :
Let Q 6 Rd be an open set, and let Q; be an open cover of Q satisfying a point-

wise overlap condition

3M 6 N Vr e 0 card {ilr e 522'} _<_ M (2.2)

12

Definition 2.2 :

Let {2b,} be a Lipchitz partition of unity subordinate to the cover {92-} satisfying
the following conditions: (i) sappwi) C (2,- for all i; (ii) 2,: 1b,- 5 1 on 9; (iii)
iWillLoo<Rd) S Coo; and (iv) [iVl/liHLoomd) S CV/diam(wi). Here Coo and
CV are two constants. Then {2b,} is called an (M, COO,CV) partition of unity
subordinate to {Qil- The partition of unity is said to be of degree k if it; E Ck(le).

The covering sets 9,- are called patches.

Assume that the space of functions V(r) = span {viz-vi} is given. Let the function
nap := 2, 1/12’112' 6 V C H162) be an approximation to u(r). The upper bounds on
the error in the function and its gradient have been derived in [45].

The choice of the functions viz-(r) and vi(r) can vary from patch to patch, and
is usually dictated by the geometry of the problem. More specifically, the functions
vi(r) can theoretically be modified to include logarithmic singularities. However, the
first step Of implementing this method is the construction of the patches Q,- Of size
h,- that form the Open cover for the domain. In general, we allow hi , for j = 1,- - - ,d
to be different in each dimension. In this dissertation, we have restricted ourselves
to rectangular domains. As the patches 9,- are to be used in a Galerkin scheme, it
is necessary for Uz- Q,- 3 Q. In other words, for any r E 52 there exists at least one

patch (2,; that contains r. The following steps is utilized to realize it [47]:

o N nodes are chosen randomly in the computation domain and h,- is chosen
initially as the size of the ith patch. With circle-shaped patch, the h,- can be
the radius Of the circle; with rectangle-shaped patch, h,- can be half Of the side
length.

0 Much more denser testing points are uniformly distributed in the computation
domain. The testing points has to be at least 10d times denser than the nodes.

Each testing point should be covered by at least one patch. If any testing point

13

is not covered by one patch, then look for the node which is closest to it. Then

increasing the size of the patch hj to make jth patch covers that testing point.

c When each testing point can be covered by at least one patch, true patch size
hi can be derived as h;- = ah,- where a > 1. Number a is used to assure
the enough overlapping between the patches so that the smoothness of the

partition of unity function 11),; can be maintained.

This will ensure that for r E (2 there exists at least one patch such that r E 92':
Next, the partition of unity functions 11),; are defined in each 52,. By definition,
2,; rbi(r) = l V r 6 52,-. To construct these functions, we use a localized version of
Shepard’s method [48] that relies on defining a function W,- (r) with respect to each
characteristic point r,- 6 9,- such that Wi(r) = 0 \7’r 6 802°. In other words, these
functions are different from zero only for r E 52,-. Several choices for these func-
tions exist; the most common are B-splines of different orders, Gaussian, regularized
Lagrangians, etc. The choice of the functions W,- depends on shape of the patch.
In this dissertation, we restrict ourselves to rectangular domains, and therefore, to
B-splines of different orders. Then in d-dimensions we construct a ftmction W2: (r)
using a tensor product of l-D functions; i.e., W;(r) = H?:1W((r — r[)/(h[)).
Next, the functions 1b,- (r) are defined by v),- (r) =—- (“SUD/(291,60,- Wk (1')), wherein
C,- := {9i 6 Cfllfli f] it; is 0}. Here defined partition of unity function 1,0,- (r) sat-
isfy the criteria 2,- r/Jz-(r) = 1Vr. Figure 2.2 shows a one dimension example of the
partition of unity function. As is evident from the figure, patches overlap with each
other, and there is one partition of unity function in each patch. The summation of
all partition Of unity functions in one patch equals one at any point. We can also
create the partition of unity function based on MLS functions [49], when high or-
der basis functions are used, the error of Shepard and MLS functions are almost the
same though MLS functions are much more complicated in terms of formulation [49].

The partition of unity functions are consistent to first order and higher order

14

7.01 42 1P3

 

0.0

w

91 :22 93

Figure 2.2. One dimensional partition of unity function.

approximations are Obtained by defining higher order local ftmctions, vi(r). As
mentioned in [45], local approximating functions can be chosen such that they bet-
ter capture the local behavior of the current. However, in this chapter, we have re-
stricted our choices to Legendre polynomials of order P. Thus, in any domain Q,- the
functions vz-(r) = span {H521 fpa) ((40 - rl’))/(h§”)) ] where pit) = 0. - -- .10.

Figure 2.3 shows a one dimension example of the construction of the scalar basis
function. The three local approximation functions are chosen as 0th, 13t and 2nd
order Legendre polynomials. Basis ftmctions, tlvi(r)v]"(r), decay to zero at the edge
of the patch, which assure the continuity of the numerical solution. Since there is
no restriction on the choice of the basis functions, they can be chosen such that
they best approximate the local physics. For instance, around the smooth geometry
plane wave can be chosen as local approximation, around the geometric singularity
eigenfunction of the singular structure can be used, etc.

It is well known that bilinear form of the differential equation (4.1) on H1(Q) is

.A(u,w) = (f, w) ‘v’ (u and w) E H1(Q) (2.3a)

A(u,w) = ~(aVu,Vw) + w2('yu,w) (2.3b)

where (~, -) denotes the standard imier product in L262). It is apparent that we

—--——_——-— ---_--W-——_-—-

 

 

 

 

 

 

 

 

1112 ”’3 $2125
\ x / I i _
/ i .-'/\x
l // \V/
* i= /
L/
$22 ‘02 1,1120%
-_-1_--_____7 _____________________
l '.
l I", if ]
i .2 .-' A f'\ A
~\ ./A
\_/ \C/

Figure 2.3. Construction of one dimensional Basis function.

have not yet incorporated the boundary conditions; the means by which various
boundary conditions are incorporated into the variational form will be dealt with
exhaustively in Section 2.3. Using the definitions of basis functions elaborated upon
thusfar, u(r) can be now approximated as u(r) = 2971,19, where N is the total

degrees of freedom. Using Galerkin method, the discrete version of this bilinear form

can be written as

A (um mm) = A (mafia-v3”) (2.4a)

(swim) = (ivy-22;") (24b)

Above equations (2.3) indicate that the Val-m, need to be evaluated. Here we use

16

 

“-1

r
I
I
I
-—--L——-—

.3.

-—-—---——_——-----—I

 

_——- —-——

——-_1

 

 

_-__J..-_-

 

1
I
I

.g

3.
I
I
l
I
T
I
I
I
—--*---d

b------—-..--------u

in--- --—--——— ----l

 

 

i-n--‘———-
I————-'————

 

Figure 2.4. Subdivision of the domain Of integration based on piecewise weight
function Wi.

the following formula to evaluate it:

V117" =V(1/Iz'v?)

 

 

 

W.
=V( 2 on)
Zakec‘, Wk ’
_VWi aneo, Wk’vf‘ - Wi aneo, Vka? (25)
(2919607; Wk)2
Wi
+ an
29,360,]: Wk 2

Note, that the equations (2.4) indicate that the interaction between any two domains
form block matrices. However, what is crucial in the process is the evaluation of the

inner products; this will be dealt with next.
2.2.2 Evaluation of inner products

The basis functions l/z-(r) that were introduced earlier are piecewise rational func-
tions. This is immediately apparent when one examines the construction of ibi(r);
a different number of functions W;(r) may contribute to 1b;(r) at different location.

This implies that the quadrature rules must be such that they respect these dis-

17

 

---r--+--+---
I l I
l I l

 

 

 

 

 

———+——+——4-—-
I I I
I I II I I
~--‘--* -4r-‘I l
l l l
I I l r r I
l I I.
T I v v | '
—-—r-~+—++-+ I
I I II I I
I .—+—L4_J.—q-J___
MI I II I I
I
l

 

 

| J [I 1

 

I I I
l l I

 

 

 

Figure 2.5. Subdivision of D];- based on the number of patches at different location.

continuities. To overcome these deficiencies, we have followed a method similar to
the one proposed by Schweitzer [50]. Here the domain of integration is subdivided
into smaller domains where the integrand can be represented as a product two ra-
tional functions. Note that the integrands comprise of products of both t/jz-(r)v]n(r)
and their derivatives. To illustrate the decomposition Of the domain, consider two
patches {22- and It; that overlap; the domain of overlap is denoted by Cij 2: f2,- 0523-.
Further assume that in each domain, the weight functions Wi(r) and Wj (r) are
polynomials of degree 1. Then in (2;, the function W;(r) is piecewise rational in
(l + l)d—domains. As the function viz-(r) is a combination of both Wi(r) and Wj (r)
for r E 9i, we can partition the intersection Cij into Ci]: = Un Dznj disjoint subcells
as in Figure 2.4, and w,- (r) is piecewise rational in each subcell. The process becomes
more involved when more than two patches overlap. Such a scenario is illustrated
in Figure 2.5. Assume Eijk := 93- fl Qj 0 9k denotes the domain of intersection of
these three domains. It follows that Um F316 Q UT, 0]}- denotes the union of do-
mains wherein the function tbi(r) is piecewise continuous. Recursive identification

of domains where the function tb.,-(r) is piecewise rational ensures the appropriate

18

 

ii

 

 

 

Figure 2.6. Integration in the sub—patch E];- that overlaps with boundary Of com-
putation domain I‘.

construction of quadrature rules.

The above procedure works well for intersection between two domains Cij : =
Q,- M23: C 9. However, if F dissects C7; -, then it can no longer be considered a union
of rectangular domains; see Figure 2.6. As before, assume that Cij 2 Un DE}, and
in each D];- the functions w, and 1b]: are rational functions. The cells D];- are be
partitioned into two sets; those that are fully contained within {2 and those that
intersect I‘. For the former, we shall use the standard quadrature rules; the latter
can be be evaluated through adaptive quadrature. In this chapter, we have restricted
ourselves to boundaries that are circular. The geometric choice permits the division
into rectangular domains and curvilinear triangles as shown in Figure 2.6. As before,
we use standard quadrature rules in the rectangles, and coordinate transformation
and higher order integration in the triangles [1]. Thus, the prescribed procedure
enables higher order evaluation of all the inner products in the variational form.

Here is the brief explanation of the higher order integration in the curvilinear

triangle. As shown in Figure 2.7, a curvilinear triangle in my-domain can be mapped

19

T

 

 

Figure 2.7. Geometry transform

to a regular triangle in fry-plane. Transform functions

it =$(€,n) (2.6)

y =y(£.n)

are known. So to the related two points (xp,yp) and (Epfllp) in cry-plane and {77-

plane. Given the definition of J acobian matrix of the transform

8.1: 3
m: ,g%
35 d
and [J] is Jacobian, then
“(551): yp) =U(€p, rip)
utxpvp) =[Jl‘1 ~ «an,» (28)

drdy =|J|d§dn

2.3 Imposing boundary conditions

It is evident from the above exposition, that (i) the basis functions V; are NOT in-

terpolatory; (ii) (2 C Un supp {Vn}. These two facts make imposition of boundary

20

conditions more difficult than those for classical h, p convergence FEM methods.
In this section, we will discuss the imposition of both the Nuemann and Dirich-
let boundary conditions, as well as the truncation of the domain using boundary
integrals.

2.3.1 Neumann Boundary Condition:

First, consider the differential equation in (2.1) with B,- {u(r)} = a-un = an-Vu(r) =
gn(r) Vr E I‘. Here, n denotes the outward pointing normal from the domain 9.

This boundary condition implies that the bilinear form is to be modified as

An, w) = (f, w) — [P dI‘gn(r)'w(r) (2.9)

From the above, it is apparent that it is sufficient for the trial and basis functions
to be in H162). There are no additional constraints, and the basis do not have
to satisfy the boundary conditions. Hence, incorporation of Neumann boundary

condition is no different than that in classical FEM.

2.3.2 Dirichlet Boundary Condition:

Next, assume that B; {u(r)} = u(r) = g(r) Vr e I‘, i.e., Dirichlet boundary condition
explicitly imposes the values of u(r) on the boundary of the domain. Alternatively,

the problem may be cast as follows: “find u(r) 6 H162) such that u(r) = g(r) V r E

I‘” and the bilinear form for w E H1(Q)
A(u,w) = — [P dra<r>w(r)%z- Vu(r) + (raw) (2.10)

This statement is not very different from that posed for standard FEM albeit with
a couple Of differences: in classical FEM (1) w 6 H662) which implies that the
integral over the boundary vanishes; (ii) the basis functions are interpolatory. Hence,
imposing the boundary conditions is tantamount to modifying the linear system. In

GFEM, both the trial and test function do not satisfy either of these properties.

21

Thus, imposing Dirichlet boundary conditions is not as straightforward. This has
been a topic of considerable discussion [51].

Several Methods attempted are introduced here. (i) to hybridize meshless meth-
ods with classical FEM; (ii) use a penalty function method; (iii) use a Lagrange
multiplier technique, and (iv) use Nitsche’s method. Of the four methods, we have
chosen to use Nitsche’s method to impose the boundary condition. While we shall
discuss this method detail, we shall also dwell briefly on other three as it will shed
light on the constraints on imposing a global boundary condition.

When hybridizing meshless methods with classical FEM, FEM meshes is gener-
ated around the boundary. When far away from the boundary, EFG basis functions
are constructed. On the interface of two basis functions, FEM basis are as usual and
in meshless methods, shape functions take care of the consistency of the approxima-
tion.

The formulation using penalty function is stated as: u(r) 6 111(9), find the

solution to

Ah(u,w) = 7-"), the e VN (2.11s)

Ah(u,w) = — (aVu,Vw) +fl/I‘dl‘ aw (2.11b)
j: = , (11‘ 2.11

h (fw)+fi/P gw < c)

This method has two advantages. The dimension of the final system is not increased,
the system is symmetric and positive definite. But the Dirichlet is imposed weakly
and results in ill-conditioned matrix system.

Solving the problem using Lagrange multipliers involves finding a solution

(a, A) E H (O) x H *1/2(I‘). The formulation can than be stated as follows: given

22

u(r) 6 H162) and A E H—1/2(I‘) find the solution to

Aha, my, p) = 5, v (w and i.) 6 111(9) x H-1/2(r) (2.12a)
Ah = — (aVa, Vw) + w2 (711,10) + (A,w) + (11,11) (2.12b)
fh = (f, w) + (gnu) (2-126)

Solution to the discrete version of this equation in the GFEM setting is not trivial. It
is well-known that for this problem to converge, it is necessary for both the interior
and multiplier spaces to fulfill the Babuska—Brezzi condition. The Babuska-Brezzi

condition can be expressed as:

(V-u,A)

I
71f uEH1(Q)HUHH-1(Q)H’\ii

AEH“1/2(I‘)

 

Sup S y > 0 (2.13)

H’1/2(I‘)

while it’s easy to find A when using classical FEM, but it is difficult to design such
multiplier space that satisfies this condition [52—54], especially for the space of ba—
sis functions being considered here. Moreover, Lagrange multiplier based techniques
lead to indefinite systems. An alternative to this could be a stabilized-Lagrange mul-
tiplier technique [55]. However, in what follows, we use the Nitsche’s technique for
imposing the Dirichlet boundary condition. This method is related to the stabilized
Largrange multiplier technique with two advantages: (i) it is relatively straightfor-
ward tO implement in a numerical scheme, and (ii) one does not need to define an

additional space of functions in H ’1/2(I‘). The method proceeds as follows: find

23

an approximate solution such that u(r) 6 VN C H162) such that

Ah(u,w) =.7-'h Vw E VN (2.14a)
Ah(u,tv) = — (aVu, V'w) + (1)2 (ya, to) +/de1‘ aun - Vw (2.14b)
+/ drr awn-Va —[3/ da: aw (2.14c)

F F
fh = (f, w) +./I“ dragn-Vw — B/Pdl‘ gw (2.14d)

where B is chosen such that it guarantees coercivity. Rigorous estimates exist for 6

[50]. As is also apparent from the above equations, the resulting system is symmetric.
2.3.3 Global Radiation Boundary Condition:

The above exposition has a significant impact on the development of a global bound-
ary condition. In standard FEM—BI expositions [1], one defines equivalent currents
on Fe, and uses the radiation boundary integral to impose either the electric field
or the magnetic field or a combination of both. The equivalent current are of both
the electric and magnetic types. As u(r) represents one field (either the electric or
magnetic field), either the magnetic or electric currents can be easily obtained. One
typically prescribes basis ftmctions for all r 6 F3 to represent currents. However,
from our preceding discussion it is apparent that basis functions prescribed on R;
have to belong to H ‘1/2(I‘e), and these functions together with those used in the
interior have to satisfy the Babuska-Brezzi condition. While such a space can be
easily developed in the case of standard tesselation, it is perhaps not possible for
GFEM. Therefore, prescribing additional unknowns on the boundary is ruled out.
The method that we use for hybridizing will depend on whether Fe = 80 or Fe C 80.
In the former case, the domain of integration encloses a volmne, and in the latter, it
is open. When the domain is closed, it is sufficient to prescribe conditions such that
the solution is unique; i.e., one does not excite the interior resonance modes. While

this is a solved problem, it is implementing this condition within the context of

24

GFEM that causes problems. If Fe C 69, then imposing conditions is considerably
more challenging, as imposing the BI is similar to prescribing a. Dirichlet boundary
condition. This is challenging on two counts; (1) the basis ftmctions V,- are NOT
interpolatory; (ii) 0 C Un supp {Vn}. Therefore, techniques that were used earlier

to overcome this problem [1] need to be modified to impose the boundary condition.

 

 

 

 

 

 

 

 

 

 

0.5x,2
l
0.5x2
—O.5x;
——(.5x,2
—(.5x,2—-0.5x,‘ 0-5xl 0.5x,2
(a) Boundary Integral do- (b) PML domain description

main description

Figure 2.8. (3.) Definition of geometry for imposing the boundary integral; (b)
Definition of the geometry for application of the PML.

Our formulation proceeds as follows: Assume that just inside the boundary, one
can define another surface I‘s that completely encloses all the inhomogeneities in

0; see Figure 2.8(a). For the purposes of discussion, assume that u(r) refers to the

25

electric field. Then Vr 6 Fe,

v(r) = amen) + r: {u(r)} (2.153.)
anus) = anuincm + 1c {u(r)} (2.15b)
£ {u(r)} := — [I‘ dI‘s (6n1u(r')g(r,r') — u(r’)6n;g(r,r')) (2.150)
c {u(r)} := [F drs (an,u(r’)ang(r, r’) — u(r’)g,(r, r’)) (2.15.1)

where an, and an are used to denote the normal derivatives with respect to the
primed and unprimed coordinates, respectively, t = —a x 2, and ui"c(r) denotes
the incident field. Denoting r = r — r' where r E R2, the Green’s functions may be
written as g(r,r') = fngflclrl) and gt(r,r’) = k2(t % t’)g(r,r’) +f’ % (VVg(r,r’)) -t
where t and f’ are the tangent vectors at r and r’, respectively. These equations are
essentially derived from surface equivalence theorems (or Huygen’s principle). Note,
the surfaces Fe and F3 can be arbitrarily close to each other; however, if they are very
close to each other, the integral operators in (2.15) may be singular/hyper-singular,
and one should evaluate these with care. Techniques for doing so are similar to those
prescribed in [56]. Next, to incorporate the global boundary condition within the
differential equation solver, one needs to specify the differential operator Be {u(r)} in
(4.1). The simplest is to specify that Be {u(r)} = u(r) = nine“) +£ {u(r)} Vr 6 Fe.

This is, of course, a Dirichlet type boundary condition, and has to be incorporated

26

using Nitsche’s method; more specifically by implementing the following

Ah(u,w) =55}, Viv E VN (2.16a)
Ah(u,w) = — (aVu,Vw) + w2 (yu,w) + [Pdl‘ craft - Vw + [Pdl‘ awn - Va
(2.16b)
— )8 A df‘ aw — Adl‘afi {u(r)} a - Vw + fl/Pdl‘wfi {u(r)}

35, =( f, w) + / draui’wa - Vw — e / dP umcw (2-16c)
I‘ 1‘

Alternatively, one can specify Be {u(r)} = 6nu(r) = anuinca) + IC {u(r)} Vr 6 Fe.
As before, this is the Neumann type boundary condition, and can be incorporated by
appropriately modifying the variational formulation. However, it is well known that
both global boundary conditions do not yield unique solutions at all frequencies.
These frequencies correspond to the null spaces of the appropriate operators [1].
Among the several methods prescribed [1] to overcome this deficiency, one is to
combine the two Operators; i.e., use a combined field formulation. Thus, the method

proceeds as follows: find an approximate solution u(r) 6 H162) such that

Ah(u,w) = E, Vw e H162) (2.173.)
Ah(u, w) = — (aVa,Vw) + w2 (yaw) —jk/I‘ dI‘w(r)n(r)+ (2.17b)
./I‘ dF w(r) [jkfi {u(r)} + IC {u(r)}] (2.17c)

a = (Lw) — (r. dI‘w(r) [jkumc(r) + anumcm] (217d)

In the above exposition, we have essentially focused on imposing the boundary
condition at Fe. As in (2.1), another condition might need to be imposed in the
interior boundaries. However, as is apparent, the bilinear form for imposing these
conditions may be derived trivially using material presented thusfar. Evaluation of

the integrals over the boundary may be carried out using several different methods.

27

In our implementation, we proceed as follows: The domain Of integration (either F3
or I}; can be partitioned into a union of subdomains over with the basis function
is piecewise smooth. Then over each such subdomain, we use a Gauss-Legendre
quadrature. While this is not truly optimal and one may construct better quadrature

rules, we did Obtain convergent results.
2.3.4 Local Radiation Boundary Condition:

Our interest in implementing the local boundary condition is purely to develop
another measure of validating the results obtained using global boundary condition.
As we will show, the results Obtained using GFEM-BI converge exponentially to
analytical solutions. However, as these solutions are available only for canonical
problems, it is of interest to know that the fields obtained using the afore-developed
scheme and those obtained using a local boundary condition agree with each other.
There is a. wealth of local boundary conditions that are available [1, 57]. Here we
implement a perfectly matched layer (PML) within GFEM. The literature on PML
is extensive, and our goal is to present a rudimentary development that can be used
as a validation modality. The approach that we use to implement this relies on the
stretched coordinate system first introduced in [58]. Using this technique results in
a slight change in a(r) in that it becomes a tensor of rank 2, i.e., 5(r) = agaij (r).

The quantity ”7(r) and the non-zero elements of c’r(r) are

 

 

81(1‘) 02 (r) = 820)

0110') = 0‘0 82(1') 31(1‘) (2.18)

’i'(r) = 7'081(r)82(r)

where r = (2:1, 2:2), 00 and 70 are appropriate constants when all stretching param-

eters are one, and

31m = 1 wig—:1) s2<r> = 1 4335?? (2.19)

28

In these equations, 31 and 32 are stretching coordinates, and ”ii (r) denotes the
conductivity tensor of the domain. As is apparent, this tensor is chosen to be
diagonal. The domain of application of the PML is denoted using (2 p M L = {20 —Q;,

where the rectangular domains 90 and Q,- are defined as {20 = 18% x 92% and $2,- =

2:] x 25%; see Figure 2.8(b). The conductivity profiles are chosen to be zero inside the
F,- and vary quadratically in the direction of the outward normal to the boundary I‘;.
011 is nonzero in domain (—0.5x§, —0.5x]) x (—0.5z§,0.5x§) and (anion?) x
(-0.5$g,0.5x%), 022 is nonzero in domain (—0.5x%,0.5z¥) x (—0.5x3, —0.5:r%) and
(—0.5:r%, 0.53%) x (05:15:12,05163). Finally, the computational domain is truncated by
imposing the condition

a - Vus + jkcos(6a)us = 0 (2.20)

on F0. Here, 60 denotes the angle of perfect absorption and is chosen to be 9a = g
in our simulation. Note that in this formulation, the unknown is scattered field us.

With these changes, the bilinear form may now be written as

Ah(u3,w) = .7}, V11) 6 VN (2.21a)

Ah(n8,w) = — (6: - Vu3,Vw) + w2 (:yus,w) — jkcos(da)/Fdl‘ usw (2.21b)

fh = (f,'w) __ (aVZuinc,w) _ “J2 (Miriam) (2210)

2.4 Numerical experiments

In what follows, we shall present a series of numerical experiments that will serve
to demonstrate the accuracy and convergence of the method presented herein. In
all examples presented below, we have chosen d = 2 merely for demonstration pur-
poses, and extension to d = 3 will be exploited in the next chapter. First, we shall
demonstrate h, p convergence for problems wherein either the Dirichlet or Neumann
boundary conditions are specified. Next, we shall demonstrate similar convergence

for our hybrid GFEM-BI scheme and also demonstrate that our scheme is free from

29

corruption due to interior resonance modes. In these examples, global boundary
conditions are used to truncate the computational domain, and Dirichlet boundary
conditions are imposed in the interior of the domain. After that, we present a set
of results that analyze electrically large problems, and compare these against ei-
ther analytical results or against GFEM augmented with local boundary truncation
schemes. Finally, we will analyze the scattering from PEC wedge with mode-based
GFEM solver. While meshless methods offer a host of advantages, one significant
hurdle/ unsolved problem is the conditioning of the resultant linear system as the or-
der of the basis function increases. This is an issue when we are trying to demonstrate
h, p convergence. In these cases, we resort to an singular value decomposition(SVD)-
based solver. However, in the analysis of electrically large Objects, wherein we are
satisfied with an error in the Lg-norm = 1e-4, we use a non-stationary iterative
solver like TFQMR [59] with block preconditioners.

In the next two examples, we demonstrate h, p convergence of this method when
imposing either the Neumann or the Dirichlet boundary condition. The domain
of analysis of both problems are the same and defined as follows: the domain of
interest 9 = (0,1)2, and the boundary r = uf r,- where {r1 : r e 0 x (0,1)},
{F2:r€ (0,1) x 1}, {F3zr e 1 x (0,1)}, and {F436 (0,1) x 0}.

In the first example, the Neumann boundary conditions are imposed on all four
walls. More specifically, B,- {u(r)} = anuirlineI‘z- = gi(r). Denoting r = ($1,352),
gi(r) = 0, —2.97c0s(4a:1), 3.0272sin(3:r2),3cos(4:rl), for i = 1,--- .4. The above
boundary conditions permit analytical solution of the (2.1). In this experiment,
uniformly distributed nodes and rectangular patches are used. The size of each patch
is 1.5 times the distance between the nodes. The weight functions W,- (r) is a product
rooftops, and the approximation is a tensor product of Legendre polynomials. Two
sets of results are shown. First, we demonstrate the error between analytical and
numerical solutions in Figure 2.9 for h = 0.11A and p = 4. We also demonstrate h, p

convergence in Figure 2.10. As is evident from the graphs, the results are excellent.

30

 

Figure 2.9. Error in the L2 norm of the numerical and analytical solutions of the
PDE with the Neumann boundary condition.

 

10.1 r .
p/

-5 9:4

 

 

 

h—Distanoe between neighbour nodes

Figure 2.10. h, p convergence of the numerical scheme applied to the solution of a
PDE with Neumann boundary conditions.

 

Figure 2.11. Error in the L2 norm of the numerical and analytical solutions of the
PDE with the Dirichlet boundary condition.

10° .

 

10.1' p.=1// 1
10-2- 1

-3
NE 10 r p: 1
1o"- -

P=3

10'5r P=4 -

.6

 

 

 

-1
h-Distanoe between neighbour nodes

Figure 2.12. h, p convergence of the numerical scheme applied to the solution of a
PDE with Dirichlet boundary condition.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . p=1 p=2 p=3 p=4
fg’c‘lcelpal h=1/4 0.0244 5.49e-4 2.61e—6 9.67e-9
TE h=1/6 0.0151 1.07e-4 2.42e—7 3.84e-10
First8 h.=1/4 0.087 8.44e-3 1.62e4 3.02e-6
modes h.=1/6 0.0533 1.67e-3 1.58e—5 1.05e—7
Principal 5:1/4 0.0567 3.78e-4 3.10e—6 7.79e—9
TM Mode h=1/6 0.0285 8.14e-5 2.59e—7 3.20e-10
First8 h=1/4 0.1557 0.0142 3.20e—3 1.99e—5
Modes h=1/6 0.1122 5.27e-3 1.85e-4 2.28e-6

 

 

Table 2.1. Error of eigenvalues in rectangular waveguide

Next, the Dirichlet boundary condition is imposed on all four walls. More

specifically, B,- {u(r)} = “(flit/ref,- = gi(r). Denoting r = (151,122), g,- (r) =
sin(3a:2),0.1411cos(4:r1),—0.653651n(3:r2),0, for i = 1, - - . .4. All the parameters
used in the computation are the same as those used for imposing the Neumann
boundary condition. Again, the Figure 2.11 plots the relative value of the error in
the entire computational domain when using h = 0.11A and p = 4. Also, the errors
for different values for h and p are shown in Figure 2.12. It is evident that Nitsche’s
method for imposing the boundary conditions shows excellent convergence.

Next, we compute the eigenmodes in a rectangular waveguide. The dimensions of
the waveguide are chosen to be 1m x 1m. The convergence data (error in L2 norm)
for both TE3 and TMz modes are presented in the Table 2.1. As is evident from this
data, the results obtained converge rapidly with increasing order and refinement.

Next, we compute the eigenmodes in a circular waveguide. The radius of circle
is chosen to be 1m. Table 2.2 show the error of eigenvalues for both TE" and TMZ
modes. Again, the results Obtained converge rapidly with increasing h and p.

Then, we compute the eigenmodes in a coaxial waveguide. The radius of in-
ner and outer circles are chosen to be 0.5m and 1m. The error of eigenvalues are
presented in the Table 2.3. The results is excellent.

Next, we examine the accuracy Of imposing the boundary integral to truncate

33

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . p = 1 p = 2 p = 3 p = 4
magma h = 1/4 00191 4.43e—4 2.19e—6 6.09e-7

TE h = 1/6 0.0127 9.14e-5 5.72e—7 4.208—7
First 5 h = 1/4 0.1137 0.0157 1208-3 3269—5
modes h = 1/6 0.0817 4.52e—3 1.12e—4 2.11e-6
Principal h = 1/4 0.0677 9.478-4 9.97e-6 5.849-7

TM Mode h = 1/6 0.0351 2.1484 4.568-7 4.426—7
First 5 h = 1/4 0.175 0.291 0.040 2.123-4
Modes h = 1/6 0.131 0.090 4.21e—4 7.41e—6

 

 

 

Table 2.2. Error of eigenvalues in circular waveguide

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . p=1 p=2 p=3 p=4
11:42:29“ h=1/4 0.0311 2.91e—3 3.09e-4 3.69e—5
TE h=1/6 0.0205 5.61e—4 2.98e—5 1.26e-6
First5 h=1/4 0.1062 0.0307 5.35e—3 8.31e-4
modes h=1/6 0.693 0.0101 7.10e-4 5.51e—5
Principal h=1/4 0.2137 0.0103 2.41e-3 2.23e—5
TM Mode h=1/6 0.1369 2.96e—3 2.73e—4 6.36e-7
First5 h=1/4 0.2169 0.0402 8.64e—3 1.91e—3
Modes h=1/6 0.1529 0.0165 2.23e—3 1.31e—4

 

34

Table 2.3. Error of eigenvalues in coaxial waveguide

 

 

 

. i
t J
i l
. P=1 ,
r J

10-

L2 error

10' —

 

 

 

1 0"“

 

10'
h-Distance between neighbour nodes

Figure 2.13. h, p convergence Of GFEM—BI (CFIE)

the domain. To do so we analyze scattering from a perfect electrically conducting
(PEC) cylinder with radius 0.1A and the truncation boundary is placed 0.1A away
from the surface. Given the configuration of the problem, we need to apply the
BI (CFIE) (2.17a) and Dirichlet boundary condition on the truncation boundary
1“; and the inner boundary, respectively. We compare these results against that
Obtained analytically to obtain the h,p convergence graphs in Figure 2.13. As is
evident, the scheme presented in this chapter demonstrates the anticipated conver-
gence characteristics. The next, we analyze scattering from a PEC cylinder over a
range of frequencies. The outer boundary is truncated using the EFIE specified in
(2.16a), and the inner boundary is truncated using a Dirichlet boundary condition.
We know that truncating the outer boundary with the EFIE would lead to unique
values for all values of [to except those that correspond to interior resonance modes

of a cylinder with a PEC wall at the outer boundary. Thus, to satisfy ourselves that

35

this is indeed the case, and to show that the variational form that imposes the EFIE
formulation presented in this chapter is valid (note: EFIE is imposed using Nitsche’s
method), we analyze scattering from a cylinder over a range ka varying from 2.0 to
8.0 in steps of 0.05. Additionally, values of Ira that are within three digits of the
resonance frequencies corresponding to TM 7' modes of the cylinder whose radius
corresponds to that of F3 are chosen. In total, we ran this simulation for a total of
128 frequency points. In all cases, h = 0.1A and p = 3. Figure 2.14 compares the
error in the field values within the computational domain to analytical data. As is
expected, Figure 2.14 shows the results obtained are accurate except at resonance
frequencies. So, the next challenge is to ensure that our results are free of corruption
by spurious modes. Again, we analyze scattering from a PEC cylinder over a range
of frequencies with combined boundary integral formula. It is apparent from Figure
2.15 that the error is fairly constant over the entire band of frequencies. It should be
noted that in all the cases mentioned thus far, we are quoting the error in the field
values at a sufiicient dense set of samples inside the computational domain. These
error are NOT in the echo width data, as they tend to be significantly smaller.
Next, we demonstrate the applicability of this technique to various scattering
problems. In all three examples described, the incident plane wave propagates along
the it direction and is polarized along 2. First, scattering from a cylinder of radius
2.7). is analyzed. The source boundary F3 and fictitious boundary F8 are at 2.85A
and 3.0A, respectively, and h = 0.11A and p = 3. The analytical and numerical data
of electric field on fictitious boundary and echo width (EW) are compared in Figure
2.16. As is evident from these graphs, the results are excellent. Next, scattering
from dielectric-coated PEC cylinder is analyzed. The radius of the PEC cylinder is
2.82A, and the thickness of the dielectric coating is 0.06A. The relative permittivity
of dielectric is 2.0. The source boundary F3 and fictitious boundary Fe are at 2.94A
and 3.0A, respectively, where A is the free space wavelength and h = 0.12A and

p = 3. The analytical and numerical data Of the electric field on fictitious boundary

36

and echo width (EW) are compared in Figure 2.17. As is evident from these graphs,
the results are excellent. Finally, we analyze scattering from a coated cylinder that
is considerably larger. The radius of cylinder PEC and dielectric coating are 3.76A
and 3.84A, respectively. The source boundary F3 and fictitious boundary [‘5 are
at 3.92A and 4.0A, respectively. As before, A is the free—space wavelength, 5,- = 2,
h = 0.15A and p = 3. Our discretization in this case is considerably coarser, however,
as is evident from Figure 2.18 (a,b), both the fields on the fictitious surface and the
echo-width agree very well with the analytical data.

In the next example, we compare the results obtained using GFEM-BI with those
obtained by analyzing the same object using GFEM-PML. As is well known, the
principal advantage of the boundary integral is that realized by a reduced compu-
tational domain. This, of course, implies that the cost of application of the BI can
be reduced to something that scales almost linearly with the number of unknowns
on the boundary. This is indeed possible by augmenting the BI with acceleration
techniques, notably by the fast multipole technique [60]. The object that we choose
for simulation is a L-shaped dielectric scatterer. The length and width of each arm
is 1A and 0.3A, respectively, and the arms are oriented along the :i: and 3} directions.
The truncation boundary Fe for the BI is conformal to the scatterer and is at a
distance of 0.25A away from the scatterer. When employing the PML, the scatterer
is embedded in an rectangular domain of size 9,; = 5.4 x 5.4A2, and the thickness
of the PML is 0.95A. In both simulations, h = 0.08A and p = 3, and the relative
dielectric constant er 2 2. The incident field is 2 polarized and propagates along
It = —1/\/2 (:i: + 3}). Figure 2.19(a,b) show the fields Obtained by both methods. As
is evident, they are identical to each other. Indeed, the relative error in the field
values at a set of points in the domain is 2.3e—3, indicating that both techniques

compare very well with each other.

37

 

 

 

 

 

 

 

 

 

 

0.02 . T
0.018]
5 Q0”
I:
3‘ 0.014 ~
5
E 0.012 - i
0.01 - -
1%
C
3 0.008 -
E amn- .
“’ 0.004 - .
0.002 M -
O l 1 l 1 L
2 3 4 5 6 7 8

ka—)

Figure 2.14. Error in E—field over a range of he with GFEM-BI (EFIE)

-2 I T I T r r I

 

 

 

 

 

 

r[—-éenflbmwa
01mmmumm
-25~ ~
-hi 1
O
000 O 00 O 000
-3: " U 9 O 30 00 C Q 0 CO 0 -
000 C 0 00 O O o O O
O O O O o O
-4- _
-45- -
'52 '3 4 5 E 7 is 0 10 1‘1 1'2 13

 

up;

Figure 2.15. Error in E-field over a range of he with GFEM-BI (CFIE)

38

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 l I I
—Andyticelresult
----- Numericalresult
___EmortnEl
1i-
3
U
“N
m
[LIN 0 f
2
III
3
«30
...1P l‘ ‘ q, ng'fu' I"\r'€w:'mthr.flfi :—‘o
‘1“th ' "1 1' '. 'n ‘ ”it
31‘5“]fiff'), 6'” . i, 1‘ |l1||[t-so
1 1 I ' --eo
_2 ' r 1 m 1 I r L '70
0 50 100 150 200 250 300 350
NW)
(a) E2 on fictitious boundary Fe
45 I v v
—-Analy1lcalrosu|t
4°] ----- Numerieelresult
35-
30-
8
"U
€25b
1% 2°'
to
15-
10-
5»
0050100150200250300350

0 (degrees)
(b) Echo width (EW)

Figure 2.16. Comparison between numerical and analytical data obtained for scat-
tering from a perfectly conducting cylinder of radius 2.7 A: (a) the electric field at
Fe; (b) Echo-width of the cylinder.

39

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2 F T T I I T
—Amlytieolresult
----- Wreath
_ _ _ElrorinE‘
7:?
m
m
N S
I.“ s
8
in
.8
[] -2.
n.“\" [v ”I‘M; 1-30
I I h a r
“ ‘ ' "“10““
I
It '50
l l l m
250 300 350
(e) E2 on fictitious boundary F e
20 I I I I I l r
——Analylice|result
13‘ ----- Numerical result
16-
14*
a 12 -
g
13
Lu 8“
6..
4..
2..
0 L 4 1 l— l _.l_ L
0 50 100 150 200 250 300 350

0 (degrees)
(b) Echo width (EW)
Figure 2.17. Comparison between numerical and analytical data obtained for scat-

tering from the coated perfectly conducting cylinder if radius 2.88 A: (a) the electric
field at 1},; (b) Echo-width of the cylinder.

40

 

 

 

 

 

1- _
53‘
~11
in
IL!” 0 E
-0 g
m
3
-1. “I , r’|\ |' :r‘hl“: "‘ “(I "20
r ’1 I, I i \
‘11’1f]}[1l5,; [ ,' will)”
1], 1'1 . 1]. -40
l .‘ '
l
2‘ A A I I I l I _%
0 so 100 150 200 250 300 350

0 (degrees)

(a) E, on fictitious boundary I}

--— Analytical result
13 ————— Numerical result

 

Echo width (dB)
F3

 

 

 

 

"o 50 100 150 200 250 300 350
9 (degrees)

 

(b) Echo width (EW)
Figure 2.18. Comparison between numerical and analytical data obtained for scat-

tering from a coated perfectly conducting cylinder of radius 3.84 A: (a) the electric
field at Fe; (b) Echo-width.

41

 

-O.2 O 0.2 0.4 .6 0.8 1 1.2

o
x (A)

(a) GFEM-BI

 

-0.2 0 0.2 0.4 0.6 0.8 1 1.2
x (A)

(b) GFEM-PML

Figure 2.19. Comparison of data obtained for scattering from an L—shaped object
using GFEM-BI and GFEM-PML.

42

 

 

 

 

 

Figure 2.20. Geometric description of the wedge

Until now, the basis functions used in all the examples have been based on
tensor products of polynomials. The power of GFEM, however, is the ability to
include analytical behavior, if known, into the approximation space. To test this
idea, To test this idea, we analyze scattering from a wedge. The analytical response
of this structure is well known [61]. Indeed, [15] developed singular basis functions
that can be incorporated into the function space of the classical FEM methods (see
references therein for other approaches to constructing singular functions). In what
follows, we will explicitly incorporate the known functional behavior. Consider a
two dimensional wedge with angle it). that is located in a rectangular computational
domain (2 = (—1.0, 1.0) x (—1.0,1.0). The tip of the wedge coincides with the origin
of the computational domain and one side lies in the y = 0 plane. A field that is
polarized along the 2 and propagating along I: = —0.9659:t — 0.25883) is incident
upon the wedge, and the wavenumber k0 = 3; this geometric configuration is shown

in Figure 3.7.

43

 

(c) Graph of error based on eigen function

Figure 2.21. Analytical result and graph of error for scalar wedge problem

44

 

 

 

 

 

 

 

 

. —e— Special approximating
- - e - Regular approximating -
-1
1O :-
g I GM
c :2 \ \ \ \
“.1 , ‘Q \ G. :
.5 , r e ‘ s : a
mu 2 0 ~ :0?
‘8 10 L' \ \
g t m: ‘0
g .
T1 .
m 10 3:" m=4
-4
10 ‘ n ‘ ‘ ‘ ‘ ‘ ‘ ‘ ‘ * ‘ ‘ ‘4 w‘
101 1O2 103
Numberofunknowns

Figure 2.22. Error plot Of scalar wedge problem

It can be easily shown that the fields are given by

 

 

477E . . .
E2 :2” —0a ngJv(kp)szn[v(ain + 77 — a)]szn[v(¢ + 77 — 0)] ( )
v 2.22
2 mm, m = 1, 2, ' ' '
217 — (1

Modes Jv(kp)sin[v(¢ + 77 — (1)] are used as scalar local approximating function. The
infinite extent of the wedge is simulated by imposing analytical solutions for the
infinite wedge as Dirichlet boundary conditions. The relative error is computed using
L2 norm of error in the electric field at 100 x 100 sample points that are uniformly
distributed along :2: and :1) direction. In this example, tensor product of polynomials
are also used as local approximating function for the purpose of comparison. Figure
2.22 show the efficiency of new basis function. From those figures, we can find
that traditional basis function creat largest error at the tip of the wedge, it also

happens when we use traditional FEM method. But with new basis function, the

45

 

 

 

 

 

 

 

 

10,7”... - . --.-...,
' ——GFEM(30°) 5
Gb —e—GFEM(15°) ‘
. Tr:— o .
10.1 3‘- “:z—Z'Q—Q—fl —"_ GFEM(7.5 )
- * - FEM(30°)
E - e -FEM(15°)
- - * -FEM(7.5°)
—2
10 r
e .
m
g .
8 .
N4 10 \
10
102 103 104

Number of unknowns

Figure 2.23. Error of FEM and GFEM solver

error aroung the tip of the wedge is reduced efficiently. Figure 2.21 show the relative
error with respect to the number of unknowns. Obviously, error with new basis
function convergence much faster than those with tradition basis functions.
Finally, we compare the error of scattering from wedge by GFEM solver and
FEM solver. The parameters of the problem are mentioned following. The angle of
wedge a is chosen as 1;, gr and 1% respectively and the angle of incident field am
is zero in all three cases. The relative error is computed by L2 norm of error of
E—field at 100 x 100 sample points which are uniformly distributed along x and y
direction separately. Figure 2.23 shows the relative error Of E field by GFEM and
FEM solver with respect to the number of unknowns when the angle of wedge is
changing. In GFEM, the number of modes in each patch is fixed to 6. We can see
that the result of GFEM is much better than FEM. Figure 2.24, Figure 2.25 are
error graph of FEM and GFEM. Fr0m Figure 2.24, The error in FEM around the

46

 

 

Figure 2.24. Error plot of FEM solver

tip of the wedge is much larger than the error at other place, at the same time, it
doesn’t happen in GFEM. The number of unknowns for FEM and GFEM is 1307

and 648 respectively. The angle a in both cases are 27%.

47

 

Figure 2.25. Error plot of GFEM solver

CHAPTER 3

VECTOR GENERALIZED FINITE ELEMENT METHOD

3. 1 Introduction

While GFEM has been used for a range of problems [22, 44] , it has not seen adap-
tation to address vector problem in high frequency electromagnetics. While some
papers exist on application of a variation of GFEM to quasi-static problems whose
governing equation is posed in terms of potentials, they do not address a fundamen-
tal difficulty: The boundary condition of the fields on the material interface require
that their tangential component be continuous, their normal component be discon-
tinuous and the normal derivatives of the tangential conditions be discontinuous as
well. Achieving all of this using functions that are continuous [across boundaries is
almost impossible. Aside from this difficulty, it is also important for the divergence
of the fields be zero (either in a strong or a weak sense). The latter was the principal
motivating factor behind the development of vector basis functions [1] in classical
FEM. As an aside, we note that this is not the first attempt at extending GFEM
to analyze vector electromagnetics problems. That honor belongs to Tsukerman
[32—37]. His treatment of the vector problem is largely inspired by that used in
classical FEM. In Chapter V of [35], vector basis function is defined as thb where
1b is barycentric coordinates, e. g. nodal basis ftmction in FEM and shape function
in GFEM, and v is polynomials. Though in FEM framework, this basis function
coincide with the interpolatory high order edge element, in GFEM solver, Vrb equals
zero in the region where there is only one patch as is this basis ftmction. Except
for this problem, he created basis functions that have zero divergence. But they
cannot handle discontinuities in the constitutive parameters. Furthermore, he did
not demonstrate the convergence characteristics of his scheme and only solved one

eigenvalue problem by one patch where the shape of patch coincides with the shape

49

of waveguide.

Another problem that bothers GFEM, as it also does any higher order solver
is the condition number of the resulting system. Indeed, this is more of a problem
here as opposed to classical FEM as the space of basis functions in GFEM is not
interpolatory.

The objective of this chapter is to extend GFEM to address vector electromag-
netics problems in both two and three dimensions. The principal contribution of

this paper is three-fold:

1. We will develop two classes of vector basis functions for homogeneous domain

and prove bounds on their convergence and the error bounds of linear operator.

2. We will develop basis functions that will enable the satisfaction of bound-
ary conditions across interfaces. These conditions will be on the tangential
components of the fields and their normal derivatives, and the normal compo-
nent of the fields. In this chapter we restrict ourselves to planar and curved

boundaries.

3. We will demonstrate h and p convergence of the proposed scheme for both
Dirichlet and Neumann boundary conditions in two and three dimensional
problems. We shall also demonstrate the accuracy of this scheme for computing

eigenvalue for cavities, waveguides, and partially filled cavities.

4. We will demonstrate the inclusion of analytic behavior in the basis function
space and the resulting improvements in convergence of fields. Our convergence

studies will be restricted to near-field calculations.

5. We will propose a SVD-based preconditioner to reduced the condition number
of the resulting system. The cost and error bound of the preconditioner will
be derived. The reduced condition number and maintained accuracy will be

demonstrated.

50

This chapter proceeds along the following lines; in the next three sections, we
will introduce the construction of vector GF EM solver for homogeneous domain,
piecewise homogeneous domain with planar interface and piecewise homogeneous
with curved interface respectively. The accuracy and convergence of the GFEM
solver will be demonstrated. In the last section, the SVD-based preconditioner will

be introduced.

3.2 GFEM solver for homogeneous domain

3.2.1 Statement of the Problem

Consider a domain 9 whose boundary is denoted by BS2 :2 I‘ = Uz- Fi- It is assumed
that this domain is embedded in a domain fie and its exterior boundary Fe is defined
as Fe := fie n O. Interior to the domain 9, the function u(r) satisfies
(v x [——1—VX] — 15270)) u(r) = f(r)
u(r)
8,: {u(r)} = g;(r) for r E F,-

(3.1)

In these equation 11 is used to denote either the electric or the magnetic field, a(r)
and 7(r) are position dependent constants, B,- are differential operators, and gi(r)
is the ftmction that is imposed on Pi- In the above equation, it is assumed that
r 6 Rd where d = 2, 3. Solution to this problem for any imposed field f (r) can be
obtained using classical vector FEM. This scheme requires tesselation of the under-
lying domain, and specifying an approximation space [1]. The approximation space
is such that it satisfies the de—Rham map, and are based on a space of polynomials.
The question that we intend posing is the development of a possible technique that
permits development of scheme that does not rely on classical tesselation, permits
the use of different kinds of basis functions in different regions of the domain 93.
As directly implementation of the scalar GFEM solver will results in the spurious

problem, spurious-free vector GFEM solver need to be developed.

51

 

 

 

 

 

Figure 3.1. Patch cross the material interface

3.2.2 Basis functions

First, we examine the constraints that we need to impose on vector basis function.
These condition emerge from Maxwell’s equations and they are: (i) in a source free
region, the field is divergence free, and (ii) at the interface between two media,
the tangential component of the field is continuous and the normal component is
discontinuous, and (iii) the normal derivatives of the tangential components are
discontinuous as well. As we shall demonstrate in the Section 3.3.3, satisfaction of
the last constraint is crucial achieving h— and p— convergence. Next, we will prove

those constraints for vector function space.

Theorem 3.1 Figure 3.1 shows a patch with a material interface cross it. The
constitutive parameters are 61 and 62 which are different from each other. Tangential
and normal directions on the interface are defined as t and it respectively. f1 and f2

are vector local approximation in each region. The criteria of the local approximation

52

function can be described as followings:

flt = f2t
fln 24 f2n
61,13, 74 6,1,9“ 1: 1,2,3...

0th.. 7% 65.121. I: 2.3...
V - f1: 0

(3.2)

V-f2=0

Proof: If electric field in both homogeneous domain is denoted by El and E2

then
V- E1 = 0
V-E2=0 (3-3)
E1t=E2t

The above properties of the field require the following properties of the vector local

approximation

f1t=f2t
V'f1=0 (3.4)
V- f2 = 0

From the condition

613111 = 6219721; (3 5)

€1ré€2

It’s easy to derive the criteria

hr; 75 f2n (3-6)

We also can obtain the criteria

(9nf1n = V'fl — Vt ' fit
= V . f2 — Vt - {2, (3-7)

= 371on

While the above criteria is important to construct high order vector local approxi-
mation, it’s not independent.

When proving the following criteria, we can provide one simple example. As
the field has the following discontinuities, the local approximation should share the

same properties.

(95,11, 7e 6,215, l: 1, 2,3...

(3.8)
6.213.. aé 65.12.. l= 2.3..-

In what follows, we start with defining basis functions in homogeneous domains
and then proceed to extend this for domains that contain inhomogeneities. We
shall also derive error bounds when using this span of functions. It is well known
that a vector function can be expressed using a set of three wave functions that are
typically denoted using L(r), M(r) and N(r) [57, 62]. These function are expressed
using L<r> = v was], M(r) = v x [amt-)1, and N(r) = I/nv x v x [we].
where (M(r) for q = l, m, n satisfy the Helmholz equation with wave number n, and
c is a pilot vector. In a source free region, functions M(r) and N(r) are sufficient
to represent the fields everywhere in the region. These fact sets the stage for the
development of vector basis function for GFEM.

Assume that in any volumina ¢q(r) z oap(r) = span {Vi’q} for i = 1, n

where the subscript ap denotes an approximation to the continuous function and

V,- : 2,3(r)vz-(r). It follows that one can envision, two different basis fimction:
113110 e V1 = span {¢(r)V X [62) gm], /)(r)V X V X [5222 Ki} and u2p E V2 =
span {V X [aunyéjfl ,V X V X [6v ’)v(r if.” From the definition of thase ap-
proximations, it follows that for the approximation to be of order p, the highest
order of the function of ’11,: should be of order p + 1 whereas vhf should be of order
p + 2. In this paper, these functions are chosen depending on the problem being
analyzed. The definition of these basis functions also dictates that 1131, has zero
divergence where as ucllp does not. To visualize these functions better, we plot two

h,p

dimensional local approximation function V X [6112. m]with c along the z direction

(or along the direction of invariance). These plots are displayed in the Figure 3.2,
and as is evident from these figures, different modes approximate curls to different
degrees.

When evaluating the terms V X u and V ' u for the first basis function, the

following formula is implemented.

 

 

 

 

 

V X 117-1 =V X vn
z ZQkECi Wk z
W- W-
=V( z ) X v" + z VX vn
ZQkECz' Wk ZQkECz' Wk 3 9
_VW Zokeo, Wk- ”’2' Zzinkec VWk X v“ ( ' )
(29,60, W1.) 2
W2-
+ V X v?
Zakeo Wk

 

 

-- -- no -- — ------ - —
h — .- - — — — - - - — - - — - h
i- - - C C ‘- - - C - - - - - — P
h o c c a - - a u - I - O O c h
4 - o n - o n u n o - - a — n d
* - - - - ‘ - - - - - - - - - G
.4 .- c- c- .- — - a. — - - - - — — all

 

.4——.——-¢———.——--——-1
—————o——o——--——--——

 

(b) Homogeneous O2 mode

1':

u.—.—.—.-.—*-—._——oc———h-

 

_————.—.—-——’———-b
~——.——o—o—.———o-—--———
—-—---~-—--—---b
—-.———.—-—.——.——--—.—--r—-
~——‘-.---—--——--_
—-—-.-—-.---.-._._.-.-i-.
_——.—.—--—.—.-—-—--—h—
~-—o——-—-—o-’--———b
~_‘~““~“‘-“--b
n—.—~—~——.—o—o~—o—o—o—.L.
H“-’-.--~—.—--.-‘-‘_

 

 

p———.-—---—————-
I-o—o‘——.—-——.-.-—-——-c.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(h) Homogeneous 22 mode

56

 

 

 

 

 

 

(a) Homogeneous 01 mode

 

 

 

-~__—P

e e «lfi4/
\\\\\\~—.a”’/’1 d IOIIIOIOII'I."||IOIIIOIOII' d \‘u. ........ .4!
\\\\\s._...”//... 0 ulllltl IIIII III: 0 ... ............ .
\\\\\--—valflll m I IIIIIIIIIIII null m ...........
\\\\s\....;;;.rl1 1 1 llllllllllllll v 0 ...............
¥\\\\~....oollll 1 r llllllllllllll 2 y ......... .

S m .....
I |||||||||||| I'D. u 1 oooooooooooooo n. .......
III|II III'I w .............. m -. ........... -
I’ll. Ill'l n .............. ... .-
I IIIIIIIIIIIIII L e uuuuuuuuuuuuuu m .. ....... -
Illlla....~\\\\\ MD A nnnnnnnnnnnnnn .. mo .. ............. A
Illalo...-\\\\\ n lllllllllllllll 1 m ................

I
ll/a:.....~\\\\\ 0 1 llllllllllllll r .0 . ........... .
I’lllro.-\\\\\\ H ffffff O- fffffff l H A: ............ \L

"""""""""" [’4 III-a-
I/rzip....a\\\\\ \l to. .. y.... ..\\
f PL: 4

zaa.....w7w~ww d ( ..-

e L-prn liq—dd e bkrhhhhfhrrrirrlr

r

TOIOIIIIIDICIIIOIDIOIDIIIDI.“ IN “- — _ v . . a — a r m j\\\ \ ..... r I If’!
I'IIIIIOIVIDIIIIIIIIIIIOIIIIIIII V~ - . . . . . . p a r 1.: I \ \ o 5 III.
Y."IOI.ICIOIOII-IIIIIIIIIIT m ‘- - . ...... . . — v r m .I I I I o I I I l
IIIIIIIIIIIIIIIOIIIIIIIOI.I.IIII 0 <\ o . ........ . c p r 2 fl - s n s L
II'IIIOIIIIIIIIIOIIIIIIIIIII' 1 ‘ s ........... . 5 r 1 v .............. A
10‘1““‘I'IOIOIOIO‘IOI. m ............... v w ..............
TIIIOIIIIIIIOIO‘IOIOIOIOI." m ............... v m ..............
T‘I‘““II‘I|I.III|IQIO‘ v .........................
I"'|'ll""“‘ m I O IIIIIIIIIIIII m llllllllllllll
‘IOIOIOI‘I.I.II-I.I‘II|IIII % v r o .......... . ~ ‘ W0 ............... fl
TICIIIO‘IO‘“IOIOIOIOIOIOII m ro . . . ...... . . ~ .1 m 1 o I I I !
I‘I‘IIOIOIOIOIOI'IOIOIOIOYO O rv — . . ...... . . - ‘ < O l I I u ......
10“““10‘ “““ v. H rra. . . . . - u - H 1 I I I 1 . u I I I I fit
IIIIIIII I'IOIOI.‘ o.-“ l'fllo a ‘\\.\‘

) fr.... \I/

Ph-—~ e ryrfrr Illbh
C a—«q I
( (

 

(g) Homogeneous 21 mode
Figure 3.2. Vector plots of homogeneous local approximation function V,-

 

Wz' n
v.
ZQkEC' Wk 2
w,- w-
V- v? + V( 2 ) -v"-2
ziszEC Wk: ZQkECz' Wk 1. (3°10)
VWz' anec, Wk - Wz' anec, VWk n

= 'V'

(ZQkECZ' Wk)2 z

V-u?=V—(

 

 

 

In the above formula, v? denotes vector local approximation from
the first kind of basis function 11,111, and belong to the function space
span {V X [62) :1 m] V X V X [621; 515]}. The above formula indicates that
only v? and V X V,- are necessary to implement the bilinear form of the equation.
Next, we derive error bounds for each of these approximating span of functions as

applied to vector electromagnetic problems.

3.2.3 Error bound

3.2.3.0.] Error bounds of basis functions: The global accuracy of the finite
element method depends on the local accuracy of the approximation to the solution.
Thus, predicting the local error bound yields the global error bound of the solver.

In what follows, we derive error bounds for both types of approximations presented

earlier: ucllp and ugp

Theorem 3.2 If

“11), ”L00 30 00(9 )

”W, c.

”Loom S diammé‘)

llVXllem-vz'f’ 516)] ,,<

“V X V X [($71 — ’02 :6) c]l|L2(QflQ.:1§E 62(‘i,h,p)

EIMEN VrEQ card{i|rEQZ-}SM

        

57

Then the error is bounded by

1/2

lL2 <\/2_A/—ICOO 261(iflhp)+262(i,,)hp (3.12)

 

Hu- nap
where Coo and CV are constant. 61 and 62 are bounds of error of local approximating

function in each patch.

Proof:

1
u — “GPHL2(Q)

= [V X (957,150+ V X V X (95715)] — 27% [V X (”5325) + V x V X (02:96)]

i

= 2.. [w ({w .m_ it cum [vxvx<{¢n—vzre}>l

2

 

 

2

 

 

 

 

142(9)

$2

 

 

swam vane»

 

 

L262)
:21),- [V X V X “(fin “ ”if: 6)]

. L262)
<2M: «1».- [V X ({¢m ”hp} 0)] Mammy)

+21%:
i

52MC§O E H (i, h, p) + 53(1, an]

i

2
+2

 

 

 

 

 

 

 

 

2” [V x V X ({‘m‘ if: 6)] “L2(nno,)

(3.13)

58

Theorem 3.3 If

”V? ”L00 5 000(9)
”Wt Loom) 5 diaZYQh)

hp
llrm-

 

 

 

 

 

L2(909i)< _

”V x [(qu— vi hm)c C]‘|L2:an;< ) — 62(i,h,p)
3M 6 N Vr E Q card{i Ir 6 0,} g M

then the error is bounded by

2 1/2

CV
<¢§T 202%(zhp)+2m1( 2.,hp) (3.15)

 

lf“'“3PlL2

where Coo and CV are constant, and 61 and 62 are error of local approximating ~

 

 

 

 

 

 

functions .
Proof:
_ 2 2
2
h .
= w [mum-u or]
i 112(5))

2

 

 

 

= Z [riv X (W -v§f;fi}é) - (W ”233} ) X W4

3

52M: *in X ((9% ‘ ”23} ‘3) ll:2(nnfl,)

+2MZH({¢m— riff; er)sz,,-2

Lfinnnu

 

L2(o) (3.16)

 

 

 

 

 

<2M 030 e i,h e (i,h
2“ p) +diam(2VQh)2:-:1( 17)]

2

59

From the above equations, it appears that the space 113p does not converge with
decreasing h. However, it should be noted that 61 scales proportional to h. This
implies that the rate convergence of this approximation w.r.t. both h and p is similar
to that of u(llp. Note, using the above proof it can be shown that the basis functions
proposed by Tsukerman will exhibit convergence characteristics similar to 1131,. The
principal deficiency of ugp is the fact that it cannot be easily extended to satisfy
boundary conditions of fields across material interfaces whereas 11‘1”, can be readily

modified.

3.2.3.0.2 Error bound of bilinear Operator: Next, we will show the error
bound of the bilinear operator. As preliminary, we will assume the error bound of

the following terms. If

I
”V‘bflhwm) S $52—12)

VX “(hm—"z m)c ‘32-]IIL2(QnQ)S 61(23hm)
V x v x [(05” will ’71:) a] 2(iahm) (3-17)

 

 

22>?

 

L00 5 000(9)

”L2(Qflflz)< —

VXVX l(¢m"”zhm)c (Ea-Sclllflmno) 3(i’h’p)

VXVXVX[(¢n—zi,’:)c] €403,147)

 

 

”L2(Qflflz)< —
EJMEN VrEQ card{i|rEQi} SM

60

Then the error of V x ucllp is bounded by

“V x u — V x ucllpllL2

g \ffifl; {5%(2', h,p) + 6303 Ml} 0520 (3.18)

2 1
CV 2

+ Z {63(i, h,p) + 63(i,hap)} M]

where Coo and CV are constant. 61, and 62 are bounds of error of local approxima-
tion flmction in each patch. 53, and 64 are error bounds of curl of local approximation
function in each patch. If the order of scalar function (pm and on are p+ l and p+2,

then we can approximate the function 61, 62, E3 and 64 as

61(23hap) ~0(hp)

6203 h, 10) ~0(hp)

63(i,h,p) wow-1) (3-19)
640', h,p) wow—1)

diammf) ~O(h)

then

lleu—quép

 

 

L2 ~ 0(hP—1) (3-20)

61

Proof:

“V x u — V x uclw 312(9)

=uv x [Z «a (V X ({2m - 2325159)]
+v x $101 (V X V X ({CD" ”hpc D] llL261)
V x [2102' (V X ({WL ‘ ”£371; (2)]

_ 2 L262)

.7. [;i.(v.vx ({M—vz ram]

2

 

 

_<_2

 

 

 

 

 

2
+2

 

 

 

 

L262)

S4

 

 

szxvx<{M—vn}e)>

 

L262)
2

+4 ZVX({¢m— 122,2} C)><V¢z'

 

 

132(9)

+4 Zoz- (VxVxVx ({an—vzlff}a))
i

2

2

 

 

132(2)

 

 

+4 ZVXVX ({on—vzf}é) XVTL'z'
'i
54M: 952' (V x V X (”m—”2% 0)) l:2(onn-)
+4MZHVX ({95m— vh’p} c)wa 2

LWQHQ)
+4MZHW (V X V x V X W" ”zhp} ))“L2(909)

 

 

L26?)

 

 

 

 

 

 

+41W:ZHVXVX({¢n- if: C)XVUi2

L2(Qflfli)
. . 2 2 2 __.__02
34M 2: {ega h,p) + 64(z,h,p)} Coo + Z {61(i,h.,P) + 62(i,h=p)} _ V h 2
i i diammi)
(3.21)

62

If we also have

”V X villL2(QnQi) S 65(i9hap)

(3.22)
. < a h,
llvzllL2(QflQz-) — 56(3. P)
Then the error of the linear operator is bounded by
1 1 2
”(V x uap)-(V xv)+uap-v— (V x u) - (V xv) —u-v”L2(m
1/2 (3.23)
S {8M 2: f1(i,h,p)65(i,h,19) + 4MC'go Z f2(i, h,P)€6(ia M10}
2' i
where
2 2 2 2 2 0%
f1“, h,p) = {6303 hp) + 640': MM} 000 + {61(iahaP) + 62(3) h,p)} m
i
f2(i, h,p) =e%(z', h,p) + 63a h,p)
(3.24)
If we assume that
€50} h,p) ~0(1/ h) (3 25)

€603,210) ~00)

Then

N(qucllp)-(Vxv)+u31.p°V—(Vxu)'(VXV)—U'V||L2(m

1/2
3 {8M 2 f1 (2‘, h, p)e5(i, h, p) + mag, 23 f2(i, h, mafia, h, 19)} (326)

Z 2

N 001094))

63

Proof:

“(VX“ital-(VXVil+uc1rp'Vi—(VXU)'(VXVil-U'Vi

5 H (v x {“‘l‘P ‘ “D ° W x Vi) _ (“‘11P _ “) William

32

+2

52

 

 

 

 

 

 

(116111, — u) - (V x Vi) :26»
(“gm " ‘1) 'Vi 22(9)

V X {uclzp - u}l|:2(9) “V X Villi2m)

2
1 _ . 2
nap ullL2(Q) llvzllL2(Q)

_<_8M Z {ega h, p) + 63a, h, p)} 030
i

02
+ Z {€¥(i, h,p) + 63(2', h,p)} —V—2 65(23 ’W)
i

+

diamm?)

4MC§O Z [5%, h, p) + 63(2', h,p)] £60312?)
i

 

 

2
L26?)

(3.27)

3.2.4 Imposition of boundary conditions and the resulting bilinear form

As is evident from the preceding discussion, the basis functions are not interpolatory,

nor do the domains conform to the boundary. While this fact does not affect the

imposition of the Neumann boundary condition, imposition of the Dirichlet bound-

ary conditions are not as simple as those for classical FEM. Methods that have been

used to impose this condition included augmenting GFEM with classical methods

at the boundary or using a Lagrange multiplier. Both these method are not stable

as they do not satisfy the Babuska-Brezzi condition [52—54]. It can be shown that

stable results may be obtained using either stabilized Lagrange multiplier technique

or the Nitsche’s method. We have used the latter with success for solving the scalar

64

Helmholtz equation as have others [55, 63]. In this paper, this technique is extended

to vector electromagnetic fields. Thus, given (3.1) with

u(r) X it = gd(r) for r E Pd
fl X (V X u(r)) = gn1(r) for r E Fn (3.28)

u(r) - it = gng (r) for r 6 Tu

where I‘ = I‘d U I}, we need to find an approximation uap(r) E Vn C H (Vx,§2).

It can be shown that the bilinear form for solving this problem may be written as

Ah(u,w) = fh(w) Vw E Vn

1 2
.A u,w =/dQ—qu-wa—w/dfl ru-w
h( ) Q 00') Q 7()

+fi1/Pddl‘(uxn)o(wxn)+a/Pddl‘(V'u)(V'w)

1 , l .
+/Pdd1'~m(uxn).(vxw)+/PddI‘—J(—IT)—(wxn).(qu)
+fi2/F dI‘(u-n)(w-n)

fh(w)=/Qde-w+/Fdd1‘ 8%?)gd-(wa)+,81/Pd dI‘gd-(wxn)

1
+ dP—g -w+,6 / dI‘ w-n
An 01(1‘) n1 2 Pn 9n2( )

(3.29)

In the above equation, ,81 is a constant that imposes the Dirichlet boundary condition
using the Nitsche’s method and makes the overall system coercive. The means
estimate the value of ,81 follows those derived by [50]. The quantities associated
with a and ,82 are regularization terms. We have found that these are not necessary
when using approximation functions of the type 113,, but make the convergence

curves smoother when used with approximation functions of the type ucllp.

3.2.5 Numerical experiments

In what follows, we shall present a series of numerical experiments that will serve
to demonstrate the accuracy and convergence of the vector GFEM presented in this
chapter. At first will focus on convergence issues related to problems in homogeneous
domains. Then we shall demonstrate the versatility and convergence properties
of this technique by studying scattering from a wedge using non-polynomial basis
functions. In some of these numerical experiments we will demonstrate the relative
convergence of both nip and ugp.

First, consider a homogeneous rectangular computation domain, where the solu-

tion to the vector wave equation is assumed to be
u(x,y) = 1.32 cos(1.76:r) cos(1.32y)§: + 1.76 sin(1.76:c) sin(1.32y)g} (3.30)

in the domain 52 = (0,1)2 with boundary I‘U21I‘z- where H] : r E O x (0,1)},
{1‘2er (0,1) X 1}, {P3 : r E 1 X (0,1)}, {F436 (0,1) X 0}.

This field is imposed on I‘ as either a Neumann or Dirichlet boundary condition.
The problem is then posed as (3.1). Given the geometry of the problem and the
nature of solution, the approximation to the field in Q can be written as u},p =
span {1,1},- (r)V x [(311251] }. In this experiment, all nodes are distributed uniformly,
and rectangular patches are used. The size of each patch is chosen to be 1.5 times
the distance between the nodes. Unless otherwise mentioned, all examples in this
chapter are built in a similar manner. The weight functions Wi(r) is a tensor
product rooftops, and the local approximation function 1},-M(r) and vim (r) are tensor
product of Legendre polynomials. Figure 3.3(a) shows the h, p convergence of GFEM
when the Neumann boundary value problem is solved. Figure 3.3(b) shows h, p
convergence for the Dirichlet boundary value problem. In the figures, the L2 norm

of relative error is evaluated at uniformly distributed 100 x 100 sample points. As

is apparent, the technique shows excellent convergence.

66

Next, we solve the same problem for both boundary conditions. However, the
approximation space will be 11,211, = span {V x [api(r)v£l#] }.All other parameters
(size of the meshes, distribution of the nodes) are identical to those used earlier. As

is evident from Figure 3.4, the error decrease with decrease in h and increase in p.

67

ErrorplotforPDEMhNeumannbomdaycondib‘on

v '7 1

 

 

—-—-P-1 ‘
——P-2 .

 

 

 

- ‘AAA-A‘ A n .-

L2 error

 

 

 

 

 

 

 

 

 

10'
h-Distance between neighbour nodes
(8.) Neumann
ErrorplotforPDEwith Diridiletbounderycondition

—— P-2 .

—— P83 1

// 1:

lg 1

0 i

N

-' 1.

1o—1

h-Distance between neighbour nodes

(b) Dirichlet

Figure 3.3. Convergence graphs for vector Neumann and Dirichlet BVPs by ucllp

0 Error plot for PDE with Neumann boundary condition

 

 

—-—P=1 ‘
—-°-—-P=2 -
—°—"P=3 1

 

 

 

A
0|
.vwvnq ' v

.5
O
I
"'1
..r

v 1'17

1 viii

L2 error

.5
O
“.1 I 2--...4

 

 

 

10'
h-Distance between neighbour nodes

(a) Neumann

Error plot for PDE with Dirichlet boundary condition

 

 

——-P=1 3
—-—P=2 .

 

 

 

L2 error

 

 

h-Distance between neighbour nodes

(b) Dirichlet

Figure 3.4. Convergence graphs for vector Neumami and Dirichlet BVPs by ugp

69

Next, we examine the accuracy of three dimensional vector basis function in
homogeneous computation domain. The computational domain is Q = (0,1)3 and

it is also assumed that the field behavior in the domain is governed by

u(r, y) = —19 * cos(0.94x) sin(1.26y) sin(1.57z):i:
—17 * sin(0.94.r) cos(1.26y) sin(1.57z)3’) (331)

+25 * sin(0.94:r) sin(1.26y) cos(1.57z)2

It is apparent from the above equations that both functions M(r) and N (r) are
necessary to represent the field. Consequently, the function space chosen is ucllp =
span {¢i(r)V x [évi’m] ,z/Jz-(r)V x V x [612237,] }. Again, all the parameters used
in the computation are the same as those in 2D problems and are only extended to
3D. The L2 norm of relative error is evaluated at 5 x 5 x 5 uniformly distributed
sample points. Figure 3.5 shows the h and p convergence of those vector basis.
The above results demonstrate h and p convergence, but the anticipated rate
of convergence th is not achieved. This said, it is worthwhile to note that the
theoretical convergence rate is obtained from dispersion analysis and is appropriate
for analyzing convergence of far fields as in RCS. But, the examples shown involve
near field quantities. To illustrate convergence, we again consider a homogeneous
domain (2 = (0,1)3. It is assumed that the field behavior in the domain is given
by u(sc,y) = exp{—j0.9425z}g}. All other parameters, viz., distribution of nodes,
overlap, etc, are identical to those used earlier and as before the L2 error norm is
evaluated at 5 x 5 x 5 uniformly distributed sample points. The convergence data

is shown in Figure 3.6, and as is evident from this figure, the convergence for p = 2

follows the predicted behavior. However, the curves for p = 1 are a shade worse.

70

hp convergence for Neunenn BCs

 

 

 

 

 

100 Y I I t I I I ‘
—-—P-1 i
—-—-p-2 1
10'1 -

//

L2 error

/

’I/

10' " ':

 

l

10-0.8 10-0] 10-0'6 104.5 10-0'4 10-03 10-02

h-Dietence between neighbour nodes

Figure 3.5. Convergence graphs for 3D homogeneous vector N eumann BVPs

 

 

—9—P=2

A

 

 

 

 

A :1
1 0
h-Distence between neighbour nodesa)

 

Figure 3.6. Convergence graphs for 3D homogeneous vector Neumann BVPs (Plane
Wave)

71

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . p = 1 p = 2 p = 3
ifgfiilpal h = 1/4 0.0616 2.88e—4 2.82e—6
TE h = 1/6 0.0297 6.178—5 2.42e—7
First 8 h = 1/4 0.1629 3.99e-3 1.65e-4
modes h = 1 / 6 0.0743 8.74e-4 1.21e-5
Principal h. = 1/4 0.0182 2.756—4 1.14e—6
TM Mode h =1/6 9.93e—3 4.506-5 1.08e-7
First 8 h = 1/4 0.1674 0.0349 5.83e-4
Modes h = 1/6 0.0992 5.276—3 1.12e—4

 

Table 3.1. Error in eigenvalues computed using ucllp

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

. . p = 1 p = 2 p = 3
333:1’31 h = 1/4 0.0206 1.17e—4 1.3297
TE h = 1/6 8.888-3 3.766—5 4.52e-9
First 8 h = 1/4 0.0697 1.3le—3 1.07e—5
modes h = 1/6 0.0263 3.81e-4 7.06e-7
Principal h = 1/4 7.13e—3 1.62e4 6.82e-8
TM Mode h = 1/6 3.76e-3 3.66e—5 7.04e-9
First 8 h = 1/4 0.1878 5.82e—3 3.97e-5
Modes h = 176 0.0420 1.89e-3 1.14e-5

 

 

 

Table 3.2. Error in eigenvalues computed using ugp

Next, we compute the eigenmodes in a rectangular waveguide. The dimensions
of the waveguide are chosen to be 1m x 1m. The convergence data (error in L2
norm) for both TEz and TMz modes are presented in the Table 3.1 and Table 3.2.
The basis functions used in these experiments are for both ucllp and 1131,. As is
evident from this data, the results obtained converge rapidly with increasing order

and refinement. Furthermore, results obtained using basis functions of type ugp

converge faster than those obtained using ucllp.
Until now, the basis functions used in all the examples have been based on

tensor products of polynomials. Here, we recomputer the scattering from a wedge

72

 

 

 

 

 

Figure 3.7. Geometric description of the wedge

in the last chapter. The problem can be solved as vector problem when magnetic
field H(r) is computed. In what follows, we will explicitly incorporate the known
functional behavior. Consider a two dimensional wedge with angle air that is located
in a rectangular computational domain (2 = (—1.0,0.5) x (—0.7, 0.3). The tip of
the wedge coincides with the origin of the computational domain and one side lies
in the y = 0 plane. A field that is polarized along the 2 and propagating along
It = —0.9659:i: — 0.25883? is incident upon the wedge, and the wavenumber k0 = 3;
this geometric configuration is shown in Figure 3.7. It can be easily shown that the

fields are given by

H(r) =V x E(r)
47rE0

 

 

__ 'U ' . _ ' ' _
E. 37, _ a g] Jv(kp)8mlv(am + r a)]sm[v(<o + M 0)] (3.32)
'0: mr m=1,2,---
27r — a

To solve this problem using GFEM, we use Jv(kp)sin[v(¢ + 7r — (1)] to create
the function space vi,m- The infinite extent of the wedge is simulated by imposing

analytical solutions for the infinite wedge as Dirichlet boundary conditions. The

73

relative error is computed using L2 norm of error in the electric field at 100 x 100
sample points that are uniformly distributed along :2: and 3) direction. Figure 3.8,
Figure 3.9 shows the analytical solution and the graph of error when the number of
patches is 12 x 9 and the number of modes in one patch is 7. From this figure, it
is apparent that this scheme captures the field behavior extremely well at all points
in the computational domain, and as is to be expected, the results are excellent at
the tip. Graph of error shows that the overall error in the domain is dominated
by that on the boundary. Finally, Figure 3.10 shows the relative error with respect
to the number of unknowns, where m is the number of modes used for the local

approximation in one patch.

74

 

(b) Error in the a“: component

Figure 3.8. Analytical result for scattering from a wedge and error plot for y com-
ponnet

75

 

(b) Error in the 6 component

Figure 3.9. Analytical result for scattering from a wedge and error plot for y com-
ponent

76

1O

 

10’

 

 

Am A AAAAJA

..--l

 

 

6
5
E 10.2: 1
O ’ :
C 1
N t
-’ m=6 ‘
10.3: 1
i m: 1
10'41 ”“2 we:
10 1o 10
Number of unknowns
Figure 3.10. h, p convergence for wedge problem

77

3.3 GFEM solver for piecewise homogeneous domain with flat interface

3.3.1 Basis function

The basis functions defined earlier are continuous as are its derivatives. While these
basis functions are ideal for modeling fields in a homogeneous domains, it is appar-
ent that these are far from satisfactory when the domain comprises of a material
discontinuity. Ideally, the approximation space should consist of functions whose
normal components are discontinuous and whose normal derivatives of the tangen-
tial components are discontinuous as well. To address this problem, consider a
domain Q,- = Q: U Qi— and I"? = 652:" F1 697;— . The dielectric properties in 93: are
assumed to be different as shown in Figure 3.11. In what follows, I‘ is assumed to
be planar, the unit normal vector to this boundary is denoted by C, and the bound-
ary is defined by the plane I‘ = (03;. Assume that the local approximation func-
tions, {V x (Cvg’mp), V x V x ((02 5)}, are augmented by two functions f+(r)
and r; (r) with supp {rim} = of and supp {ri (1)} = a; . Next, we derive the
characteristics of these functions Vr E Fit.

Prior to doing so, we define the following: r = p + (i, V = Vp + 8C6 and

f§'=(r) 2p fa+ (fl-CE for a = :l:. Operations on these functions are defined as:

Tar) = of (rage) — fi;(r)) Vk = 1,2,...
N0(r) = fggc) — 1,; (r) (333)
Nk(1‘) = of (fit-(r) — f2.2(r)) Vic = 2,3, - .-

Using these definitions, the conditions that the local approximation functions fi

78

 

 

 

 

bi

 

 

 

Figure 3.11. Interface of two media

(and the fields) must satisfy are

T0(r)lVreI‘;t == 0 (3.34a)
Tk(r)|VreI‘:.t # O k = 1,2, - .. ,p —1 (3.34b)
N0(r)lVreI‘;.t ¢ 0 (3.34c)
Nk(r)|vrerii # 0 k = 2,3, . u ,p — 1 (3.34d)

V . ff (r) = 0 (3.34e)

Note, the partition of unity 212,- is continuous across the interface of two media.
Therefore if the local approximation functions satisfy the boundary requirements
exactly, then the basis functions will also satisfy those requirements exactly.

Next, the additional local approximation function may now be defined in a man-

79

ner similar to that used earlier.

f,-+(r) =V >< (v3£(p,é)€) + V x V x (1133300, ()5)
=Vp >< (vflmxfi) + Vp x Vp x (v33(p,C)C) 1' E 923
+6 x V, x (€6,123, 31p. (1)
f,- 01- — vp x (1133(0C0X) + V, x vp x [C(c— €016,213, 300. CO1]

+C x Vp x ((6411,: (p, C0)) r E Q;
(3.35)

It can be verified from the above definition of local approximation functions that
. h, “

the functions Vp x Vp x (112.;(p, C)C) and Vp x Vp x [C (C— C0)6C11£;l p(p, C0)] are

normal to the interface, and functions Vp x (112’ m(p,C )C ), Vp x (112- ;n(p,CO)C),

C x Vp x ((6,137, p(p, C)) and C x Vp x (Cacvgfn (p, (0)) are tangential to the

th

interface. It also follows that to obtain a p order approximation, the functions

hp

112-” m pand 112 7fShould be of order p + 1 and p + 2, respectively. Next, we use these

definitions of nfz-ifi) to verify that they satisfy the conditions (3.34). Substituting

these functions into (3.33) and evaluating these Vr 6 F23, one obtains

Tk<p.<o1= [631,+ (r1 — 631,111] |C=C0
=vpx (<63v3,3,(p.101)+éxvp >< (Ca3+1v 3,300.10 1)
N0(p1Co) = [13,01 — 1,; (r1] | C: ,0
=vp x V, x (113311. (016) -c‘ (336)
Nap, <01 = [831,32 - 831,2] I,=CO
v, x v, x (gagegjfwm) 1; > 1
0 k = 1

It is apparent from the above equation that the functions Tk(r), NO(r) and N 110')

are all not equal to zero. Thus, the additional local approximation functions defined

80

in (3.35) has discontinuities in (i) its normal component and (ii) normal derivatives

of its normal and tangential components. Also,

v 131 =v - v x (v3;3,(p,C)é) +V - V x V x (133(1), ()6)
1 ’ (3.37a)
=0

and

v - 1; =v - V, x (112%(01Colé) + V - Vp x Vp >< [C(C — (016413311, (0)]
+V - {C X Vp >< (533112301, (0))}
:0 + V, x W x (66,1331), (0)) — C ° (W x Vp >< (Cacv333(p£o))]
=0
_ (3.37b)

Thus, these basis functions, together with the basis functions ucllp presented in Sub-
section 3.2.2, satisfy all the criteria laid out in (3.34). This is pictorially depicted
in Figure 3.12. Here vector plots of f i are displayed for a patch that straddles the
material interface for different orders. From these figures, it is apparent that the
tangential components of the span of approximation functions are continuous across
the interface whereas the normal components of the span of approximation func—
tions are discontinuous. In some figures, the approximation function in Q— is zero,
but that does not effect the discontinuity created at the interface. Consequently,
the basis function space V162) is augmented by the functions 132-(rm.i and denoted
by film). It can be easily verified that the function space {11(9) c H (Vx,o).
The method presented here, defining additional functions, is one way that can be
used to impose appropriate conditions on the fields. It is also conceivable to define
these jumps via other means. For instance, representing discontinuities using surface

integrals along the boundary or appropriately defining the vector potential L.

81

 

 

 

 

 

 

 

 

__.p. ........
_._.. .. . . .
—_—.. .. ... .
——... .... . .
——-. ...... .
—-—-. .u.... .
—--. .o o...
—._-. ....... .
-.... ........
__-.. .. .....
-—... ........
———.. ........
—-_.. ........
—‘_.. . .. .
-_.. ... .....
_-PL

__hb_ __~__h—_
_.__. ___...~_
.____ ___..._.
_.___ ____..._
..._. .______.
-_.__ .__.._..
_..__ __...__.
_._._..__.__.
_____ .__._.._
___...__.____
_._._ _.___~..
.__._ .._____.
_____ ..___...
___...______.
_.___ _____-.
b_VbP_nb—p_~—

 

 

(b) Inhomogeneous 02 mode

(a) Inhomogeneous 01 mode

\NNNNKKN

 

\

\

\

\
\\\\\\h
.\“‘~

"IIII

 

\\\\\\\\\
\\\\\\\\\
\\\\\\\\\
\\\\\\\\\

‘\‘\\\\\‘

u “.l‘\“u|\

I“-l‘-||

Illllllll
I’ll’ll’l
III/Illll
lll/l/ll/
[Ill/Ill!
(I’ll/III

 

III“\\\

fIIIIr

 

 

Ill/AXII

 

 

 

 

 

(d) Inhomogeneous 11 mode

(c) Inhomogeneous 10 mode

 

llllll

 

 

 

 

h I h I l h I P

 

--------

uuuuuuuu

tttttttt

........

oooooooo

llllllll

 

hhbp
‘s....
“...
‘s—..
Vs...
.s...
i .....
r

7.

r.
yo.
yo...
7’...
Z...
7'....
a—_.

 

 

44444441

 

(f) Inhomogeneous 20 mode

(e) Inhomogeneous 12 mode

 

......

 

IIIIIIIII

IIIIIIIII

oooooooo

.......

........

IIIIIIIII

nnnnnnnnn

 

""""

 

 

kkk‘tx‘L

 

 

ru\\\\u‘o\|I-\\It\|\\
v. ||||| ‘ “““““
r IIIII t """"""
T IIIIIIIIIIIIII
lllllllllllllll
.............
.............
.............
x ooooooooooooo A
nnnnnnnnnnnnnnn L
L ..... 1 IIIIIIIII
L ..... L """"" l
I ''''' .1", """ 0.
I’I'l Ill"”"'l
F'V' rr'r'rrr

‘

 

 

 

(h) Inhomogeneous 22 mode

(g) Inhomogeneous 21 mode

Figure 3.12. Vector plots of inhomogeneous local approximation function v,-

82

3.3.2 Proof of the completeness of the vector basis function

In what follows, we present an alternative proof that the proposed basis functions
are complete. Assume that the field in domain 52+ and domain Q" , shown in Figure
3.11, are u+ (r) and u- (r), respectively. Two fields are independent of each other
except for the continuity of the tangential components of u(r) and V x u(r) at the
interface. To prove completeness, we redefine the fields in the domain 9 = (2+ U0“.
To this end, let ut (r) denote the field for r 6 Q, and ur(r) = u+ (r) — ut(r) for
r 6 9"". The fields in 52+ and 9" can be rewritten as
u+(r) = ut(r) + u7~(r) r E 9+

u(r) = (3.38)
u-(r) = ut(r) r E Q—

where fields ut(r) and u7-(r) are assumed to be smooth as are their derivatives.

They are independent of each other, and ur(r) and V x ur (r) have zero tangential

components at the interface. Because u(r) is divergence free, we can express it as:

ut(r) =V X (1720+ V X V X (Até)
(3.39)

ur(r) =V x (FTC) + V x V x (Aré)

where Ft, FT, At and Ar are independent quantities with constraints Fr((0) = 0
and 64A,.(C0) = O. The constraints are from the fact that the tangential compo-
nents of ur (r) and V x ur(r) are zero at the interface. Imposing a requirements,
that any approximation to this function be pth order, the basis functions (or local
approximation functions) should be able to approximate Ft and Fr to p + 1 order
and approximate At and Ar to p + 2 order.

Next, we will demonstrate that this requirement can be satisfied by the basis
functions defined in this paper. Let fo (r), fi(r) be the local approximation func-
tions that represent the field locally. In the patch across the interface, we approx-

imate the fields by approximation function in homogeneous domain and additional

83

approximation function that creates the discontinuity satisfying constraints at the
interface. Then, the functions for homogeneous domain and inhomogeneous domain

is denoted by f0 in whole domain and fi in 0*, respectively. More specifically:

f0(1‘) =V X (Um,0(r)C) + V X V X (1172,0005)

f+(r) =V x (vm,a(r)é) + v x v x (vn,a(r><‘)

f‘ (r) =vp x (vm,a(P, Cox“) + vp x Vp x [(c - (0)64vn,a(p.<o)é] (3.40)
+5 >< Vp >< [act'n,a(P,C0)C]

A

2v x (vm,a(p, (0)0 + V x V x [(C — (0)04vn,a(p,Co)§]

where we omit the superscript h and p and subscript i and add subscript 0 and a
for homogeneous basis and additional basis for brevity. Obviously ”Um,0, vmfl are
p + 1 order and vn,0, 'Un,a are p + 2 order. Let Fr,app, Ft,appv Ar,app and Atflpp
be local approximation to Fr, Ft, Ar and At. Compared with (3.38), the local

approximation can be written as:

Ft,app =vm,0(l‘) + ”Um,a(Pa C0)
Fr,app =’Um,a(r) — ’Um,a(P, CO) (3.41)
At,app =vn,0(r) + (C — (0)3(Un,a(P: C0)

Ar,app =Un,a (r) — (C _ (0)3(vn,a(P, (0)

Obviously, Ftfipp and Fr,app can approximate Ft and F,- to 19 +1 order, and Atflpp
and Ar,app can approximate At and Ar to p + 2 order. Constraints Fr((0) = 0 and
BCATKO) = 0 can also be satisfied. Therefore the basis functions are complete to

order p.
3.3.3 Numerical experiments

In the next series of examples, we consider problems wherein the medium is piecewise

homogeneous. We will examine both boundary value and eigenvalue problems. First,

84

consider a 2D inhomogeneous domain. The computational domain is defined by
Q = 0+ + {2" with (2+ = (0,1) x (0.5,1) and S)— = (0,1) x (0,0.5). The relative
permittivity in (2+ and {2‘ are chosen to be 5r = 1,3, respectively. The fields are as
follows: Vr 6 (2+ the electric field is E = (-:i: — J39)exp{—j [V355 — (3] —1/2)]}, and
‘v’r E Q“ the electric field is E = 1/\/3(—\/3:i:— y“)exp{—j[\/3:r— 3(y—1/2)]}. Using
this data, we derive the apprOpriate Dirichlet and Neumann boundary conditions
that are to be imposed.

The basis functions used belong to V(Q). These basis functions create a dis-
continuity in the normal component and the normal derivatives of the tangential
component. A tensor product of Legendre polynomials is used as a generating func-
tion, the L2 norm of relative error is evaluated at uniformly distributed 100 x 100
sample points. This is depicted in Figure 3.13, and as is readily apparent, the
method shows h, p—convergence when either the Dirichlet or Neumann boundary
condition is used.

Next, we examine the performance of these basis function when applied to a 3D
problem. The computational domain is defined by Q = {2+ +9" with 9+ = (0,1) x
(0, 1)x(0.5,1) and fl- = (0, 1) x (0, 1) x (0, 0.5). The relative permittivity in 52+ and
(2" are chosen to be a,» = 1,3, respectively. The fields imposed are as follows: Vr 6
9+ the electric field is E = g)exp{—j[\/3:r—(z—0.5)]}—3)0.56:rp{—j[J3m+(z—O.5)]},
and Vr e {2— the electric field is E = QO.Sexp{—j\/3[a: — x/Z’Kz — O.5)]}. This data
is used to impose the Dirichlet boundary conditions. Figure 3.14 demonstrates the
h and p convergence. These results are obtained by sampling the fields at 5 x 5 x 5
points in domain and computing the L2 error norm. As is apparent, the results are

excellent.

85

 

1o . - -LL,,

 

—-—p=1s
—e—p=2:
—l-—P=3-
+P=4

10'1

 

 

 

1o-

LZError
—A
o:
L m
........ .........
. . . .......l . .......l . A

 

 

10

”’5: i
10...; .3
10'7: - - . A . - 1 ‘

 

-1
10
h-Distanee between neighbor nodesao)

(a) Neumann

 

 

 

 

 

v v . v I ‘
—'—P-‘|3
+P32.
+P‘4:

I

E 1

0

I:

I.”

3 '

3 .
t
“J

g 1

N .
.J

I.

'i

A A A -_ A ‘ l

 

 

 

-1
10
Distance between neighbor nodesOD)

(b) Dirichlet

Figure 3.13. Convergence graphs for 2D inhomogeneous vector Neumann and Dirich-
let BVPS

 

3 :1//
10'1

§ _
I.th 10 2E
N r
'J l
b L/
10.3? :3

 

 

 

10'
h-Distenoe between neighbor nodesoo)

Figure 3.14. Convergence graphs for 3D inhomogeneous vector Dirichlet BVPs

Next, we compute the eigenmodes in a brick cavity. The computational domain
is defined by a = {2+ + $2- with mm, 1) x (0,0.1) x (0.5, 1) and $2- = (0,1) x
(0, 0.1) x (0, 0.5). The relative permittivity in {2+ and 52‘ are chosen to be 57- = 1, 2.
This numerical experiment is identical to the one designed in [64] to validate the
classical vector FEM code. The convergence data (error in L2 norm) for eigenmodes
are presented in Table 3.3. As is evident from this data, the results obtained converge

rapidly with increasing order and refinement.

3.4 GFEM solver for piecewise homogeneous domain with curved inter-

face

3.4.1 Basis function

Until now, the vector GFEM solver has been developed for homogeneous domain

and piecewise homogeneous domain with flat material interface. When solving large

87

 

 

 

 

 

 

 

 

 

 

 

p=2 p=3

. . h=1/2 6.2e-4 1.88-5

magma h=1/4 2.094 9.1e-6
h=1/6 5.6e-5 6.1e—6

. h=1/2 0.019 1.7e-3
233:: h=1/4 4.9e-3 7.3e-4
h=1/6 1.833 3.7e—4

 

 

Table 3.3. Error in eigenvalues computed using nap

problem such as antenna design, air plane design and so on, the curved material in-
terface should be analyzed. When high order accuracy of solution is required, high
order accuracy in both basis functions and geometry description are necessary. In
classic FEM, there are two ways to obtain satisfactory geometry description. First,
we can descritize the regular mesh fine enough to describe the geometry, but it limit
the geometry description to first order accuracy. Second, we can utilize the high
order mesh based on the geometry transform which realize the high order accuracy
of the geometry description, but the property of vector basis function, especially the
divergence free property, is spoiled. Those drawbacks come from the continuity re-
quirments of basis function space on the edge of the mesh in classic FEM. In GFEM,
as the continuity requirement of the basis function can be realized only by the parti-
tion of unity function, we have more freedom in constructing the local approximation
function on the material interface. More specifically, the construction of the basis
function needn’t depend on the geometric transformation. This property makes it
possible to obtain high order accuracy of geometry transform and divergence-free
vector basis function together. In this section, we will propose a GFEM solver to
deal with the curved material interface and demonstrate the convergence of those
basis functions.

Figure 3.15 show a patch Q,- on the curved interface. The constitutive parameters

are different on two side. The tangential direction and normal direction to the

88

 

(22 (21

\..

Figure 3.15. Figure of patch on the curved structure

 

 

 

interface are denoted by t and 7‘1 on the interface. If the fields in both homogeneous
domain is denoted by E1 and E2, then on the interface, those fields should satisfy

the following criteria:

Etl = Et2

Enl # En2
BnEti 7f anEtZ (3-42)
BnEnl 7E BnEnz

V-E=0

The above criteria show that the tangential component of the field are continuous
and the normal component of the field are discontinuous. The normal directives of
both tangential and normal components are discontinuous as well. The fields are
divergence free in both homogeneous domain. As the vector basis function for homo-
geneous domain is continuous everywhere, the additional vector local approximation

flmction fadd should be defined in the patch to create the appropriate discontinu—

89

ities. The additional function fadd is not necessary to be defined on both sides of
the interface, we can just define fadd in one side such as (21 in Figure 3.15. Then

the criteria for the fadd are the following:

ft=0

fn 7'é 0
anft 75 0 (3.43)

anfn ¢ 0

V-f=0

That means, continuous component of the fields result in the zero component of
fadd on the interface, discontinuous component of the fields result in the nonzero
component of fadd- The addition function fadd is also divergence free. The last
criteria is very important to get the high order approximation of the divergence free
field.

The most straightforward way is to construct the addition functions based on
geometry transform. By doing so, the curvilinear structure 91 is transformed to
a regular structure. After defining the additional basis function in that structure.
we transform it back to the original structure. Unfortunately, we found that when
the additional function satisfies all the continuity requirements on the interface, it’s
not divergence free. In GFEM, as mentioned before, the continuity requirements of
the basis function on the edge of the patch can be realized only by the partition of
unity function. That means the construction of the addition function fadd can only
depend on material interface. Next, the construction of the addition function will
be introduced in detail.

We defined the local transverse and longitudinal coordinate as [0 and 2, then the
material interface can be defined as 20 = P(p) Vp 6 Pp. I‘p denotes the projection

of the interface in patch Q,- on the p plane, and we can define the fadd at any place in

90

 

 

Figure 3.16. Local coordinate of patch on the curved structure

the domain I‘p x (— inf, P(p)) as shown in Figure 3.16. The fadd can be formulated

using

v =V x (qsmé) +V x V x (457,2)
f(p, 2) =V(p, 2) — [V(p, P(p)) — fi(V(p, P(p)) - 73)] (344)

+ 2Vp- [V(p, P(p)) — fi(V(p, P(p)) - fi)](z - P(p))

where the it represents the unit vector normal to the interface. This vector can be
expressed mathematically by it = 2‘: — VpP(p). Scalar function aim and qbn can be

any function. For pth

order accuracy, function rim and dm should be of order p + 1
and p+2. The basis function defined above can satisfy all the criteria, including the
divergence-free criteria. Figure 3.17 shows 2D examples of the addition function. In
these figures, tangential component of the addition function is zero on the interface,

and the normal component is nonzero.

91

 

 

  
      

l|\\\e\\\\
| \ \\\\\\
.\\\\\\\
...-“‘

...,a’»»‘
::I//////
I III/////

 

 

 

 

 

l‘l“ll‘l‘

    

.........

Il'el‘“‘1

‘ “‘i.

 

 

 

(b) Inhomogeneous 02 mode

(a) Inhomogeneous 01 mode

 

llllllll

[riff/,1

 

 

 

 

““““

\\\\\\§

 

 

 

£52157,

 

III? III

 

 

 

 

 

 

 

(d) Inhomogeneous 11 mode

(c) Inhomogeneous 10 mode

 

 

 

 

 

 

 

  

I:

..\\‘\\~

..,....

 

 

 

(f) Inhomogeneous 20 mode

(e) Inhomogeneous 12 mode

 

 

:..:
- - ‘ \\\\\\
-\\\AAAA\
.\\\\\XR\
...\\\\x‘
..,,,,,/
.,//////
all/IIIZ

| III/[’1

 

 

 

 

 

 

(h) Inhomogeneous 22 mode

(g) Inhomogeneous 21 mode

Figure 3.17. Vector plots of inhomogeneous local approximation function V,- for

curved interface

92

The zero tangential component and nonzero normal component of the fadd is

proved mathematically as follows:

ft(p. P(p)) =Vt(Pa P(p)) - [Vt(p,P(p)) - 0]
+ [23 - (2 - 7073le - [V(p. P(p)) - fi(V(p, P(p)) - fi)l(P(p) - P(p))
=0
fn(p. P(p)) =Vn(p,P(p)) - [Vn(P:P(P)) - 1101(1), P(p)) - 11)]
+ (5 - fiW-Vp - [V(p, P(p)) - 7"Iii/(p, P(p)) - 71)] (P(p) - P(p))

=fi(V(p, P(p)) - a)

750
(3.45)
The following is the proof of divergence free property of fadd-
V - fn(p, P(p)) =V - V(p, 2)
- (V(p, P(p)) 41)]
(3.46)

+ Vp - [V(p, P(p)) - rib/(p. P(p)) 41)]

=0

The functions introduced herein demonstrate means to develop additional basis
function for curved surface. For completeness, we first introduce basis function
developed with geometry transformation. A curvilinear structure in cry-domain can

be mapped to a regular structure in én-plane. Transform functions are

3: =$(€,n) (3 47)

y =y(€,n)

93

 

..........
o

a
'o
'-
yo

 

 

 

 

 

 

Figure 3.18. Geometric transform

Then Jacobian matrix is defined as:

6:1: 6
[n= 335% as)
6x 2],:
35 n
and L] | is J acobian. If we use u and ii to represent the vector function in curvilinear

structure and regular structure, then we can construct the addition function in

curvilinear structure by the following transform formula:

11 =[J]"1fi
an)

VXHfl§UFVXfi

This formula can be easily extended to three dimension.

In both of the above two methods, the material interface is described by some
sample points on the interface. Locations of other point on the interface is computed
by interpolation of those sample points with Lagrange polynomials as interpolatory

functions.

94

3.4.2 Numerical Experiments

Next, several numerical results will be shown to demonstrate the accuracy of this vec-
tor GFEM solver. First, we shall demonstrate the convergence of the local approxi-
mation function. Next we shall demonstrate the convergence of the basis function.
Finally, one large two dimension problem will be used to demonstrate the feasibility
of this solver. At first, convergence of the local approximation will be demonstrated.
Figure 3.19 shows one patch on the curved interface, the curved interface is a seg—
ment of the circle with radius 0.48/\ and cross the center of the rectangle which can
be described as x2 + (y + 3)2 = 32 with I: = 2MB. The rectangular computation
domain is (—0.5,0.5) x (-0.5,0.5). The analytical results in the patch is the total
fields when a TE plane wave propagating to the dielectric coated PEC cylinder with
radius of the PEC 0.24/\. The patch is on the boundary of the dielectric coating
with the side of the patch parallel to the line between the center of the patch and the
center of the PEC. When we change the size of the patch, we keep the center of the
patch on a fixed point on the boundary of dielectric coating. The relative L2 error
is computed on the 100 x 100 sample points uniformly distributed on the patches.
The additional scalar function is chosen as tensor product of Legendre polynomials.
The Figure 3.20 show the error based on both of two methods introduced above.
The dashed line show the error of local function based on geometry transformation.
From the figure we can find that when the order of the function is high, the accuracy
of the approximation is spoiled by the geometric transform. This results shows the
importance of divergence-free criteria in creating the additional fimction. The solid
line show the additional function based on [3.44]. It shows the h— p- convergence of

the additional function. The error is excellent

95

 

 

 

 

R = p.42»

 

Figure 3.19. Figure of patch on the curved structure

 

    
   
 

   

 

 

 

 

 

 

10“ . . . . . . ,
i
-2
10 1‘
10.3 0" . I 1,
1
E :
2 4
"—1 10 ‘i
3
— P81 Cunt. lnte. .
“'5 —P-2 ourv. late. 1
—P83 Cunt. lnte.
-O-P-1 Geo. Trans. ‘
.6 - O - P-2 Geo. Trans.
10 -O-MGeo.Trans.-§
10"

Figure 3.20. h—, p— convergence of the local approximation

96

Next figure shows different errors when the center of the PEC circle is moving
around the patch, e.g., in Figure 3.21 the a is changing from 0 to 27r. In this case, the
size of the PEG and dielectric coating are the same with before. The center of the
patch is still on one point of the dielectric coating’s boundary. In this experiment,
the size of the patch is 0.05). and p = 3, The relative L2 error is computed on the
100 x 100 sample points uniformly distributed on the patches. The Figure 3.22 show
that the error is always smaller than 10-5.

Next, we will demonstrate the convergence of the vector GFEM solver in a simple
problem. As shown in Figure 3.23, there is a dielectric coated PEC with the radius of
the PEC and dielectric as 0.06) and 0.13/\. The computation domain is bounded by
two circles with radius 0.095A and 0.16A. Relative permittivity of dielectric coating
is 2. Dielectric boundary condition is imposed. Incident field is a TB plane wave and
propagates along —:f: direction. The L2 error is computed on the 100 x 100 sample
points uniformly distributed in the domain. The Figure 3.24 shows the convergence
of the basis function. The fields with h = 0.055). and p = 4 are also shown in
the Figure 3.25,Figure 3.26. Figures show the analytical result and error of Ep and
E45 respectively. As the Ep is normal to the interface, it is discontinuous on the
interface, when Egg is continuous on the interface. The figure demonstrates that
property and shows that the error is small and distributed uniformly throughout

the computational domain and is not restricted to the surface.

97

 

 

 

 

 

\

Figure 3.21. Figure of rotated curved structure

 

270

Figure 3.22. Error with different curved surface

98

 

 

 

Figure 3.23. Geometric description of dielectric-coated PEC cylinder

 

 

 

 

 

 

 

10"
l
10.3 1,
E 4
O
C
N-l ‘ ‘
1o 1
10'5 —P-1 1:
+922 3
—o—P-3 I
10" - 1,,
10
110.)

Figure 3.24. Convergence graphs of GFEM solver with Dirichlet BC

99

AWMMEP

.8 llvlllm
41111111

"unmet” ~
1’1““““f.......n9“
. fl

 

(b) Error for E,

Figure 3.25. Analytical and numerical result of p component of the field

100

 

(8) Analytical result for E¢

Enorot P

 

(b) Error for 13,,

Figure 3.26. Analytical and numerical result of ()3 component of the field

101

Last, a 2D large problem will be analyzed. As shown in Figure 3.27, there is also
a dielectric coated PEC with the radius of the PEC and dielectric coating as 2.7/\
and 2.85A. The relative permittivity of dielectric coating is 2. The computation
domain is truncated by a circle with radius 3.0A. The incident field is TE plane
wave propagating along —:E direction. Dirichlet boundary condition is imposed on
the boundary. The distance h between the neighboring nodes is 0.43/\, the order
p = 4, numerical results are computed on two circles inside and out of the interface
and close to it with radius 2.82/\ and 2.88/L The results are shown in the Figure
3.28, Figure 3.29. Results show that the approximation is excellent and the relative

error is 2 x 10‘3.

 

Figure 3.27. Figure of patch on the curved structure

102

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

x 10
_E
4 E
- - - Mm
3.5 r .
a
r.-
3 3' ‘
g 2.5 r '*
.2 2 ‘
'o
C
-Q 1.5) «
Ill
1 l' ..
0.5
o l I
0 1 2 3 4 5 0
9
(a) E, in the interface
0.018 V I V I L: I
_ Em
0.010 * E "
- - - W
0.014L 1
a
g o 012 - .
3
0 0 01
o , -
f
3 0.008 t ‘
a
I:
-, o.ooe- -
m
0.004 P ‘
0.002 " ‘
o l J
0 1 2 3 4 5 6

(b) E? in the interface

Figure 3.28. Analytical and numerical result of large 2D inhomogeneous problem in
the interface

103

x10

 

 

 

.

I

I

I -4
m "1
ii

 

5“
on

u

N

E out of dielectric coating
d N
in in

p

 

O
on .5
.1);

 

 

 

0°
.s
N
u.
A
0|
0

(8.) Ep out of the interface

 

0.018 V U I I I _

0.016 r E

 

 

0.014 '

0.012 '

0.01 -

0.008 r

0006*
I

E 9 out of dielectric coating

0.004 '

0.002 '

 

 

 

(b) E4, out of the interface

Figure 3.29. Analytical and numerical result of large 2D inhomogeneous problem
out of the interface

104

 

 

3.5 Overcoming linear dependence and condition number issues

Another problem that plagues GFEM, as it does any higher order method is the
condition number of the resulting system. Indeed, this is more of a problem here as
opposed to classical FEM as the space of basis functions is not interpolatory. The

vector basis functions in ith patch ¢i(r)V x [621513;] and ¢i(r)V X V x {59?}?

l a...
not orthogonal. In fact, they may even be dependent on each other. Consequently,
a good starting point to obtain a handle on the condition number is to redefine
basis functions in each patch so that they are orthogonal to each other. To this end,

we compute a symmetric matrix of inner product of basis functions {S}, in the ith

patch as defined by
{L}..- = [a dam-{S}? (3.50)

2

where {S},- is a vector whose elements are all basis functions that are defined in
the ith’ patch. Representing the singular value decomposition of the above matrix
as {L},;i = UzTDiUz- where U, is an orthogonal matrix, and D2- is a diagonal
matrix that contains the singular values. We can then redefine the basis functions
as {S}; = Ui{S}z-. It is apparent from this redefinition that the inner product
of any two basis in the domain (2,- are orthogonal. The basis functions related to
small singular values in D,- are dropped. These basis functions can be readily used
in all the quantities derived earlier and has been found to considerably reduce the
condition number of the resulting system. More specifically, consider the bilinear
form Ah(u, w) = .77}, (w) in (3.29). Using the span of basis functions defined earlier,
one can rewrite (3.29) as a matrix equation of the form Zj [Alij {x} j = {f }z where
a {x}j is a vector of unknowns coefficients of basis functions in Qj, sub-matrix
[Alij = Ah(u,w) for uap E Qj and w E (22-, and the corresponding right hand
side, {f l2 = Fh(w). It follows that given the redefinition of basis functions, one can

alternatively redefine the vector of unknown coefficients as {17:}z = [U],{x},. Using

these definitions, matrix equation arising from the inner product can be written as

[Aijl =[UlilAlz'lelg

{fli =lUlz'{f}i

(3.51)

It is apparent from the preceding description that the SVD-based preconditioner is
developed on a patch by patch basis. This implies that the number of operations
scales as proportional to the number of patches in the domain. For instance, in a
2D problem, if there are N patches and pth order basis functions are being used,
then the total cost of performing the SVD scales as 0(Np6) and imposing required
transformations scales as 0(Np4). This cost is derived using the fact that the
matrices used in the transformation are block diagonal and not full. Given that
p = 0(1), the overall operation is typically proportional to the number of patches
used in the domain. Thus, it is not very expensive to retrofit existing codes with
this preconditioner and, in this case, considerably reduces the cost.

The error bound of the new set of basis functions can be obtained. Assuming
that the threshold of singular value in [D]m, used to remove the dependent basis
functions, is defined as ed. The approximated field can be represented by two set
of basis functions as uap = 2k akS;c = Zj chj. In each patch, we have {a}?w =
{C};T[U]z where {a},- and {c},- are vectors representing the coefficients of original

basis functions and reformatted basis functions in ith

patch respectively. In all
N reformatted basis functions, we assume that N1 of them are accounted with
the singular values larger than ed, then rest N2 =-- N - N1 of them are those
which singular values are smaller than ed which are removed as dependent basis
functions. Number dz- denotes the singular value related to the basis function Si-

For convenience, we assume that the first N1 basis functions has singular value

larger than ed. Constant cmam represents the largest coefficient of reformatted basis

106

functions. The error bound can be derived as follows:

 

 

 

 

N1
“‘1 — fiap'lLZm) = u - Zcisz'
2:] L262)
N ~ N ~
= u - (Z 0282' - Z 6282')
N ~ N ~
= u - ciSz' + Z Cisz'
i=1 [9(9) i=N1+1 Lgm) (352)

N
g “u — nap” 'l" 2 102' “Silleml

z=N1+
N
s llu-uapllL2m, + :3 we?
i=N1+1

5 H“ _ uaPHL2(Q) + NZOmanx/Ecl

Next, several results will be shown to demonstrate the efficacy of this SVD-based
preconditioner. Figure 3.30 shows the condition number of the final system with
different order of basis functions and different threshold ed. Figure 3.30 shows that
when the order of basis function increases, the condition number of the original
matrix (ed = 0) increases very fast, but with different threshold, the condition
number increases slowly. Figure 3.31 shows the error of the solver with different
threshold of preconditioner. It demonstrates that the errors are almost unchanged
while the condition number of the system is reduced significantly. Figure 3.32 shows
the condition number of the system with order of basis fimctions changing from 1

to 10. With different threshold, the condition number slowly increases.

107

 

 

 

 

 

 

 

 

10 Y I I I' 1
—— liveshoid 10 "0
10 — 0mm .10 ‘9

10 - — W :10 '“G
x —— threshold :0
E I/
g 10a» // .1
3 /
‘5 ,/

fl /

E 10 l/’ '
3 /’
.5 ”
=96 108 [’1/ "
8 /
o

106 D- 4

104 l L l l I

1 1.5 2 2.5 3 3.5 4
OrderoiGFEMsoiverm)

Figure 3.30. Condition number of system with differen order of basis functions(p =
1,...,4)

 

 

 

 

 

L2 Error

 

 

 

 

1 1.5 2 2.5 3 3.5 4
Orhr of GFEM solver (p)

Figure 3.31. Error of the solver with differen order of basis functions

108

 

102°

 

 

 

 

 

 

—-Thmho|d: 101° [‘74:
——Thmhold: 10’12 ,,I// /
‘5 ——Thneshold: 10'“ /// /
5 1° —-——Threehold:0 / / / “
_ - - — BeelfitlineforThreehoid:O //
.2. ./
‘5 /
E 1010 / ’ / -
g /
c /
E 10‘5 ~ / / -
/ /
10° - -
10° 101
OrderofGFEM solver(p)

Figure 3.32. Condition number of system with different order of basis functions(p

= 1,...,10)

109

CHAPTER 4

TIME DOMAIN GENERALIZED FINITE ELEMENT METHOD

4.1 Introduction

Tfansient analysis techniques that are based on differential equation based tech-
niques have been popular, but this is largely due to finite difference time domain
technique [65]. However, recent work in time domain finite element methods (TD—
FEM) [1] has revived some of the interest in further development of this direction of
research [66—68]. Recent developments in this area are truly exciting as they offer a
viable and efficient alternative to FDTD. They inherently offer some advantages: (i)
they conform to the object and (ii) it is possible to use higher order basis fimctions.
An apparent disadvantage seems to be the fact that a straightforward approach im-
plementation of TDFEM seems to be the cost required to obtain the solution at
each time step. A novel approach to solving this problem was proposed in [1], and
other recently developed methods permit more efficient time stepping.

The generalized finite element methods (GFEM) provide an intriguing analy-
sis methodology, in that, it is a general framework for differential equation based
analysis and permits the use of non-polynomial space of basis functions. Further-
more, the basis function space does not rely on an underlying mesh. This said, the
enthusiasm has to be tempered largely due to the fact that while the potential of
this method is considerable, progress still needs to be made to realize the kind of
efficiency and flexibility that current day finite elements possess. The aim of this
chapter is two-fold; (i) methods for time domain analysis using generalized finite
element; and (ii) analysis that is easily portable to standard finite elements as well.
Some of the inspiration for this work is drawn from our earlier work on time domain
integral equations [69].

Classical TDFEM typically employs either central difference or the Newmark

110

method for time stepping; the latter is unconditionally stable, and both are accurate
to second order. However, an alternative to these methods, is exploiting the fact
that in any of these solver, we seek a response to a bandlimited source. This implies
that one can exploit this to create both an alternative matching scheme, and well as
a more efficient iterative scheme. Thus, the principal contributions of this chapter
are as follows: (i) we will develop a time domain generalized finite element solver
and prove/ demonstrate convergence; (ii) we introduce a marching scheme based on
approximate prolate spheroidal wave (APSW) functions and analyze stability; (iii)
other marching schemes are analyzed for stability within the context of the TGFEM
is provided; (iv) higher order convergence (time and space) is demonstrated. This
chapter is organized as follows: Next, we will propose two high order time stepping
schemes with error analysis. Then the stability analysis of the high order time
stepping schemes will be provided. After that, we will propose one iterative solver.
Finally, we will demonstrate the accuracy and stability of the technique with a host

of numerical examples.

4.2 Formulation for Time domain GFEM solver

4.2.1 Discretization of time domain wave function

Consider a domain 9,- whose boundary is denoted by 69 z: I‘ = U2 Pi- The electric

(or the magnetic) field by u(r, t) satisfies the equation

 

2 r
(V X [fiVx] + “0% + 30);) u(r,t) = _8f(at,t)

B,- {u(r,t)} = gi(r, t) for r E F,-

(4.1)

where f (r, t) denotes the electric (or magnetic) source that is bandlimited to an
angular frequency wmax and is approximately quiescent for t < 0. In the above
equations it is assumed that r 6 Rd for d = 2, 3, B, are differential operators, and

gi(r,t) is the function that is imposed on I”). In what follows, we will describe the

111

manner in which solution to these equations may be obtained via the GF EM.
As is done in all discrete solution to PDEs, the function is represented in terms

of space time basis functions of form

Nt N3
u(nt) = Z Z ujnm —jAt)sn(r) (4.2)

j=1n=1
where Nt and N3 represent the number of temporal and spatial unknowns, T(t) and
Sn (r) are the temporal and spatial basis functions, At is the time step size, and ujn
are the unknown coefficients need to be determined. The spatial basis fimctions are
formulated in terms of GFEM space.

Next, we address the temporal basis functions or equivalently time stepping
approaches. Three different types of time stepping approaches are investigated: (i)
backward Lagrange; (ii) Newmark, and (iii) those resulting from using approximate
prolate spheroidal wave (APSW) functions. As an aside, all time stepping schemes
can be mapped into equivalent temporal basis functions. In what follows, we will use
time stepping and temporal basis functions interchangeably. The APSW functions
used here are those introduced by [70] and have been used extensively in the solution

to time domain integral equations [69]:

 

T0) = Sin(mea$t) Sin [(1)/(“71%)2 — I] (4.3)

meaxt sinh(a)\/(N;1—t)2 — 1

 

 

In time domain, APSW function can be used as interpolatory function to the uniform
time step size. x is oversampling factor such that the sampling frequency ws =
xwmag; and At = a/ws. APSW function also has compact support. p is the half-
width of the prolate, and T(t) = 0 for It] > (Np +1/2)At. a = «Np(x — 1)/X is the

time bandwidth product. The error of approximation based on the APSW function

112

can be represented as

N
000— Z f(kAt)T(t—kAt)Is
kz—N

 

32,73,101) ~ 0(e-NpflX-D/X) (4.4)
This fimction is band-limited to wmax Figure 4.1 shows the APSW function in time
and frequency domain. While these functions are bandlimited and approximately
time-limited, they do not satisfy the march condition, i.e., they do not results in
causal system that would enable a marching-on—in-time process. Means to overcome
this have been prescribed [69], and will be summarized in the next subsection.
Lagrange polynomials can also be used as high order temporal function. Follow-
ing definition shows that it can also be used as the interpolatory function in time

domain.
(t — At)(t — 2At) - - - (t — pAt)
MN)"

 

T(t) = (4.5)

But when order p is increased, the function is getting more and more oscillatory,
e.g. the function is not band-limited. This makes the final system unstable. This

property will be shown in the result section.

113

APSW function in time domain
12 I . I

 

0.8 *- .
0.6 ~ .
0.4 ~ -
0.2 r -

VVV Vv -
-02 . 1 .

'—2 —1 0 1 2

‘ 1110-10

 

 

 

 

 

(a) APSW function in time domain

APSW function in frequency domain

1.5 I I T M l I I I I

 

 

 

 

 

 

 

0.5 ~ ~
0
_05 1 1 1 1 1 1 1 1 1
-5 -4 -3 -2 -1 O 1 2 3 4 5
0) 9

x10

(b) APSW’ function in frequency domain

Figure 4.1. APSW function in time and frequency domain

114

4.2.2 Causal system for time domain wave equation

Finally, to derive the MOT system, we will use the approximation space defined
earlier in (4.1). It will be assumed that the boundary condition on different portions
of are defined to be either of the Dirichlet or Newmann type, i.e., [V x u(r, t)] X ft =
g1(r, t) and h - u(r,t) = g2(r, t) are imposed on I‘n, u(r,t) x 11(r) = g3(r,t) is
imposed on I‘d. Using Gelerkin testing in space results in the following bilinear

form:

Ah(u, (Sm(r), t)) = fh
.Ah(u, (Sm(r),t)) = [fidnfidlr—NV x Sm(r)) - (V x u(r,t))
62u
+ [a d01<r)sm(r) - 571m)
Bu
+ [a d9fi(r)sm(r) . amt)
+ Cl f9 dQ(V - Sm(r))(V . u(r,t))
1 ..
'1” [I‘d drm(v X Sm(r)) ' (Ll(l‘,t) X Ti)
+f dF-l—(S (r) x n) - (V X u(r t))
pd u(r) m ’ (4.6)
+ C2/P dF(Sm(r) - h)(u(r, t) -h)
+ C3 fr‘ddI‘(Sm(r) x 7‘7.) - (u(r, t) x h))
f}, = — Ldflsmh) - %(r, t) + A” dI‘Sm(r) - g1 (r, t)

+ c2 fr dI‘gzc. t)(sm(r) ~71)

1 ..
+ I‘d drag-$3039 ' [V X Sm“) X n]

+ C3 Ad ng3(r,t) ~ (Sm(r) x h)

where cl, c2 and C3 are constants to include penalty terms and additional bound-

ary condition respectively (we have shown in [71] that while these terms are not

115

absolutely necessary for h—, p— convergence, the convergence is smoother if these
conditions are imposed).

Substituting (4.2) to (4.6) and using point testing in time, results in an MOT
equations. For different time stepping schemes (equivalently temporal different basis

functions), these equations may be written as 232- 1 Zin— j = Vi- In this equa-

K
tion, [13]” = ujn, le'lm = fmltzti: and Z,- = Altzti, 111 + 102 + 1 is the number
of time steps included in the temporal support of the time basis function. When
using backward Lagrange or Newmark time stepping schemes, the resulting system
of equations is causal, e.g. 101 = 0. This is not the case when using APSW flmc-
tions. The solution to this problem maybe obtained using an extrapolation scheme
prescribed by either [69] or [72] to make the MOT system causal. More specifically,

extrapolation coefficients kg, is applied to approximate the unknown vector 123.. j

With j < 0 based on the Nsamp solved VBCtOT IZ-j With j Z 0 and Nsamp S 52 +1.

0
1+1) = Z hglj+q (4-7)

(1:1 "Nsamp

I

Explicit expression for the matrices can be easily obtained. As is evident from the
above exposition, these equations can be used in a stepwise manner to obtain the
fields for all time. Next, the error bound of the time domain GFEM solver will be

analyzed.
4.2.3 Error bound of Time domain GFEM solver

The error bound of this time domain solver can be analyzed. We assume f (:13) and
9(1) are approximation of linear operators f (.13) and g(r). Given the error of the

approximation

f (I) =f($) + 6(f (17))

90‘) =90?) + 64905))

(4.8)

116

Then the error of approximation of operator f (g(r)) is

f(g(rv)) =f(g(x) + eats)»
41002)) + f(e<g(x>))
=f(0(x>) + was» + f(e(g(rc))))
=fa<x>> + O(e(f(§(:v)))) + 002190)» (4.9)
=f(§z(:c)) + 0020002)) + cam») + 00190:»)
=fa(a:)) + O(e(f(g(x)))) + O(e(gc») + O(e(f(e(g(x)))))

#1100?» + 0(6(f(9(-’II)))) + O(e(9($)))

The difference between both sides of z is small when e(g(:r)) << g(r) and
e( f (33)) << f (:r). The theorem shows that the error of two operators can be ap-
proximated by the summation of the error of two operators. So in (4.6), errors of

the .Ah and f), are

6(Ah) =0(e(V X u)) + O(e(u)) + O(e(V - u)) + O(e(u X f1))+ O(e(u - u))

2
meagre)» +o<e<§,—f(t») +0(e(f(t))) (4-10>

e(f'h) =0(e(u)) + O(e(u ~ 71.)) + O(e(u x 71)) + O(e(V X 11 x 71))

In [73], the order of the error of GFEM basis functions O(e(u)) is the same as
error of local approximation function. In the above formulas, O(e(u», O(e(u - 13.)),
O(e(u X 73.)) have the same order accuracy and O(e(V x u x 71)), O(e(V x u)) and
O(e(V - u)) are one order lower.

In time domain, magi—22m»), O(e(%f(t))), 0(e(f(t))) are used to represent
the error of the numerical approximation of the time domain operators. The follows

are the errors of different time stepping schemes:

117

 

 

Backward Difference:

mag-f0») =0(At)

O(e(g—t-m)» =0<At>

0(6(f(t))) =0
Forward Difference:

O(e(gm») =0<At)

meg—m» =01At>

O(e(f(t))) =0
Central Difference:

wagers») =0((At)2)

O(e(gfc») =0((At)1)

O(e(f(t))) =0
p—order Lagrange polynomials interpolation:

$22—10)» =0<<At)P-1)

angst») =0((At)p)

O(e(

O(e(m») =0<<1At>P+1n

(4.11)

(4.12)

(4.13)

(4.14)

Approximate prolate spheroidal wave(APSW) functions with width parameter

118

AT

O(e(g— 22:10)» O(A—tze—“ZN/wom)

O(e(g, ft())) =0(Aite—"2N/woet) (4.15)

O(e(m» =0(e-“2N/w0“>

4.3 Analysis of stability

The analysis of stability is a very important issue in this time—marching scheme.
The criteria of stability in discrete signal system is utilized in this paper to analyze
whether and when the system is stable. The bilinear form equation (4.6) can be

expressed as:

[Sl{u(ts)} + [Q]{ii(ts)} = f(ts) (4.16)

where (z' — 1)At < t3 _<_ z'At and the current time can be represented as t; = iAt+t’.
Approximating the time domain result by temporal basis function T(t) where T(t) =

0 for |t| > (N + %)At, the equation can be discretized further as:

i+N i+N
Z [S]{I}jT(-J'At+ts)+ Z [Q]{1}jT(-jAt+ts)={V(ts)} (4-17)
jzi—N j=i—N

With extrapolation algorithm in (4.7), this equation can be rewritten as:

Z ([Sl{1}j0j—i+N + [Q]{I}j0j—1+N) = {V(ts)} (4-18)
j=i—N
where
T((N — am + t’)
02' = (2' = 0,. .. ,N — Nsamp) (4.19)

T((N — 1))At + t’) + 211V“ T((N — j)At + t’)h§:{,vv
(i- =N+1—Nsampa---:N)

119

T((N — 2')At + t’)
52' = (z' = 0,. . . ,N — :Nsamp) . (4.20)
T((N — 1')At + t’) + 2].]; _N+1 T((N — j)At + t’)h;.:11‘[

(2=N+1—N3amp,...,N)

The period time At should be extracted from 6.,- and 5;, so by stretching the temporal
.. .. I

function T (t) = T(th) and using t’ = i—t, The following equation with normalized

form is obtained:

21181a.+[§-—]—’)’,>{I}.+j = {van} (4.21)

where
T((N—1H?)
i=0,...,N—N
(1,: ~ ( 2N famp) ~ ,N (4.22)
T((N—0+t’)+Zj=N+1T((N—j)+t’)h;_N
(i=N+1—Nsamp,...,N)

f((N — i) + i’)
b,= .. (i=0,,... N‘stamp) ~ , (4.23)
T((N ) +t’)22£N11T((N -j> + ”021%
(i=N+1—N3amp,...,N)

The z-transform of the left side of the above equation result in the equation:

 

{gases [21:9 )2le i): (4.24)
_Z__2N;_— —0 biz i —-1

———.{I},- =<At> [Q] [S]{I}j (4.25)
Ego“ ’i‘z

120

 

 

 

 

 

 

i=0 i=1 i=2 i=3 i=4 i=5 i=6
4 bi -6.0le—5 1.048—3 -8.93e—3 -0.973 3.964 -4.976 1.994
Nsamp = az- 0.0 0.0 0.0 0.0 0.0 0.0 1.0
5 bi 1.0-6.016—5 1.04e-3 0.912 -4.658 9.489 -8.660 2.915

 

 

 

 

 

 

 

 

 

 

Table 4.1. The coefficients of characteristic equation

It’s easy to find that ‘42in bizi)/ (211:0 aZ-zi) is the eigenvalue of matrix
(At)2[Q]_1[S] and written as A, then

N

2(1),- + Aa,)zz = 0 (4.26)

i=0
The condition of stability for the above system, or called the A-system, is that
all the roots of equation (4.26) are in the unit circle. If there is a number
Amaa; when VA 6 [0,Amam], all roots of z are in the unit circle, e.g. A—system

is stable, then the time stepping scheme is conditionally stable with the criteria

At < \/()1mag;) / (p([Q]‘1[S])). Analysis of the stability by characteristic equa-

 

tion (4.26) w.r.t 2: can be transferred to the characteristic equation w.r.t. m using
2 = [1%, where the stability criteria is that all roots of 11) should lies at the left side
of the complex plane. Then Routh criteria can be utilized to analyze the stability
of the system.

An experiment will be presented to demonstrate the feasibility of the above
scheme. In the experiment, APSW function is chosen as temporal basis function
with testing point t’ = 0, the parameter chosen is At = 2.33 x 10‘113, Nsamp = 4,
N = 6. The coefficients of characteristic equation is shown in the Table 4.1. When
/\ 6 [0,100], Figure 4.2 shows the stability of A-system based on both Rough criteria
and the numerical solution of A-system. When Rough criteria is implemented, value

1 shows that the A-system is unstable and value 0 represents that the system is

121

stable. It shows that VA 6 [0, 100], the A-system is stable. The appropriate At can
be chosen to make the approach stable. In the same figure, with different A, the
largest amplitude of roots of z in Figure 4.2 is computed by intrinsic subroutine in
Matlab. From the figure, we can see that they match with each other. Another
experiment is also presented with N sump = 5 and all other parameter unchanged.
The coefficients of characteristic equation is also shown in the table. From the Figure
4.3, We can obtain the conclusion that there is no Amag; can be chosen so that
VA 6 [0, Amax], the system is stable. So the time-marching system is unconditioned
unstable. Figure 4.4 and Figure 4.5 show the same results with smaller region of A

which make the details of the figures around the origin clearer.

4.4 Iterative Solver

The main bottleneck in using TD(G)FEM as opposed to F DTD is the fact that the
resulting systems are not explicit. This problem is exacerbated if the basis functions
are higher order. One way to overcome this problem in TDFEM is to develop basis
function that are orthogonal in each subdomain [1]. While this approach is efficient,
it does sacrifice higher order spatial convergence. Using a SVD approach to define
orthogonal basis functions in TDGFEM permits higher order convergence. Another
fact that can be exploited is the fact that the signals analyzed are bandlimited. This
implies that one can make use of extrapolation schemes described earlier, but within

the iterative scheme. We shall dwell on the latter in a little more detail.

122

 

 

 

 

 

 

 

 

 

1.5 T I r I I I T I I
-————4aqnntnxusofz
-—-—-snxmnybyrouflhuflena
1 i .
(15> .
O
_0.5 l 1 l l 1 l l l
O 60 70 80 90

10 20 30 40 50
A

 

 

 

 

 

 

 

 

 

 

1.5 I l I I l I I I I
————— nxnsofz
‘——--sunmflybyrouduuflene
119 "
(15- -
0
_0.5 I 1 l 1 l 1 l 1 1
0 1O 20 30 40 50 60 70 80 90

Figure 4.3. Stability of system with Nsamp = 5

123

 

1.5 l I I T I W I

 

— largest roots of:
-—-— stability by routh criteria

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

0.5 e -
0
_o.5 l L l l 1 1 l 1 M
1 2 3 4 5 6 7 8 9 10
A
Figure 4.4. Stability of system with Nsamp = 4
1.5 r I l I I j I I I
largest roots of z
—-— stability by routh criteria
1- ‘
0-5 I. l!
0 - fi
_0-5 I 1 1 l l l
1 2 3 5 6 7 8 9 10
A

Figure 4.5. Stability of system with Nsamp = 5

124

 

At any instant of time ti, the MOT system of equations can be written as
[Z]0{I},- = {V}. It follows from the above arguments that {I},- can be estimated in
terms of {I}, forj = i— Neg; —1,...,i— 1. Next, the matrix [Z]0 can be subdivided
into submatrices [Z]0,mn where the subscript {mn} indicates an interaction between

the mt

h and nth patch. In terms of sub-matrices, the iteration can be rewritten
as [Z]mm{I}i = {V},- — an[Z]mn{I}i where {I}; is an extrapolated estimate of
{I},-. This provides the first step of the iterative solver. One can build upon this to

provide the next guess, and so on. This iterative process is rapidly convergent for

an L2 error norm of 10‘s. It is not as robust for lower error requirements.

4.5 Numerical Experiments

In what follows, we shall present a series of numerical experiments that will serve to
demonstrate the accuracy and stability of the time domain GFEM presented in this
paper. These fall into three classes: First, we will focus on the stability issue when
different boundary conditions are imposed in homogeneous domain. Second, we will
demonstrate h—, p— convergence of spatial basis functions. Then, the convergence of
TDGFEM for different temporal schemes in two— or three-dimensions, homogeneous
or inhomogeneous domain is demonstrated. Finally, we will shows the results of
TDGF EM with iterative solver.

First, consider a homogeneous rectangular computational domain 0 = (0,1)2
with boundary F = U? Pi where {F1 : r E 0 X (0,1)}, {F2 : r 6 (0,1) x 1},
{F3 : r E 1 X (0,1)}, {F4 : r 6 (0,1) x 0}. In this domain, it is assumed that the

electric field is of the form

:c—rro

 

E(r, t) = g2(t + )/728:I:p{—(t + iicflfl/TZ} (4.27)

where 2:0 = 5.5 and 7' = 3.15x10‘93. The solution to (4.1) is solved by imposing the

analytical solution as either Neumann, Dirichlet or Impedance boundary condition.

125

In this experiment, all nodes are distributed uniformly, and rectangular patches are
used. The size of each patch is chosen to be 1.5 times the distance between the nodes.
The weight functions Wi(r) is a tensor product rooftops, and the local approxima—
tion function 112-M(r) and vi,n(r) are tensor product of Legendre polynomials. The
parameters of spatial basis functions are h = 0.38A and p = 3. APSW functions are
used for as temporal basis. The half width of prolate function Np and the number of
sample points for extrapolation Nsamp are 6 and 4, respectively. The time step size
is At = 2.33 x 10—115. Linear equation is solved by SVD techniques. Figure 4.6,
Figure 4.7, Figure 4.8 shows analytical and numerical results at location (0.44, 0.44)
with Dirichlet boundary condition, Neumann boundary condition and impedance
boundary condition. From these figures, we can see when Dirichlet boundary condi-
tion and Neumann boundary condition are imposed, the solutions are corrupted by
interior resonance modes. As is expected, unique solutions are obtained when using
the impedance boundary condition.

Next, we will demonstrate the h— and p— convergence of spatial basis functions
in time domain GFEM solver. In this experiment, the rectangular computation
domain and other parameters about meshes (size of the meshes, distribution of the
nodes) are exactly the same with the last example, so as the electric field solution
of the time domain wave equation. Impedance boundary condition is imposed.
When temporal basis functions are the same APSW functions with those in the
last example, spatial basis functions are changing with h = 0.76A, 0.38A,0.25A and
p = 1,2. Temporal basis functions and time step size are identical to those used
earlier. The total number of time steps is 6000. The results are L2 norm of relative
error evaluated at uniformly distributed 10 x 10 sample points at all 6000 time steps.
Linear equation is solved by SVD techniques, too. Figure 4.9 demonstrates the h-,

p- convergence of the spatial basis functions.

126

D Dirichlet BC is imposed.

 

10 I I r r f

 

 

 

 

     

 

 

 

——Analyticalresuit
—— -Numerieelresult
«AW-\N f‘ \ \ fix I“ A
105 - 1.1"”,3 1 [‘1',’ l‘,’1’1H]T,“,“( [(1, ,I(\|’\‘rl;‘|£
l l
.1] , l|1.|[|lll ""'1':'.lllll
m>100 ‘
10'5 ~
10.” 1 1 L M 1
4 6 8 10 12
«5) x10-8

Figure 4.6. Analytical and numerical result with Dirichlet BC

127

101°

Figure 4.7. Analytical and numerical result with Neumarm BC

Neumann BC is imposed.

 

I I l I

 

 

 

 

 

 

 

 

128

.,
”\f:‘l\f{;‘\\” Ami: \,\"/. Mfi\r,n{,r (WIN; 'N,‘r"\,,,\f fr“!

l!!!" i ‘I'I‘IIII,||Hf|ll'lll‘

1 1 l
l
-—Anaiytimlresun
- - -Numen'calrasult
6 8 10 12

«5) x10-8

lnundbnaaBCHshnpoaui

 

 

   
 

 

 

 

 

 

 

1010 I I I I I
--—-Amaflficah15uk
- — -hhunedamiasun
105 -
>~ .1
LU 100 I ”1‘11"“"11'1
in, 11111,I1,I,l1'1i',‘m
I 11lllI‘lIIlIl‘VU“MfrW
1I1. "'l'“"1'3111W11'1.111
-5 1 '11 il|l11.,111,{111~, 11
1o 1 1 l1 I111'l1'l 111,111. -
- .11-1l111 1111111111
,1 ll ,1 11(11‘1'11
ll i ll 1 “1‘,
I 1 lll‘l‘
10-1C l l l 1 l
4 6 8 10 12
“5) x10-8

Figure 4.8. Analytical and numerical result with Impedance BC

129

_ A 1 = 2.33e-11
1o . .

 

 

 

 

 

 

 

 

: —‘_P31 1
1 —9—P=2 .
’ i
1
10-2:
2 1
1?
Lu
‘5
E
2
“.1 _3
10 r '1
’ 1
10'4 - 00
1

Distance between nodesa)

Figure 4.9. h— p convergence of the spatial basis function

In the next series of examples, we will demonstrate the convergence of differ-
ent temporal schemes in two- or threedimension, homogeneous or inhomogeneous
domain. Three kinds of temporal schemes are tested: scheme based on APSW func-
tions, scheme based on lagrange polynomials and N ewmark scheme. First, in 2D
homogeneous domain, the electric field solution of the time domain wave equation,
the computation domain and parameters for the mesh constructor are the same as
in the last example, the spatial basis fimctions are also defined in the same way with
h = 0.38A and p = 3. Impedance boundary condition is imposed. When APSW
functions are used as temporal basis functions, the half width is varied from 3 to
6 and the number of sample points for extrapolation Nsamp = 4. When Lagrange
polynomials are used, the order changed from 3 to 4, and Newmark -16 method is
used, fl is chosen as 0.4, 0.6, 0.8. The time step size is At = 2.33 x 10’113. The re-

sults are L2 norm of relative error evaluated at uniformly distributed 10 X 10 sample

130

points at all 6000 time steps. Figure 4.10 shows the errors of all three schemes. Horn
the figure, we can see the error is reduced with the increasing of the Lagrange poly-
nomials’ order and the half width of the APSW functions. The errors of Newmark
13 method are almost the same with different 3. When we increase the number of
sample points for extrapolation Nsamp beyond 4, the system is unstable. This un-
like what is observed in time domain integral equation. This is the principal reason
for the lack of convergence of this scheme.

Next, we will demonstrate the similar thing in 2D inhomogeneous media. Con-
sider a 2D inhomogeneous domain, the computational domain is defined by Q =
12+ + $2- with 52+: (0,1) x (0.5, 1) and W: (0,1) x (0,05). The relative permi-
tivity in 9+ and Q- are chosen to be Er = 1,3, respectively. The fields are as follows:
Vr 6 52"“ the electric field is E(r, t) = :i:{2(t + y—jgflflrzexfl-(t + 3:6-y-QV/7'2] +
2R(t — (31—605) + y1yO)/T2exp[— (t— (y—C.0 5)+ y 1y0)2/72]}, and ‘v’r E Q" the elec-
tric field is E(r, t) = 1&2T(t2 + y;cQ._‘/§_,_ —C—Q)/T2exp[—(t + wfi+ KEV/T2]
where R: figm— 1+1F ,yo— _ 5. 5 and r— _ 3.15 x 10 93. The solution to (4.1)

is solved by imposing the analytical solution as impedance boundary condition. The

 

 

mesh parameters (size of the meshes, distribution of the nodes), the basis functions
in spatial domain and temporal basis in time domain are identical to those in the
last example. The results are L2 norm of relative error evaluated at uniformly dis-
tributed 10tz'mele sample points at all 6000 time steps. Figure 4.11 shows that the

result is similar to those obtained in two dimension homogeneous simulation.

131

 

1f : I T WI I I I I I I
: prolate function
1 —9— lagrange polynomials 1
—*— Newmark Method

 

 

LA-AA

 

 

 

10.1:

. ; LA..1

l
N V

v vw—vvvr
A A ; .JAL

Relative error in L2 norm
—L
O

I
U

10

'I

-LLI

 

 

 

 

10'4 1
2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7

N for prolate and lagrange

 

Figure 4.10. Convergence of different temporal scheme in 2D homogeneous problem

132

 

 

18 . I I I I I I I
.
’

 

 

 

 

 

 

 

 

prolate function 3
I -—*— Newmark Method ‘
10'1 :- -
g I
O >
c 1
NJ 1
5 2
E 10 1' ‘1
0 : 1
0 1 1
s 1 -
.9 b
o
m 1
10-31 ‘1
10.4 l l l l l l L I 1
2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7

N for prolate and lagrange

Figure 4.11. Convergence of different temporal scheme in 2D inhomogeneous prob-
lem

133

Next, we will demonstrate the similar thing in 3D homogeneous domain. Con-
sider a homogeneous computation domain 0+ = (0,1) x (0,1) x (0,1). In this
domain, it is assumed that the electric field is of the form E(r,t) = Q2(t +
#Q)/T2€1L‘p{—(t + EZ—xQP/TZ} where :60 = 5.5 and 7' = 3.15 x 10-93. All other
parameters, viz., distribution of nodes, overlap, definition of spacial basis functions
and temporal basis functions, time step size, etc, are identical to those used earlier
and just extend to 3D case. and the L2 error norm is evaluated at 5 x 5 x 5 uni-
formly distributed sample points. The results shown in Figure 4.12 demonstrate the
convergence of TDGFEM solver in 3D.

Next, we will demonstrate the feasibility of the solver in 3D inhomogeneous
domain. Consider a inhomogeneous computation domain (2 = 9+ + 9" with 52+ =
(0,1) x (0,1) x (0.5, 1) and Q" = (O, 1) x (0,1) x (0,05). The relative permittivity
in (2+ and {2" are chosen to be 61- = 1,3, respectively. In this domain, it is assumed
that the electric field is of the form E(r, t) = 3}{2(t+Z—_Ezfl)/T2exp[— (t+z—_c—zQ)2/72]+
2R(t— (Z105) + 2120)/72e$p[—(t— @795) +2220)2/T2]}, and Vr E (2" the electric
field is E(r, t) = @2T(t+ #x/g+5%0)/T2e$p[—(t+§:co—‘§\/§+y—Tz0)2/T2] where

1-— ,/§ _ _ —9 -
R: 1——3+\/—,T— -1—3+2\/—, 20:5.5andT—3.15>< 10 3. All other parameters, v12.,

 

 

distribution of nodes, overlap, definition of spatial basis functions and temporal basis
functions, time step size, etc, are identical to those used earlier and as before the
L2 error norm is evaluated at 5 x 5 x 5 uniformly distributed sample points. The
results shown in Figure 4.13 also demonstrate the convergence of TDGFEM solver

in 3D piecewise homogeneous domain.

134

 

1d). 7 I I I I I 7 V Y

—— prolate function
. —e— lagrange polynomials I
—l—- Newmark Method

 

 

 

 

-1

A A

Relative error in L2 norm
3

 

 

 

1 l l l l J l l

2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7
N for prolate and lagrange

 

10

Figure 4.12. Convergence of different temporal scheme in 3D homogeneous problem

 

10

 

 

 

 

 

 

 

 

 

E —~— prolate function 1
. —e— lagrange polynomials l
' -N-— Newmark Method *
1°"? 1
g I i
o . l
C . d
N
_l , i
.E 2
g 10 ? ‘3
¢ I I
o r i
> r d
'5
_«2 > ,
o J
a:
10-33 .. ':
1 0.4 l 1 1 L l l l 1 1

2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7
N for prolate and lagrange

Figure 4.13. Convergence of different temporal scheme in 3D inhomogeneous prob-
lem

136

Next, we will demonstrate the feasibility of the time domain solver for scattering
from PEC sphere. The radius of the PEC is 0.45 and the open domain is truncated by
a concentric sphere with radius 0.65. In this domain, it is assumed that the incident
electric field is of the form E(r, t) = g2(t + Egfl)/T2erp{—(t + fl—CfiQ)2/T2} where
100 = 7.5 and T = 2.57 x 10—93. The solution is solved by imposing the analytical
total electric field as Dirichlet boundary condition. The parameters of spatial basis
are h = 0.57A and p = 3. The time step size is At = 9.52 x 10—113. All other
parameters, viz., distribution of nodes, overlap, definition of spatial basis functions
and temporal basis functions, etc, are identical to those used earlier. The L2 error
norm is evaluated at 13 x 13 x 13 uniformly distributed sample points. The relative
error is 1.1e—2. ?? shows analytical and numerical results at location (0.3, 0.3, 0.3).
As shown, the result is good.

Next, iterative solver is implemented in the time domain GFEM solver. Consid-
ering a homogeneous rectangular computational domain 9 = (0,1)2 with boundary
1‘ = Ugr, where {1‘1 ; r e 0 x(0,1)}, {1"2:r 6 (0,1) x 1}, {r3 : r e 1 x (0,1)},
{1‘4 : r 6 (0,1) x 0}. In this domain, it is assumed that the electric field is of the
form

E(r,t) = {12(t + 55—2330)/726$p{_(t + 2299/72} (4.28)

where 11:0 = 5.5 and T = 3.15 x10‘93. The solution to (4.1) is solved by imposing the
analytical solution as Neumann boundary condition. In this experiment, as before,
all nodes are distributed uniformly, and rectangular patches are used. The size of
each patch is chosen to be 1.5 times the distance between the nodes. The weight
functions W, (r) is a tensor product rooftops, and the local approximation function
vim (r) and vi=n(r) are tensor product of Legendre polynomials. The parameters
of spatial basis flmctions are h = 0.38A and p = 3. APSW functions are used for
as temporal basis. The half width of prolate function p and the number of sample

points for extrapolation N samp are 6 and 4, respectively. The time step size is

137

At = 2.33 x 10‘113. Iterative solver is used and

llZlofllg " {WI < 16—8 (4.29)
|{V}|

 

is used as the ending criteria of the iterative solver. Figure 4.14 show the numerical
results with relative error 1.4—3. Figure 4.15 show the number of iterative steps
in each time step. From the figure we can see that at the first few time steps, the
numbers of iterative steps are large and then decrease as expected. For the first few
time steps, the error of initial guess is large. When the initial guess is getting more

and more accurate, the number of iterative steps decreases.

138

 

 

 

 

 

 

 

 

 

0 500 1000 1500
Time step

Figure 4.14. Error with iterative solver

 

4O . w

0‘
0|
1

8

N
U!
L

 

.3
0|

 

Number of iterative steps
N
O

.5
O

 

 

 

Index of time step

Figure 4.15. Number of iterative steps in each time step

139

CHAPTER 5

CONCLUSION AND AREAS FOR FUTURE WORK

Generalized finite element method (FEM) is a generalization of the classical finite
element method. As a partition of unity-based solver, GFEM solver simplifies mesh
generalization and enlarges the local approximation function space. At the same
time, it maintains the h—, p— and hp— convergence in the computation domain.
Vector GFEM solver can be developed in inhomogeneous domain to simulate the
field which is divergence free and satisfy the continuity requirements on the material
interface. This research has proved the convergence of the solver mathematically and
demonstrated it in two dimension and three dimension simulations. We can solve
wave equation with every boundary conditions encountered in the EM simulation
for close domain and open domain. We also have developed an extension of GFEM
solver to simulate the transient field. SVD-based preconditioner is utilized to solver
the problem of the high condition number of the linear system when the high order
numerical solver is utilized.

A summary of most important contributions in this research follows:

0 Implementation of N itsche’s method to impose the essential boundary condi-
tion in GFEM solver. Convergence of the scalar and vector GFEM solver with

this method is demonstrated in both 2D and 3D simulation.

0 Development of the first hybrid GFEM-BI technique that enables the trun-
cation of the computational domain, and demonstrated its h, p convergence
properties. We have specified the manner in which different types of boundary
integrals may be incorporated into the simulations. Validating the proposed

technique for scattering from more complex shapes necessitated the develop-

ment of GFEM-PML techniques.

140

0 Development of an extension of CF EM to enable application to vector electro-
magnetic problems in piecewise homogeneous domain with smoothly curved
interface. Basis functions can simulate the divergence free field which satisfy
all the continuity requirements on the material interface. The convergence of
this spurious-free solver is proved mathematically and demonstrated in 2D and

3D simulation.

0 Application of GFEM solver to the problem where the geometry has singular-
ity. The eigenfunction of the singular structure is used in local approximation.
The advantage of GFEM solver is demonstrated compared with classical FEM

solver.

0 Development of time domain GFEM method for solving time domain differ-
ential equations. The method uses different kinds of high order time stepping
scheme including Approximate prolate spheroidal wave functions. The ac-
curacy of the solver is demonstrated in EM simulation in both 2D and 3D,

homogeneous and inhomogeneous domain.

0 Stability analysis of general high order time stepping solver are proposed and

demonstrated with APSW-based time stepping schemes.

0 Development of the svd-based preconditioner to high order numerical solver.
The error bound and cost of this method is derived. We demonstrated that
when the condition of the resulting system is improved dramatically, the ac-

curacy of the solver are almost the same.

It also found that GFEM solver has the following advantage/ disadvantage com-
pared with classical FEM solver

Advantage:

e The mesh generation of GFEM is simpler than the classical FEM. The nodes’s

location and patches’s size and shape can be chosen randomly. The edge of

141

patch is not necessary to conformal to the boundary of computation domain.

0 A broader class of functions can be implemented to simulate the field in GFEM.
When the same accuracy is required, GFEM solver results in smaller linear
system. When simulating the fields, regular basis function in classic FEM is
not as efficient as special basis functions in GFEM such as the eigenfunction-
based basis function around the geometric singularity and the planewave-based
basis function in a very smooth region. Though in FEM, singular basis function
is also developed around the geometric singularity, there is limitation in the

size of the mesh of singular basis function.
Disadvantage:

0 It’s more expensive to integrate the basis function and its derivative in GFEM
when evaluating the matrix elements. Compared with the analytical integra-
tion in FEM, numerical integration in GFEM results in an additional burden.
But the burden can be elevated when the nodes are distributed uniformly with

the same shape and size of patches.

0 Implementation of geometry information is more expensive in GFEM. In FEM,
implementation of information of mesh is straightforward. In GFEM, more

procedure is necessary to utilize the mesh information to GFEM solver.

Until now, vector GFEM solver has been developed for electromagnetic simula-
tion with simple structure. Application of the technology to realistic problem implies

the further research along the following which is necessary:

0 Planewave—based basis function is much more efficient than the polynomials in
the region where fields are smooth which has been demonstrated by two dimen-
sion homogeneous results. Augmenting the function space within planewave

may be practical.

142

0 When vector GF EM solver for curved material interface has been developed
and demonstrated in 2D problem. Convergence of this solver should also be

demonstrated in 3D.

0 Vector GFEM solver has been developed for piecewise homogeneous domain
with smooth curved material interface. Considering the material interface with
geometric singularity such as vertex or edge, appropriate local approximation
function has not been developed. This is a challenging problem. Based on the
continuity requirements, the additional function for this material interface can

be constructed as the eigenfunction of PEC with the same structure.

0 Given meshes that describes the structure of the electromagnetic problem, the
geometric information easily utilized by GFEM solver should be extracted from
the mesh information. Those information should be able to be implemented

by GFEM solver automatically.

0 High order integration in the computation domain, e.g. volume integration,
and on the boundary, e.g. surface integration, need to be implemented for
high order solver. The quadrature rules of such integration can be derived by

geometric transformation.

143

APPENDICES

144

Vita

Chuan Lu was born in Tianjin, BBC, on September 16, 1979. He received the
BS. degree in electrical engineering from Tsinghua University, Beijing in 2002. He is
currently pursuing the Ph.D. degree in the Department of Electrical and Computer
Engineering at the Michigan State University. Since 2002, he has been a Research
Assistant at the Michigan State University. His research interests are in area of com-
putational electromagnetics, especially differential equation-based numerical solver
such as generalized finite element method.

A list of Chuan Lu’s recent articles is given below:

Journal papers:

1. C. Lu and B.Shanker, “Analysis of complex structures using Generalized Fi-
nite Element Method,” IEEE Transactions on Antenna and Propagations, In
preparation.

2. O. 'Ihncer, C. Lu, B. Shanker and L. C. Kempel, “Dispersive Analysis of
Vector Generalized Finite Element Method,” IEEE Transactions on Antenna
and Pmpagations, In preparation

3. C. Lu and B.Shanker, “Transient Analysis of Electromagnetic Fields using
the Generalized Finite Element Methods,” IEEE Transactions on Antenna

and Prapagations, Submitted.

4. C. Lu and B.Shanker, “Generalized Finite Element Method for Vector Electro-
magnetic Problems,” IEEE Transactions on Antenna and Propagations, vol.
55, Issue 1, pp. 1369-1381, 2007.

5. C. Lu and B.Shanker, “Development of Hybrid Boundary Integral-Generalized
(Partition of Unity) Finite Element Solvers for the scalar Helmholtz Equation,”
IEEE Transactions on Magnetics, vol. 43, Issue 3, March 2007 Page(s):1002 -
1012

6. Haixiang Cao, Lv Chuan, Jianguo Jiang and Tamasi, “Analysis of Power Loss
and Choice of Parameter to over voltage-limited filter for PWM AC Inverter
controlled motor,” Transactions of China Electrotechnical Society,2003

Conference papers:

1. O. Tuncer, C. Lu, N. V. Nair, B. Shanker and L. Kempel, “Analysis of Error
Propagation in Vector Generalized Finite Element Methods,” Intemation Con-

ference on Elecromagnetics on Advanced Applications, Torino, Italy, Septem-
ber 2007

2. N. V. Nair, C. Lu, B. Shanker, “Differential and Integral Equation Solvers
based on Generalized Moments and Partitions of Unity,” Intemation Confer-

ence on Elecr‘omagnetics on Advanced Applications, Torino, Italy, September
2007

145

03

10.

11.

12.

13.

C. Lu and B.Shanker, “Analysis of complex structures using Generalized
Finite Element Method,” North American Radio Science Meeting URSI-
CNC/USNC’, Ottawa, ON, Canada, 2007

C. Lu and B.Shanker, “Transient Analysis of Electromagnetic Fields using the
Generalized Finite Element Methods,” North American Radio Science Meeting
URSI-C'N C/ USN 0, Ottawa, ON, Canada, 2007

O. Tuncer, C. Lu, B. Shanker and L. C. Kempel, “Dispersive Analysis of
Vector Generalized Finite Element Method,” North American Radio Science
Meeting URSI-CNC/USNC, Ottawa, ON, Canada, 2007

C. Lu and B.Shanker, “Development of Generalized Finite Element Method
for Vector Electromagnetic Problems,” IEEE AP-S International Symposium
on Antennas and Propagation and USNC/URSI, Albuquerque, New Mexico,
2006

C. Lu, J. Villa and B. Shanker, “Application of Generalized Finite Element
Method to Special Scattering Problem,” IEEE AP—S International Symposium
on Antennas and Propagation and USNC/URSI, Albuquerque, New Mexico,
2006

C. Lu, B. Shanker, E. Michielssen, “Development of "lime Domain General-
ized Finite Element Method for Vector Electromagnetic Problems,” IEEE AP-
S International Symposium on Antennas and Propagation and USNC/URSI,
Albuquerque, New Mexico, 2006

H. Huang, C. Lu and B. Shanker, “Cartesian Harmonics and Fast Compu-
tation of Potentials of the Form,” IEEE AP-S International Symposium on
Antennas and Propagation and USNC/URSI, Albuquerque, New Mexico, 2006

Z. Zeng, K. Arunachalam, C. Lu, B. Shanker, and L. Udpa, “Element-free
Galerkin method in modeling microwave inspection of civil structures,” Di-
gest Book of the 12th Biennial IEEE Conference on Electromagnetic Field
Computation- CEFC', Florida, USA, Apr.-May 2006

C. Lu and B.Shanker, “Development of Hybrid Boundary Integral-Generalized
(Partition of Unity) Finite Element Solvers for Scalar Problems,” IEEE AP-

S International Symposium on Antennas and Propagation and USNC/URSI,
Washington, DC, 2005

C. Lu and B.Shanker, “Solving Boundary Value Problems Using the General-
ized (Partition of Unity) Finite Element Method,” IEEE AP-S International
Symposium on Antennas and Propagation and USN C/ URSI, Washington, DC,
2005

C. Lu and B.Shanker, “Towards the Development of Vector Generalized (Par-
tition of Unity) Finite Element Method,” IEEE AP—S International Sympo-
sium on Antennas and Propagation and USNC/URSI, Washington, DC, 2005

146

14.

15.

L. Chuan and B.Shanker, and L. Kempel, “Partition of Unity methods in
electromagnetics,” Proceedings of the URSI-EM Theory Symposium, p. 331,
2004.

Deng Wen; Lv Chuan; Jiang Jianguo; Huang Lipei; Tamai Shinzo, “Filter
design of restraining over-voltage and bearing current in induction motor fed
by PWM inverter,” TENCON ’02. Proceedings. IEEE Region 10 Conference
on Computers, Communications, Control and Power Engineering, 2002

147

Bibliography

[1] J. M. Jin, The Finite Element Method in Electromagnetics. New York: Wiley,
2002.

[2] I. H. Shames and C. L. Dym, Energy and finite element methods in structural
mechanics. Second edition. Hemisphere Publishing, New York, NY, 1984.

[3] D. K. Gartling and J. N. Reddy, The Finite Element Method in Heat Transfer
and Fluid Dynamics. CRC Press, 2001.

[4] A. Konrad, “Vector variational formulation of electromagnetic fields in
anisotropic media,” IEEE Trans. Microwave Theory Tech, vol. MTT-24,
pp. 553-559, 1976.

[5] B. M. A. Rahman and J. B. Davies, “Penalty function improvement of waveg-
uide solution by finite elements,” IEEE Trans. Microwave Theory Tech,
vol. MTT-32, pp. 922-928, 1984.

[6] J. C. Nedelec, “Mixed finite elements in T3,” Numer. Meth., vol. 35, pp. 315—
341, 1980.

[7] M. L. Barton and Z. J. Cendes, “New vector finite elements for three
dimensional magnetic field computation,” J. Appl. Phys, vol. 61, pp. 3919—
3921, 1987.

[8] J. P. Webb and B. Forghani, “Hierarchal scalar and vector tetrahedra,” IEEE
Trans. Magn., vol. MAG-29, pp. 1495—1498, 1993.

[9] Z. J. Cendes, “Vector finite elemetns for electromagnetic field,” IEEE Trans.
Magn., vol. MAG-27, pp. 3958-3966, 1991.

[10] R. D. Graglia, D. R. Wilton, and A. F. Peterson, “Higher order interpolatory
vector bases for computational electromagnetics,” IEEE Transactions on An-
tennas and Propagation, vol. 45, pp. 329—342, 1997.

[11] W. Rachowicz and L. Demkowicz, “A three-dimensional hp—adaptive finite ele—
ment package for electromagnetics.” TICAM Report 00-04, 2000.

[12] M. Costabel, M. Dague, and C. Schwab, “Exponential convergence of hp—fem for
maxwell’s equations with weighted regularization in polygonal domains,” Math-
ematical Models and Methods in Applied Sciences, vol. 15, 2005. To Appear: pa-
per available at http://perso.univ-rennes1.fr/monique.dauge/core/index.html.

148

[13] D. Boffi, M. Costabel, M. Dague, and L. Demkowicz, “Discrete compactness
for the hp version of rectangular edge finite elements.” ICES Technical Report
04-29, 2004.

[14] D.-K. Sun, L. Vardapetyan, and Z. Cendes, “Two-dimensional curl-conforming
singular elements for fem solutions of dielectric waveguiding structures,” IEEE
Transactions on Microwave Theory and Techniques, vol. 53, pp. 984—992, 2005.

[15] R. D. Graglia and G. Lombardi, “Singular higher order complete vector bases
for finite methods,” IEEE Transactions on Antennas and Propagation, vol. 52,
pp. 1672—1685, 2004.

[16] J. Jirousek and A. Venkatesh, “Hybrid-trefftz plane elasticity elements with
p—method capabilities,” Int. J. Num. Meth. Eng, vol. 35, pp. 1443—1472, 1992.

[17] J. S. Hesthaven and T. Warburton, “Nodal high-order methods on unstructured
grids,” Journal of Computational Physics, vol. 181, pp. 186—221, 2002.

[18] I. Babuska and J. M. Melenk, “The partition of unity method,” International
Journal for Numerical Methods in Engineering, vol. 40, pp. 727—858, 1997.

[19] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl, “Meshless
methods: An overview and recent developments,” Comp. Meth. Appl. Mech.,
1996.

[20] L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,”
The Astronomical Journal, vol. 82(12), pp. 1013—1024, 1977.

[21] G. T. B. Nayroles and P. Villon, “Generalizing the finite element method:diffuse
approximation and diffuse elements,” Computational Mechanics, vol. 10,
pp. 307—318, 1992.

[22] Y. Y. L. T. Belytschko and L. Gu, “Element-free galerkin methods,” Inter-
national Journal for Numerical Methods in Engineering, vol. 37, pp. 229—256,
1994.

[23] I. Babuska, U. Banerjee, and U. Osborn, Meshless and generalized FEM: A
survey of some major results. Springer Verlag, 2002.

[24] C. A. M. Duarte and J. T. Oden, “Hp clouds—an hp meshless method,” Numer-
ical Methods for Partial Differential Equations, vol. 12, pp. 673—705, 1996.

[25] S. N. Atluri and S. P. Razek, The Meshless Local Petrov-Galerkin{MLPG)
Method. Tech Science Press, 2002.

[26] L. Xuan, B. Shanker, and L. Udpa, “A meshless element free galerkin method for
nde applications,” in Proceedings of Appliced Computational Electromagnetics
Society, 2002.

149

[27] L. Xuan, B. Shanker, Z. Zeng, and L. Udpa, “Element-free galerkin method in
pulsed eddy currents,” The International Journal of Applied Electromagnetics
and Mechanics, vol. 19, pp. 463—466, 2004.

[28] L. Xuan, Z. Zeng, B. Shanker, and L. Udpa, “Development of a meshless finite
element model for nde applications,” IEEE Transactions on Magnetics, vol. 40,
pp. 12—20, 2004.

[29] L. Xuan, Z. Zheng, B. Shanker, and L. Udpa, “Meshless method for numerical
modeling of pulsed eddy currents,” IEEE Transaction on Magnetics, vol. 40,
pp. 3457 — 3462, 2004.

[30] L. Chuan, B. Shanker, and L. Kempel, “Partition of unity methods in electro-
magnetics,” in Proceedings of the URSI-EM Theory Symposium, p. 331, 2004.
See the technical report series @ECE, Michigan State University for a full ver-
sion of the paper or the lead PIs website.

[31] Z. Zeng, L. Xuan, B. Shanker, L. Kempel, and L. Udpa, “Development of vector
basis function for meshless efg method.” Presented at URSI, 2003, 2003.

[32] I. Tsukerman, “A flexible local approximation method for electro— and magne-
tostatics,” IEEE Transactions on Magnetics, vol. 40, pp. 941 — 944, 2004.

[33] A. Plaks, I. Tsukerman, and G. F. B. Yellen, “Generalized finite-element
method for magnetized nanoparticles,” IEEE Transactions on Magnetics,
vol. 39, pp. 1436 — 1439, 2003.

[34] I. Tsukerman, “Spurious numerical solutions in electromagnetic resonance prob-
lems,” IEEE Transaction on Magnetics, vol. 39, pp. 1405 — 1408, 2003.

[35] I. Tsukerman, “General tangentially continuous vector elements,” IEEE Trans-
actions on Magnetics, vol. 39, pp. 1215 — 1218, 2003.

[36] L. Proekt, S. Yuferev, I. Tsukerman, and N. Ida, “Method of overlapping
patches for electromagnetic computation near imperfectly conducting cusps and
edges,” IEEE Transaction on Magnetics, vol. 38, pp. 649 - 652, 2002.

[37] L. Proekt and I. Tsukerman, “Method of overlapping patches for electromag-
netic computation,” IEEE Transaction on Magnetics, vol. 38, pp. 741 — 744,
2002.

[38] Z. Zeng, B. Shanker, and L. Udpa, “Modeling microwave nde using the element
free galerkin method—development of perfectly matched layers for domain trun-
cation,” in Proceedings of the Tenth International Workshop on Electromagnetic
Non-destructive Evaluation, 2004.

[39] T. F. M. Hara, T. Wada and F. Kikuchi, “A three dimensional analysis of
rf electromagnetic fields by the finite element method,” IEEE Trans. Mag-19,
vol. 6, pp. 2417-2420, 1983.

150

[40] J. P. Webb, “The finite element method for finding modes of dielectric-loaded
cavities,” IEEE Trans. M TT—33, vol. 7, pp. 635—639, 1985.

[41] I. Bardi and O. Biro, “An efficient finiteelement formulation without spurious
modes for anisotropic waveguides,” IEEE Transactions on Microwave Theory
and Techniques, vol. 39, pp. 1133—1139, 1991.

[42] M. E. K. Hayata, M. Koshiba and M. Suzuki, “Vectorial finite element method
without any spurious solutions for dielectric waveguiding problems using trans-
verse magnetic field component,” IEEE Trans. M TT-34, vol. 11, pp. 1120—1124,
1986.

[43] A. Konrad, “On the reduction of the number of spurious modes in vectorial
finite element solution of three-dimensional cavities and waveguides,” IEEE
Trans. M TT-34, vol. 2, pp. 224—227, 1986.

[44] S. N. Atluri and S. Shen, The Meshless Local Petrov Galerkin Method. Tech.
Sci. Press, 2002.

[45] J. M. Melenk and I. Babuska, “The partition of unity finite element method:
Basic theory and applications,” Comput. Methods Appl. Mech. Engr., vol. 139,
pp. 289—314, 1996.

[46] G. C. Hsiao and R. E. Klienmann, “Mathematical foundations for error esti-
mation in numerical solutions of integral equations in electromagnetics,” IEEE
Trans. Antennas Propag., vol. 45, pp. 316—328, 1997.

[47] M. Griebel and M. A. Schweitzer, “Particle-partition of unity methods for
pdes,” SIAM J. Sci. Compt, vol. 22, pp. 853—890, 2002.

[48] D. Shepard, “A two-dimensional interpolation function for irregularly spaced
data,” in Proceeding of ACM, p. 517?24, 1968.

[49] C. A. Duarte, “The hp cloud method,” Dissertation, 1996.

[50] M. Griebel and M. A. Schweitzer, Meshless Methods for Partial Difl‘erential
Equations. Springer Verlag, 2002.

[51] S. Mendez and A. Huerta, “Imposing essential boundary conditions in mesh-free
methods,” Computer Methods in Applied Mechanics and Engineering, vol. 193,
pp. 1257—1275, 2004.

[52] J. Pitkdranta, “Boundary subspaces for the finite element method with lagrange
multipliers,” Numer. Math, vol. 33, pp. 273—289, 1979.

[53] J. Pitkdranta, “Local stability conditions for the babuska method of lagrange
multipliers,” Math. Comput., vol. 35, pp. 1113-1129, 1980.

[54] J. Pitkc’iranta, “The finite element method with lagrange multiplier for domain
with corners,” Math. Comput., vol. 37, pp. 13—20, 1981.

151

[55] R. Stenberg, “On some techniques for approximating boundary conditions in the
finite element method,” Journal of Computational and Applied Mathematics,
vol. 63, pp. 139—148, 1995.

[56] L. F. Canino, J. J. Ottusch, M. A. Stalzer, J. L. Visher, and S. M. Wandzura,
“Numerical solution of the helmholtz equation in 2d and 3d using a high-order
nystrm discretization,” J. Comput. Phys, vol. 146, pp. 627—663, 1998.

[57] W. C. Chew, Waves and fields in inhomogeneous media. New York: IEEE
Press, 1995.

[58] W. C. Chew and W. H. Weedon, “A 3d perfectly matched medium from mod-
ified maxwell’s equations with stretched coordinates,” Microw. Opt. Technol.
Lett., vol. 7, pp. 599—604, 1994.

[59] Y. Saad, Iterative Methods for Sparse Linear Systems. New York: PWS Pub—
lishing Company, 1996.

[60] V. Rokhlin, “Rapid solutions of integral equations of scattering theory in two
dimensions,” J. Comput. Phys, vol. 86, pp. 414—439, 1990.

[61] R. F. Harrington, Time Harmonic Fields in Electromagnetics. IEEE Press
Series on Electromagnetic Wave Theory, Wiley-IEEE Press, 2001.

[62] J. A. Stratton, Electromagnetic Theory. New York, New York: McGraw-Hill,
1941.

[63] C. Lu and B. Shanker, “Development of hybrid boundary integral-generalized
(partition of unity) finite element solvers for the scalar helmholtz equation,”
technical report, Michigan State University, 2005. Submitted to IEEE Magnet-
ics.

[64] J. A. Chatterjee and J. L. Volakis, “Computation of cavity reasonances us-
ing edge-based finite elements,” IEEE Transactions on Microwave Theory and
Techniques, vol. 40, pp. 2106—2108, 1992.

[65] K. S. Kunz and R. J. Luebbers, The Finite-Diflerence Time-Domain Method
in Electromagnetics. CRC, 1993.

[66] L. E. R. Petersson and J. M. Jin, “Analysis of periodic structures via a time-
domain finite-element formulation with a floquet abc,” IEEE Transactions on
Antennas and Propagation, vol. 54, pp. 933—944, 2006.

[67] Z. Lou and J. M. Jin, “A novel dual-field time-domain finite-element domain-
decomposition method for computational electromagnetics,” IEEE Transac-
tions on Antennas and Propagation, vol. 54, pp. 1850—1862, 2006.

[68] Z. Lou and J. M. Jin, “Modeling and simulation of broad-band antennas using

the time-domain finite element method,” IEEE Transactions on Antennas and
Propagation, vol. 53, pp. 4099—4110, 2006.

152

[69] D. S. Weile, G. Pisharody, N.-W. Chen, B. Shanker, and E. Michielssen, “An
exponentially convergent scheme for the solution of the time-domain integral
equations of electromagnetics,” IEEE Trans. Antennas Propagat, 2004.

[70] J. J. Knab, “Interpolation of band-limited functions using the approximate
prolate series,” IEEE Trans. Inform. Theory, vol. IT-25, pp. 717—720, 1979.

[71] C. Lu and B. Shanker, “Generalized finite element method for vector electro-
magnetic problems,” tech. rep., Michigan State University, 2005.

[72] J. A. Cadzow, “Extrapolation procedure for band-limited signals,” IEEE Trans.
Acoust, Speech, Singal Processing, vol. 27, pp. 4—12, 1979.

[73] C. Lu and B. Shanker, “Generalized finite element method for vector electro-
magnetic problems,” IEEE Transaction on Antennas and Propagation, vol. 55,
pp. 1369—1381, 2007.

153

   

[lllllllllllill]