GENERALIZED METHOD OF MOMENTS: A NOVEL DISCRETIZATION
TECHNIQUE FOR INTEGRAL EQUATIONS
By
Naveen V. Nair

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY
Electrical Engineering
2010

ABSTRACT
GENERALIZED METHOD OF MOMENTS: A NOVEL
DISCRETIZATION TECHNIQUE FOR INTEGRAL EQUATIONS
By
Naveen V. Nair
Integral equation formulations to solve electromagnetic scattering and radiation problems have existed for over a century. The method of moments (MoM) technique to
solve these integral equations has been in active use for over 40 years and has become
one of the cornerstones of electromagnetic analysis. It has been successfully employed
in a wide variety of problems ranging from scattering and antenna analysis to electromagnetic compatibility analysis to photonics. In MoM, the unknown quantity
(currents or fields) is represented using a set of basis functions. This representation,
together with Galerkin testing, results in a set of equations that may then be solved
to obtain the coefficients of expansion. The basis functions are typically constructed
on a tessellation of the geometry and its choice is critical to the accuracy of the final solution. As a result, considerable energy has been expended in the design and
construction of optimal basis functions. The most common of these functions in use
today are the Rao-Wilton-Glisson (RWG) functions that have become the de-facto
standard and have also spawned a set of higher order complete and singular variants.
However, their near-ubiquitous popularity and success notwithstanding, they come
with certain important limitations.
The chief among these is the intimate marriage between the underlying triangulation of the geometry and the basis function. While this coupling maintains continuity of the normal component of these functions across triangle boundaries and
makes them very easy to implement, this also implies an inherent restriction on the
kind of basis function spaces that can be employed. This thesis aims to address this
issue and provides a novel framework for the discretization of integral equations that

demonstrates several significant advantages.
In this work, we will describe a new umbrella framework for the discretization
of integral equations called the Generalized Method of Moments (GMM). We will
show that within this framework it is possible to choose from a very wide variety
of functions to form the basis space. While the choice of basis functions can be
completely arbitrary, the thesis will describe in detail one particular scheme that
utilizes a surface quasi-Helmholtz decomposition and provides a mechanism for the use
of any available knowledge of the physics of the problem. The thesis will describe the
theoretical framework required for the design of this scheme and develop error bounds
for the approximation of the unknown quantities using the proposed basis functions.
Following that, we will delve into the implementation of the GMM starting from
commonly available triangular tessellations, carefully identifying and solving several
key issues along the way.
The GMM, developed here, will be applied to a plethora of examples, validated
against analytical solutions and compared against the standard RWG discretizations.
It will further be shown that the use of the quasi-Helmholtz decomposition will result
in a scheme that is well conditioned over a wide range of frequencies. We will then
extend the GMM to both the analysis of scattering from dielectric bodies using the
PMCHWT and M¨ller operators and to the analysis of transient scattering from
u
conducting objects. Several results will be presented in either case that demonstrate
the advantage of the proposed method.

DEDICATION

To my wife, my parents and my sister

iv

ACKNOWLEDGMENT
This thesis would have been impossible without the unfaltering support and advice of
Shanker Balasubramaniam, my primary advisor. His understanding and enthusiasm
has helped me navigate many a tough corner during the last four years. He has taught
me so much more than computational electromagnetics that it would be an injustice
to even try to list it here. He has gone, and continues to go, far beyond the calls of
duty as an advisor to be a friend, guide and reality check to me and his focus, energy
and work ethic will be constant sources of inspiration throughout my life and career.
During the first three years of my graduate career I was supported and guided by
my co-advisor Satish Udpa. The experience I gained under his tutelage will prove
invaluable to me in the future. I have been indeed blessed with two excellent PhD
advisors; either of whom any student would be lucky to get alone. I would also like to
thank each member of my doctoral committee: Leo Kempel, Ed Rothwell and Jian
Liang Qian for their support and guidance.
Several people have played key roles behind the curtain to make my graduate
student life work. Chief among them is Linda Clifford, who has patiently proofread many a document for my research and cooked many a meal for thanksgiving
dinners throughout the seven years I have been in graduate school. I would also
like to acknowledge Pauline Vandyke for all her support when I neeeded a teaching
assistanship to push through a semester or when I needed to steal some coffee from
the faculty kitchen.
I would like to acknowledge the support of my friends and lab mates: Vikram
Melapudi, Andrew B., Andrew P., Ozgur Tuncer and Dan Dault for being excellent
touchstones for my ideas over the past several years. In particular I would like to
thank Andrew Pray for generating several of the results in the latter part of this
thesis. A special thanks needs to go out to Vikram Melapudi and Mayur Mudigonda
for their invaluable friendship over the years.
v

With more than two decades of schooling nearly behind me, I would like to thank
my teachers over the years who have played a major part in making me who I am today: Poonam Lal, Nakku and Kusum Chandragupta, back in highschool, for teaching
me how to learn; A. Venkatesh at the Indian Institute of Technology Madras for first
getting me interested in heat transfer, and therefore, partial differential equations and
finally, Krishnan Balasubramaniam for being the first advisor I ever had and instilling
in me a deep rooted passion for rigor in reasearch.
I would like to thank my brother in law, Saju, for appearing in my family’s life
exactly around the time of my comprehensive and distracting my entire family enough
to give me an extra year to finish this thesis.
And, finally my deepest gratitude goes out to my family: my wife Aneesha, my
mother Sobhana, father Vijayakumaran and my sister Indu whose constant and unwavering support has helped me keep my sanity through my (long) tenure as a graduate
student. It is to them that I dedicate this thesis.

vi

TABLE OF CONTENTS

List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

x

1 Introduction . . . . . . . . . . . . . . . . . . . . .
1.1 Differential Equation Methods . . . . . . . . . . . . . . .
1.2 Integral Equation Methods for Electromagnetics . . . . .
1.3 Basis Functions for Electromagnetic Integral Equations .
1.4 Higher order basis functions . . . . . . . . . . . . . . . .

.
.
.
.
.

.
.
.
.
.

. .
. . .
. . .
. . .
. . .

.
.
.
.
.

.
.
.
.
.

1
2
3
4
6

2 Mathematical Preliminaries . . . . . . . . .
2.1 Maxwell’s Equations . . . . . . . . . . . . . . .
2.2 Vector and Scalar Potentials . . . . . . . . . . .
2.3 Integral Equation Representation of Solutions .
2.4 EFIE and MFIE for PEC structures . . . . . . .

.
.
.
.
.

.
.
.
.
.

. .
. . .
. . .
. . .
. . .

.
.
.
.
.

.
.
.
.
.

10
10
13
16
20

3 Numerical Solution of Integral Equations . . . . . . . . . . . .
3.1 Method of moments solution for integral equations . . . . . . . . . . .
3.2 Approximation errors . . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23
25

4 The
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8

30
33
35
42
46
51
52
52
56
60
60

4.9

.
.
.
.
.

.
.
.
.
.

. .
. . .
. . .
. . .
. . .

Generalized Method of Moments . . . . . . . . . . . . . .
Partition of Unity functions and Overlapping Patches . . . . . . . . .
Definition of Basis Functions . . . . . . . . . . . . . . . . . . . . . . .
Design of the Vector Approximating Function . . . . . . . . . . . . .
Error Bounds for the Discrete Integral Operators . . . . . . . . . . .
Benefits of the GMM framework . . . . . . . . . . . . . . . . . . . . .
Implementation of the GMM . . . . . . . . . . . . . . . . . . . . . . .
Construction of Overlapping Patches . . . . . . . . . . . . . . . . . .
Computation of Matrix Elements . . . . . . . . . . . . . . . . . . . .
r and r lie on the same triangle . . . . . . . . . . . . . . . . .
r and r lie on triangles that share an edge . . . . . . . . . . .
r and r are not in neighbouring triangles but are close (|r−r | ≤
0.15λ) . . . . . . . . . . . . . . . . . . . . . . . . . .
r and r are sufficiently separated from each other . . . . . . .
A transformation for the near singular integration of
g(R)j(r ) for r → r . . . . . . . . . . . .
4.8.1 r and r are not in neighbouring triangles but are close (|r−r | ≤
0.15λ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.8.2 r and r are sufficiently separated from each other . . . . . . .
4.8.3 Integration with acceleration techniques . . . . . . . . . . . .
Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . .
vii

60
60
60
61
61
62
62

4.10 Quality of Approximation . . . . . . . . . . . . . . . .
4.10.1 Approximation of Polynomial Functions . . . .
4.10.2 Approximation of sinusoidal functions . . . . . .
4.10.3 Approximation of non-smooth and near-singular
4.11 Scattering Results . . . . . . . . . . . . . . . . . . . . .
4.11.1 Analytical Validation and p-convergence . . . .
4.11.2 Comparison with standard RWG . . . . . . . .
4.11.3 Comparison with measured data . . . . . . . . .
4.12 Flexibility of the GMM scheme . . . . . . . . . . . . .
4.12.1 Mixing orders of basis functions . . . . . . . . .
4.12.2 Mixing types of basis functions . . . . . . . . .
4.12.3 Use of non-conformal meshes . . . . . . . . . .
4.13 Condition Number of the System . . . . . . . . . . . .
4.14 Summary . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .
. . . . . .
. . . . . .
functions
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.

62
63
64
64
69
70
71
74
78
78
79
83
83
91

5 Application of the Generalized Method of Moments for Dielectric
Bodies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.2 GMM Discretization of the PMCHWT and M¨ller Integral Equations 97
u
5.2.1 M¨ller Operator . . . . . . . . . . . . . . . . . . . . . . . . . . 97
u
5.3 Implementation Results . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.3.1 Validation of GMM PMCHWT and M¨ller Solvers . . . . . . 99
u
5.3.2 Salient features of the GMM scheme . . . . . . . . . . . . . . 100
Mixtures of Basis functions . . . . . . . . . . . . . . . . . . . 105
Mixed Order Bases: . . . . . . . . . . . . . . . . . . . . 105
Mixed Types of Bases: . . . . . . . . . . . . . . . . . . 106
hp-adaptivity: . . . . . . . . . . . . . . . . . . . . . . . 108
Non Conformal meshes . . . . . . . . . . . . . . . . . . . . . . 111
5.3.3 Low frequency stability and multiply connected scatterers . . 111
6 Analysis of transient scattering using the Generalized Method of
Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.1 The Fourier transform pair . . . . . . . . . . . . . . . . . . . . . . . . 118
6.2 Time Domain EFIE and MFIE . . . . . . . . . . . . . . . . . . . . . 119
6.3 The Marching on in Time (MOT) scheme . . . . . . . . . . . . . . . 122
6.4 Temporal basis functions . . . . . . . . . . . . . . . . . . . . . . . . . 124
6.5 GMM for transient scattering . . . . . . . . . . . . . . . . . . . . . . 126
6.6 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . 129
6.6.1 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
6.6.2 Low frequency stability . . . . . . . . . . . . . . . . . . . . . . 130
6.6.3 Mixtures of basis functions . . . . . . . . . . . . . . . . . . . . 135
6.6.4 Non conformal meshes . . . . . . . . . . . . . . . . . . . . . . 138
viii

7 Conclusions and Future Work . . . . . . . . . . . . . . . . . . 143
7.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
7.1.1 Mathematical proof for the well conditioned nature of the GMM 145
7.1.2 Generalized higher order patch construction . . . . . . . . . . 146
7.1.3 Further development of the transient GMM . . . . . . . . . . 147
7.1.4 Application to multi-scale geometries . . . . . . . . . . . . . . 147
7.1.5 Domain Decomposition . . . . . . . . . . . . . . . . . . . . . . 148
7.1.6 Application to other integral equations . . . . . . . . . . . . . 148
Bibliography

. . . . . . . . . . . . . . . . . . . . . . . . . . 150

ix

LIST OF FIGURES

General representation for a domain D+ = V with inner and outer
boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2

General description of a scattering problem . . . . . . . . . . . . . . .

20

3.1

A schematic showing the definition of the RWG basis functions on an
edge shared by two triangles . . . . . . . . . . . . . . . . . . . . . . .

28

Currents on a PEC cube showing the singular nature of the currents
near the edges and the smooth currents towards the center. For interpretation of the references to color in this and all other figures, the
reader is referred to the electronic version of this dissertation. . . . .

32

4.2

Construction of GMM patches (Ωi ) shown as shaded region. . . . . .

34

4.3

Examples of linear PU functions in (a) one and (b) two dimensions .

36

4.3

Examples of linear PU functions in (a) one and (b) two dimensions .

37

4.4

Construction of scalar GMM basis functions. (a) the GMM approximation function and (b) the product of the approximation function
and the partition of unity function. . . . . . . . . . . . . . . . . . . .

38

Vector GMM basis functions based on the quasi-surface Helmhotlz decomposition basis functions. (a) The quasi curl-free basis functions
and (b) the quasi-divergence free basis functions. Both functions are
of order p = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

44

Construction of GMM patches (Ωi ) shown as shaded region. The
patches are constructed as a set of triangles joined at a node. . . . . .

52

2.1

4.1

4.5

4.6

x

4.7

Construction of (a) GMM Patches (Ωi ) shown as shaded region and
(b) Projection Planes (Γi ) from standard triangulations. The patches
are constructed as a set of triangles joined at a node. . . . . . . . . .

55

Basis functions defined on (a) the projection plane and (b) the true
geometry (surface of a tessellated sphere). The figure (b) shows that
there is no discontinuities introduced by the projection. . . . . . . .

57

Quality of approximation of (a) polynomial function and (b) sinusoidal
functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

4.10 Quality of approximation of discontinuous function: (a) without using
a partition of unity; and (b) with a partition of unity. . . . . . . . .

66

4.11 Comparison of approximation error with and without PU. . . . . . .

70

4.12 Validation studies: Scattering from PEC spheres of different radii. Fig.
(a) shows the RCS comparison and figure (b) shows the convergence
of solutions with basis order p. . . . . . . . . . . . . . . . . . . . . . .

72

4.13 Validation studies: scattering from (a) a PEC plate and an open cavity.
Plot of bistatic RCS (φ = 0). Fig. (b) bistatic scattering RCS from a
NASA thin almond and a PEC conesphere. The RCS is compared at
φ = 0 for the almond and θ = π/2 for the conesphere. Inset figures
show the surface current on the almond and the conesphere. . . . . .

75

4.14 RCS comparison for a PEC box (0.7λ). The figure compares RCS due
to scattering of a −ˆ directed plane wave polarized along −ˆ . The
x
y
bistatic RCS is compared at φ = 0 against that obtained using an
RWG discretization and against measured data. . . . . . . . . . . . .

77

4.15 RCS comparison for regular GMM and mixed order GMM: (a) comparison of bistatic RCS (φ = 0) for a PEC box. Fig. (b) shows the
comparison of bistatic RCS (φ = 0) for a PEC thin almond. In each
ˆ
case, the incident field is directed along −ˆ and polarized along x and
z
bistatic RCS is evaluated at φ = 0. . . . . . . . . . . . . . . . . . . .

80

ˆ
4.16 RCS obtained due to scattering of a x polarized plane wave directed
ˆ
along −ˆ from a PEC box (λ × 0.5λ × 0.5λ) directed along x obtained
z
using regular GMM and mixed basis GMM with polynomials and plane
waves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

82

4.8

4.9

xi

4.17 Backscatter RCS comparison for scattering from a non conformally
meshed PEC plate (shown in figure (a)) RCS obtained by using the
GMM discretization on a non conformal mesh compared against that
obtained by discretizing on a (standard) conformal mesh for a PEC
plate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

84

4.18 EFIE Condition number comparisons for PEC sphere (0.2), square
plate (1m) and cube(1m) using GMM, RWG and loop star basis. . .

85

4.19 EFIE Condition number comparisons for PEC sphere (0.2m) and strip
(1m × 0.2m) using higher order (p = 3) complete GMM and RWG basis. 86
4.20 EFIE Condition number comparisons for non-uniformly discretized
square plate (1m) and thin almond(1 : 3 : 9) using GMM and RWG
and loop star basis. The plate is discretized at average edge length
ranging from 0.1λ to 0.025λ and the almond is discretized from 0.1λ
to 0.02λ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

4.21 Comparison of RCS (GMM vs. loop star) at 30M Hz, 30kHz and 3kHz
for constant h for a Plate: side 1m. The RCS is compared against that
obtained using Loop Star basis functions. . . . . . . . . . . . . . . .

88

4.22 Condition number vs. frequency for a PEC box: The figures show the
condition numbers using the regular GMM basis set and using just the
quasi-divergence free basis functions from the GMM basis set. In each
case; the RCS is computed to good accuracy as shown inset. . . . . .

90

4.23 Bistatic RCS (φ = 0) comparison for a dielectric toroid (outer radius
0.75λ0 and inner radius 0.5λ0 , εr = 4.0). This result demonstrates
the advantage of using an local quasi-surface Helmholtz decomposition
basis function set. The inset shows direction and polarization of the
incident field along with an image of the surface current density on the
toroid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

5.1

RCS obtained using GMM basis functions used to discretize the PMCHWT amd M¨ller equations for the solution of scattering from (a)
u
dielectric spheres (radii 0.1λ0 ,0.3λ0 and 0.5λ0 and εr = 4.0) and (b)
a cube (side λ0 , εr = 4.0) due to plane wave incidence. The bistatic
RCS (at φ = 0) are compared against that obtained using an RWG
discretization and analytical results (for spheres). . . . . . . . . . . . 101
xii

5.2

RCS obtained using GMM basis functions due to scattering from (a)
NASA thin almond of εr = 4.0 that fits in a box of size 0.4λ0 ×1.2λ0 ×
ˆ
3.0λ0 due to a plane wave incident along −ˆ and polarized along x and
z
(b) a cone-sphere of radius 0.5λ0 , cone height 2.6λ0 and εr = 4.0, due
ˆ
to a plane wave polarized along −ˆ and incident along x. Bistatic RCS
y
is computed at φ = 0 for the almond and at θ = π/2 for the conesphere
and used to validate the GMM scheme against RWG. . . . . . . . . . 103

5.3

RCS comparison of standard GMM vs. Mixed-order GMM (a) NASA
almond discretized with GMM basis functions complete to order p = 4
near edges and p = 1 elsewhere. (b) Dielectric box using p = 1 near
edges and p = 2 in larger patches near center. The bistatic RCS is
ˆ
computed at φ = 0 and the incident field is polarized along x and
incident along −ˆ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
z

5.4

RCS comparison of standard GMM vs. Mixed-Basis GMM (5.0λ0 ×
2.0λ0 × 0.5λ0 ). Electric field is incident along −ˆ and polarized along
z
ˆ
x. The bistatic RCS (φ = 0) is compared for GMM using polynomial
basis functions against GMM using polynomials mixed with plane wave
basis functions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

5.5

RCS comparison for scattering from a dielectric arrow (2.0λ0 ×3.0λ0 ×
ˆ
λ0 ). Electric field is incident along −ˆ and polarized along x. The
z
bistatic RCS (φ = 0) is compared for four different combinations of h
and polynomial order of basis function p. . . . . . . . . . . . . . . . 110

5.6

RCS comparison for scattering from a non conformally meshed dielectric box. Bistatic RCS (φ = 0) obtained by discretizing the M¨ller
u
integral equation on a non conformal mesh (shown inset) compared
against that obtained by discretizing on a (standard) conformal mesh
for a thin dielectric (εr = 4.0) box (2.0λ0 × 0.2λ0 × 0.5λ0 ) . . . . . . 112

5.7

Condition number comparisons. All objects are discretized at 0.1λm at
30M Hz and the operating frequency is swept from 30M Hz to 30Hz
and condition number of the impedance matrix obtained on GMM
discretization is compared against that obtained using a standard RWG
discretizaion. (a) Sphere (radius 0.1m, εr = 4.0) and (b) Cube (side
1m and εr = 4.0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

5.9

Bistatic RCS (φ = 0) comparison for a dielectric toroid (outer radius
0.75λ0 and inner radius 0.5λ0 , εr = 4.0). Figure 5.9(a) shows direction
and polarization of the incident field along with an image of the surface
current density on the toroid. . . . . . . . . . . . . . . . . . . . . . . 117
xiii

6.1

General description of a scattering problem . . . . . . . . . . . . . . . 120

6.2

Temporal basis functions for solution of the TDIE. (a) shows the Lagrange polynomial basis functions which are not band limited as opposed to (b) approximate prolate spheroidal functions that are band
limited but require extrapolation to maintain causality. . . . . . . . . 126

6.3

Backscatter RCS comparison for (a) two spheres of radii 0.2m and
0.5m and (b) a thin NASA almond (0.3m × 1.15m × 3m) compared
against RWG discretization. The inset figure shows snapshots of the
magnitude of the current density on the almond. . . . . . . . . . . . 131

6.4

Backscatter RCS comparison for (a) a thin box of size 20m×1m×20m
compared against RWG discretization and (b) a PEC cube of side
0.7m compared against both RWG and experimental result. The inset
figures show snapshots of the magnitude of the current density on the
scatterer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

6.5

Backscatter RCS comparison for (a) thin box (20m × 1m × 20m) and
(b) a toroid (outer radius 0.75m and inner radius 0.5m). The inset
figures show a snapshot of the magnitude of the current density and
direction of propagation and incidence. . . . . . . . . . . . . . . . . 135

6.6

Backscatter RCS comparison of standard GMM vs. Mixed-order GMM
(a) Box using p = 1 near edges and p = 2 in larger patches near center
(b) NASA almond discretized with GMM basis functions complete to
order p = 2 near edges and p = 1 elsewhere and (c) dielectric box using
p = 1 near edges, p = 3 near the center and p = 2 in the transition
region. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

6.7

Backscatter RCS comparison for scattering from a non conformally
meshed PEC box. RCS obtained by using the GMM discretization on
a non conformal mesh (shown inset) compared against that obtained
by discretizing on a (standard) conformal mesh for a thin box . . . . 142

xiv

Chapter 1
Introduction
The preface to a collection of seminal papers on computational electromagnetics,
compiled in 1973, says “Computer electromagnetics is a relatively new and rapidly
growing field” [1]. Today, almost four decades later, computational electromagnetics
is a field that is still growing equally rapidly. The growth of computational electromagnetics has provided the community with simulation tools that solve a wide array
of problems ranging from electrical to electronic to computer to optical. For most of
the initial part of this growth, computational tools were largely used to complement,
and not supplant, experiments. The design of most electromagnetics structures such
as antennas, gratings etc. were reliant on the designers’ intuitive understanding of
the underlying physics. This, in turn, implied that non-intuitive designs were automatically precluded. However in the last two decades, computational techniques have
grown by leaps and bounds in terms of speed and accuracy, and the simulation of
complicated, realistic structures is well within reach. Today, therefore, computational
techniques are a viable alternative to experiments in a wide variety of applications.
The idea of using computational techniques for the solution of partial differential equations historically comes from the field of fluid dynamics [2–5] where the most
problems, even for reasonably “canonical” geometries, are analytically intractable ow1

ing to the inherent non linearity of the Navier Stokes equations. Maxwell’s equations,
by contrast, are linear and therefore, permit a wide class of analytical and asymptotic
approaches. However, the solution of realistic electromagntic problems on complex
geometries and the rapid advancements in computing power strongly motivates the
use of computational techniques. The spate of techniques that exist today in computational electromagnetics generally fall into two categories: differential equation
based and integral equation based. In the next two sections, we will visit the salient
features of both these methods.

1.1

Differential Equation Methods

Differential equation based methods [6–8] have the advantage of a simpler formulation
since they directly discretize the PDE. As a result, they lend themselves easily to theoretical analysis viz-a-viz existence of solutions, stability and convergence and there is
a plethora of work that has been devoted to these studies [7, 9–17]. Differential equation based techniques also yield matrix systems that have O(Ns ) entries, where Ns is
the number of degrees of freedom. The direct discretization of the PDE also implies,
however, that boundary conditions need to be explicitly enforced in the numerical
scheme. Extensive work has been done on the incorporation of boundary conditions
including the radiation boundary condition [18] which has seen considerable mathematical work. Differential equation methods have seen an explosion in application
specific research with the development of hierarchical basis functions [8, 19], discontinuous Galerkin methods [20] and vector generalized finite elements [21], and domain
decomposition methods [22] among others.
2

1.2

Integral Equation Methods for Electromagnetics

Integral equation methods [23–30] offer more advantages than differential equation
based methods for scattering problems primarily because the number of unknowns is
reduced significantly, for surface scattering problems. This does come at the expense
of a theoretically more involved formulation due to which these methods have eluded
complete theoretical analysis for a very long time. While several recent efforts have
been made to fully understand the theoretical nature of integral equations [31–37],
several questions about the existence, uniqueness and convergence of solutions remain
open. The second problem is the full matrix that needs to be solved as a result of
the numerical discretization of the integral equation, which yields a matrix system of
2
size 0(Ns ); correspondingly leading to higher demands on computational complexity
3
O(Ns ), and memory.
Recently however, the state of art of integral equation (IE) based methods has
received a shot in the arm with the development of fast methods [38–40] which reduced
2
the computational complexity and memory scaling from O(Ns ) to O(Ns log α (Ns )),
where 0 ≤ α ≤ 2. Steady advances in this technology (a woefully incomplete list being
Refs. [38,39,41–47]) has made the analysis of radiation and scattering from electrically
large structures possible; see [48] for a more comprehensive list of applications and
recent advances. In a similar manner, the development of the plane wave time domain
method (a time domain analogue of the fast multiple method) has enabled transient
analysis of electrically large and complex geometries [49–58]. This method reduces
2
the computational complexity from O(Nt Ns ) to O(Nt Ns logα Ns ), where Nt is the
number of temporal degrees of freedom.
While electromagnetic problems have been analysed using integral equations since
1897 when Poklington analysed wave propagation along wires, the most commonly
3

used numerical technique for the “discretization” of integral equations is the method
of moments (MoM), proposed for electromagnetics by Harrington [24]. This technique
starts out by expanding the unknown quantity in a set of basis functions and testing
the resulting integral equation using a space of testing functions. While it is possible
to define “whole domain” basis functions for some canonical problems; most realistic
problems are solved using a simplicial tessellation of the geometry. The basis functions
are then defined local to each patch of this tessellation. It can be shown that the
acuraccy of the integral equation scheme relies heavily on definition of the basis
functions. As a result, there has been a significant amount of work that has gone into
the design of appropriate basis functions. In the rest of this chapter, we will focus on
the design of these basis functions as this is the main theme of the dissertation.

1.3

Basis Functions for Electromagnetic Integral
Equations

Historically, the design of basis functions, and indeed the solution of the integral
equation itself has been motivated by the design of the tessellation for the geometry. The first method utilized for the solution to scattering problems using integral
equations was to model the surface with a wire mesh, called a wire-grid model. This
method with its ease of implementation and ability to solve a rather large variety
of scattering problems, remained popular in the community for almost two decades
starting with [59]. However the method is inherently inaccurate in the model for the
computation of near field quantities and input impedance. Furthermore, the method
also suffered from ill conditioned matrices, numerical resonances [60] and other theoretical issues [61]. Surface patch modelling, developed as an alternative to wire-grid
modelling was seen to preform better on all these grounds and therefore, much of
the research gravitated towards developing basis functions for surface patch models.
4

The first set of surface patch models started with both planar [62] and non-planar
quadrilaterals [63] for the magnetic field integral equations (MFIE). This was followed
by planar rectangular patches for the electric field equations by [64]. Sinusoidal basis functions were then extended to surface patches from their wire-grid equivalents
for the electric field integral equation (EFIE) [65]. This was followed by what was
probably the first attempt to model surfaces using triangular patches [66]. While the
approach in [66] modeled square plates, the authors astutely pointed out the ability
of triangular patches to model arbitrary surfaces and their work was also one of the
first to model the entire integral equation using an equivalent of the Galerkin scheme
in use today. Planar triangular patches along with the MFIE and phase variation into
the definition of the basis function for the MFIE was first introduced by [67, 68]. In
what was the first work that directly implemented surface patch models for arbitrary
surfaces, the authors used non planar triangular patches to model the MFIE [69]. The
subsequent development of a special set of basis functions for triangles by Rao and
Glisson for arbitrary triangulations, analogous to roof-top functions for rectangular
patches [70] was a critical contribution to the research. The authors in [70] followed up
this work with what is arguably the most influential paper in the method of moments
research where they developed the Rao-Wilton-Glisson (RWG) functions to model
vector surface currents on non-planar triangulations of arbitrary surfaces [71]. These
basis functions also support an important property called divergence conformance.
This means that the divergence of the basis functions do not have any discontinuities
across internal patch boundaries. It must be mentioned that a curl conforming equivalent of these basis functions were also independently introduced by Nedelec [72] et
al. and implemented by Bossavit et al. [11], Barton et al. [14] and others and the
divergence conforming functions were also independently introduced by Raviart and
Thomas [73].

5

1.4

Higher order basis functions

While we will revisit the RWG functions in great detail elsewhere in this thesis, suffice
it to say that it has become the de-facto standard for the method of moment analysis of electromagnetic integral equations. Several higher orders and singular variants
have been proposed in literature [74,75] and the RWG method has become the corner
stone for all frequency and time domain integral equation analysis over the past three
decades. However, while higher order basis functions have been in use for a few years,
the ability to include arbitrary functions in the basis space or to switch from one basis
function to another is still limited. This is a direct consequence of the historically
close relationship between the basis function and the underlying tessellation. Both
the original RWG and its higher order/singular variants have been designed to satisfy continuity properties across internal boundaries of the tessellation. It is widely
accepted that any basis function designed to model continuous quantities, must satisfy some constraints on continuity across patches. As a result, the incorporation of
different classes of functions in the basis space is only possible within a framework
that ensures continuity [76–78].
The ability to include arbitrary functions may be critical when some information is
known about the physics of the problem, and it is desirable to include this behavior in
the definition of the basis functions. For example, [79–81] incorporate Meixner series
based functions in the basis space to model singular currents near edges. However,
if the marriage with the tessellation could be removed, it would be possible to mix
arbitrary classes of basis functions within the same geometric region. Unfortunately,
all existing methods are specific to the underlying tessellation. Indeed, inclusion of
a completely arbitrary basis function space might not even be possible within the
framework of current basis function definitions.
Along these lines, a question has recently been posed whether the explicit enforce6

ment of continuity should be required [82, 83]. Indeed [82] makes arguments that, in
two dimensions, the use of a moderately high order complete basis function set eliminates the need for explicit enforcement of continuity conditions. Likewise, it has been
demonstrated in two dimensions, that the continuity constraints can be replaced by
a partition of unity function to maintain the h and p convergence [84].
The problems of enlarging functions spaces has been addressed in the finite element
community. Melenk et al. [76] developed the scalar GFEM as an umbrella framework
for a class of mesh-based and mesh-less methods that allowed for the inclusion of a
larger class of basis functions. He demonstrated that the global errors were bounded
by the local errors in approximation. The scheme was extended to vector finite
elements for electromagnetics by Lu et al. [77]. The GFEM uses a product of a
partition of unity function and a higher order approximation function to define basis
functions on overlapping domains. However, extension of this method to integral
equations is not trivial and will form the central focus of this thesis. In GFEM,
overlapping patches are defined by combining canonical patches. However, unlike in
the GFEM framework, the patches for surface integral equations need to conform to an
arbitrary scatterer surface and can not be expected to map on to canonical geometries.
This, in turn, makes the definition of vector basis functions on these non-canonical
surfaces challenging. The structure of these patches may lead to functions that are
piece-wise differentiable (for piecewise continuous patches). Therefore, computation
of matrix elements involves evaluating singular and hyper singular integrals that need
to be carefully handled. Furthermore, one needs to define vector basis functions on
these surfaces to represent currents. This poses additional challenges.
In this dissertation, we will describe the design and implementation of a new
umbrella framework for the discretization of integral equations called the Generalized
Method of Moments. We will show that within this framework it is possible to choose
from a very wide variety of functions to form the basis space. While the choice of
7

basis functions can be completely arbitrary, the thesis will describe in detail one
particular scheme that utilizes a surface quasi-Helmholtz decomposition and provides
a mechanism for the use of any available knowledge of the physics of the problem. The
thesis will describe the theoretical framework required for the design of this scheme
and develop error bounds for the approximation of the unknown quantities using the
proposed basis functions. Following that, we will delve into the implementation of the
GMM starting from commonly available triangular tessellations, carefully identifying
and solving several key issues along the way.
The GMM will be applied to a plethora of examples, validated against analytical
solutions and compared against the standard RWG discretizations. It will further
be shown that the use of the proposed basis functions will result in a scheme that
is well conditioned over a wide range of frequencies. We will also provide a host
of examples showcasing the advantages of the GMM scheme, including its ability to
incorporate multiple types and orders of basis functions and its ability to handle non
conformal meshes. We will then extend the GMM to both the analysis of scattering
from dielectric bodies using the PMCHWT and M¨ller operators and to the analysis
u
of transient scattering from conducting objects. Several results will be presented in
both cases that demonstrate the advantage of the proposed method over traditional
RWG based discretization.
The remainder of this thesis will be organized as follows. Chapter 2 will lay out
the mathematical preliminaries required for the rest of this thesis. We will start from
the Maxwell’s equations and derive an integral formulation for the solution of electromagnetic scattering from surfaces. We will specialize the derivation to conducting
structures in preparation for Chapter 3 where we describe the method of moments
technique for the discretization of integral equations and describe the Rao Wilton
Glisson basis functions that are most commonly in use today. In Chapter 4 we will
describe the GMM in detail. This chapter will first describe the theoretical founda8

tions of this work and then proceed to discuss its implementation details. It will also
present a variety of results to validate the technique and demonstrate its advantages.
The following chapter will discuss the extension of the GMM to the analysis of scattering from dielectric bodies. The formulation will be discussed followed by specifics
of the implementation and results. Chapter 6 will describe the GMM applied to the
analysis of transient scattering. This chapter will start with the definition of the time
domain integral equations and proceed to describe the standard temporal discretization techniques. This will be followed by the application details of the GMM to time
domain integral equations and a gamut of results describing the advantages of the
new technique. Finally Chapter 7 will provide a summary and describe possibilities
for future work.

9

Chapter 2
Mathematical Preliminaries
In this chapter, we will present the mathematical preliminaries from which the rest of
this thesis will follow. The equations and derivations presented herein can be found
in a number of references including [29, 85, 86].

2.1

Maxwell’s Equations

To lay the foundation for the integral equations for electromagnetics, we start from
the Maxwell’s equations in their point form. Consider a linear, isotropic, homogeneous region in space having relative permittivity εr and relative permeability µr .
Electromagnetic fields E(r, t), H(r, t) and corresponding flux densities D(r, t) and
B(r, t) in this region satisfy Maxwell’s equations:
∂B(r, t)
− M(r, t)
∂t

(2.1a)

∂D(r, t)
+ J (r, t)
∂t

(2.1b)

× E(r, t) = −

× H(r, t) =

· D(r, t) = (r, t)

(2.1c)

· B(r, t) = m (r, t),

(2.1d)

10

where J (r, t) and (r, t) are the electric current and charge density and M(r, t) and
m (r, t) are their (fictitious) magnetic equivalents. Further, the current and charge
satisfy the continuity equations given by
∂ (r, t)
=0
∂t

(2.2a)

∂ m (r, t)
=0
∂t

(2.2b)

D(r, t) = ε0 εr (r, t) ∗t E(r, t)

(2.3a)

B(r, t) = µ0 µr (r, t) ∗t H(r, t),

(2.3b)

· J (r, t) +

· M(r, t) +
In addition, the constitutive relations

and

∞
where ∗t indicates a convolution in time defined as a(t) ∗t b(t) =

dτ a(t − τ )b(τ ),
0
relate the flux densities to the corresponding fields. The assumptions of linearity
allow for the application of the principle of superposition and thereby, allows for the
application of the Fourier transform. Equivalently, if we assume that all the time
.
dependent quantities are time harmonic of the form A(r, t) = Re A(r)ejωt , where
ω is the frequency in radians per second (and A represents any of the vectors above),
we obtain the familiar time harmonic form of equations (2.1a)-(2.1d) given by

× E(r) = −jωµ0 µr (r)H − M(r)

(2.4a)

× H(r) = jωε0 εr (r)E(r) + J(r)

(2.4b)

ε0
µ0

· εr (r)E(r) = ρ(r)

(2.4c)

· µr (r)H(r) = ρm (r),

(2.4d)

11

and the equations of continuity transform to

· J(r) + jωρ(r) = 0

(2.5a)

· M(r) + jωρm (r) = 0.

(2.5b)

and

For simplicity, we have used the same symbols for εr and µr in frequency and time
domain. At the interface between two media 1 and 2 with material parameters µ1(2)
and ε1(2) , the fields satisfy the boundary conditions given by
ˆ
n(r) × (E1 (r) − E2 (r)) = −Ms

(2.6a)

ˆ
n(r) × (H1 (r) − H2 (r)) = Js

(2.6b)

ˆ
n(r) · (ε1 (r)E1 (r) − ε2 (r)E2 (r)) = −ρs

(2.6c)

ˆ
n(r) · (µ1 (r)H1 (r) − µ2 (r)H2 (r)) = −ρm,s ,

(2.6d)

and

ˆ
where n is the normal to the material boundary from medium 2 to 1.
To derive the vector wave equations for electric field, we take a curl of (2.4a) and
substitute from (2.4b) to obtain

×

1
µ(r)

× E(r) − ω 2 ε(r)E(r) = −jωJ(r) −

×

1
M(r)
µ(r)

(2.7)

.
.
where µ(r) = µ0 µr (r) and ε(r) = ε0 εr (r). A similar equation can be derived for the
magnetic field:

×

1
ε(r)

× H(r) − ω 2 µ(r)H(r) = −jωM(r) +
12

×

1
J(r) .
ε(r)

(2.8)

An important simplification is obtained upon assuming that the material is homogeneous, which results in

×

× E(r) − k 2 (r)E(r) = −jωµJ(r) −

× M(r)

(2.9)

× H(r) − k 2 H(r) = −jωεM(r) +

× J(r).

(2.10)

.
where k = ω 2 εµ, and

×

2.2

Vector and Scalar Potentials

While it is possible to solve the above equations to obtain the fields for a given source
distribution, it is often easier to proceed by defining a set of auxiliary potentials.
To this end, first consider Maxwell’s equations in a medium free of magnetic sources
(ρm (r) = 0 and M(r) = 0). Since

· B(r) = 0 under this assumption, we can define

the magnetic vector potential, A(r) such that B(r) =

× A(r). Combining this with

(2.4a), we obtain the following equations

× A(r)

B(r) =

× (E(r) + jωA(r)) = 0.

(2.11)
(2.12)

Equation (2.12) implies that E(r + jωA(r can be written as a gradient of a scalar
Φe , called the electric scalar potential, as

E(r) + jωA(r) = − Φe (r).
13

(2.13)

To obtain a unique solution for the field from these potentials, we need an additional
constraint. One such constraint is the Lorentz gague, given by

· A + jωµεΦe = 0.

(2.14)

Substituting the definition of A into (2.4b) and assuming, again, a homogeneous
medium, we obtain
1
µ

×

× A(r)

= jωε[−jωA(r) −

Φe (r)] + J(r)

(2.15a)

1
= jωε[−jωA(r) + jωµε
· A(r)] + J(r)
1
· A(r) + J(r)
= ω 2 εA(r) + µ

(2.15b)

2 A(r) + k 2 A(r) = −µJ(r).

(2.15d)

(2.15c)

leading to

Substitution of the Lorentz gauge directly into (2.13) will result in an equation for
Φe as
2 Φ (r) + k 2 Φ (r) = − ρ(r) .
e
e
ε

(2.16)

Either by inspection, or by completing a similar derivation for a medium free of electric
sources (but containing magnetic sources) result in the definition of the electric vector
potential F(r) and the magnetic scalar potential Φm (r) and the following related
equations.

D(r) = −

× F(r)

H(r) = −jωεF(r) −

(2.17)
Φm (r)

(2.18)

2 F(r) + k 2 F(r) = −εM(r)

(2.19)

2 Φ (r) + k 2 Φ (r) = − ρm (r)
m
m
µ

(2.20)

14

By superposition, in a medium having both magnetic and electric sources, the field
equations in terms of the potentials can be written as
1
ε

× F(r)

(2.21)

× A(r) − jωF(r) −

Φm (r)

(2.22)

E(r) = −jωA(r) −
H(r) =

1
µ

Φe (r) −

In Cartesian coordinates, the splitting of the vector Laplacian into its components
motivates the analysis of (2.15d) and (2.16) (or equivalently, (2.19) and (2.20)) using
scalar wave equations of the form
( 2 + k 2 )ψ(r) = f (r),

(2.23)

where ψ denotes each component of the field and f denotes the corresponding source
component. As with any partial differential equation, obtaining a unique solution to
the above equation requires boundary conditions.

For scattering problems, which will form the bulk of the problems that we deal
with in this thesis, we solve the above equation subject to the Sommerfield radiation
boundary condition. For scalar waves this condition is given by

lim |r|
|r|→∞

∂ψ(r)
+ jkψ(r)
∂|r|

= 0.

(2.24)

In order to solve the above partial differential equation, one can define a Green’s
function for the differential equation and its boundary condition, given by

g(|r|) =

e−jk|r|
.
4π|r|

15

(2.25)

The solution to (2.23) is given by

ψ(r) = g(|r|) ∗s f (|r |),

(2.26)

.
where ∗s denotes the spatial convolution defined by a(r)∗s b(r) =

dr a(r−r )b(r ).

Ω
Using the Green’s function, the electric and magnetic potentials can be now written
as

A(r) = µg(|r|) ∗s J(r),

(2.27)

F(r) = εg(|r|) ∗s M(r),

(2.28)

1
Φe (r) = g(|r|) ∗s ρe ,
ε
1
Φm (r) = g(|r|) ∗s ρm .
µ

2.3

(2.29)
(2.30)

Integral Equation Representation of Solutions

The above preamble sets the stage for the derivation of the electric and magnetic field
integral equations for the fields. We follow the derivation of [29], which in turn is
derived from [87] and [88] but extended to non-smooth surfaces. We start with the
vector Green’s identity for any two vector functions P(r) and Q(r) with continuous
first and second derivatives.

dr (Q(r) ·

×

× P(r) − P(r) ·

×

× Q(r))

V
ˆ
dr n(r) · (P(r) ×

=

× Q(r) − Q(r) ×

× P(r)) .

(2.31)

∂V
Assume that the volume V is represented by D+ and its inner boundary by Ω and
the outer boundary by Ω0 . Let the volume bounded by the inner boundary be D− .

16

ˆ
ns

D+
Ω

ˆ
ns
D−

Ω0

Figure 2.1: General representation for a domain D+ = V with inner and outer
boundaries

ˆ
Since Q(r) is arbitrary, we choose an arbitrary pilot vector a and define
e−jkR
.
ˆ
ˆ
Q(r) = ag(R) = a
,
4πR

(2.32)

.
where R = |r − r | and r and r are the observation and source locations respectively.
.
Further, choice of P(r) = E(r) results, after some detailed vector manipulation [87]
17

in

D+
dr

dr

jωµJ(r )g(R) + M(r ) ×

g(R) −

ρ(r )
g
ε

=

(2.33)

ˆ
ˆ
jωµ n(r ) × H(r ) g(R) − (ˆ (r ) × E(r ))g(R) − n(r ) · E(r )
n

g(R)

Ω+Ω0
(2.34)
For simplicity we will assume at this point that the observation location r does not
lie on Ω0 and that Ω is smooth. The case where these assumptions fail happen can
be dealt with explicitly and the reader is referred to [29] for an excellent description.
However, under this assumption, (2.34) reduces to

T E(r) =

D+

dr

−−

dr

jωµJ(r )g(R) + M(r ) ×

g(R) −

ρ
ε

g(R)

ˆ
jωµ n (r ) × H(r ) g(R)

Ω+Ω0
−ˆ (r ) × E(r ) ×
n
ˆ
− n (r ) · E(r )

g(R)
g(R) ,

(2.35)

where primed coordinates relate to the quantities (and operations) defined on the
source position r and − indicates that the integral is performed by excluding the
principal value contribution from the singularity r = r . T is a parameter such that
.
.
T = 2 if r ∈ Ω and T = 1 otherwise. By observation (and applying the duality of

18

Maxwell’s equations) a similar equation for H(r) can be written as

T H(r) =

D+

−jωεM(r )g(R) + J(r ) ×

dr

+−

g(R) +

ρm
µ

g(R)

jωε[ˆ (r ) × E(r )]g(R)
n

dr
Ω+Ω0

+ˆ (r ) × H(r ) ×
n
ˆ
+ n (r ) · H(r )

g(R)
g(R) .

(2.36)

Further, allowing Ω0 to move out to ∞, it can be seen that the integral over Ω0
vanishes as a result of the radiation boundary conditions resulting in
T E(r) = Ei (r)+

dr

jωµJ(r )g(R) + M(r ) ×

− − dr
Ω

−ˆ (r ) × E(r ) ×
n
−[ˆ (r ) · E(r )]
n

g(R) −

ρ
ε

jωµ[ˆ (r ) × H(r )]g(R)
n

D+

g(R)

g(R)
g(R) ,

(2.37a)

and
T H(r) = Hi (r)+

D+

dr

+ − dr
Ω

−jωεM(r )g(R) + J(r ) ×

g(R) +

ρm
µ

g(R)

jωε[ˆ (r ) × E(r )]g(R)
n
+ˆ (r ) × H(r ) ×
n
+[ˆ (r ) · H(r )]
n

g(R)
g(R) ,

(2.37b)

where we have denoted all the fields outside ∂V0 by Ei and Hi and will call them
the incident fields resulting from sources at or beyond infinity.
These equations are called, respectively, the electric and magnetic field integral
equations (EFIE,MFIE). We will now specialize these equations to the first case that
19

will be studied in the next chapter. For the rest of the thesis beyond chapter 4, we
will revisit this formulation and specialize it to other cases as necessary.

2.4

EFIE and MFIE for PEC structures

Consider a perfect electrically conducting body occupying a domain D− ⊂ R3 in free
space. Let Ω denote the bounding surface of D− and let the outward pointing unit
ˆ
vector normal at any point r on Ω be denoted by n(r). Assume that a time harmonic
field {Ei (r), Hi (r)} is incident on this body. This field induces currents on the surface
Ω, that radiate scattered fields {Es (r), Hs (r)}. Thus, the total electric and magnetic
fields in all of D+ = R3 /D− , denoted by E(r) and H(r), can be decomposed into
the incident and scattered fields by
E(r) : R3 → C3 = Ei (r) + Es (r)

(2.38a)

H(r) : R3 → C3 = Hi (r) + Hs (r)

(2.38b)

and

{Ei (r), Hi (r)}
ˆ
ns
D−

{Es (r), Hs (r)}
Ω

Figure 2.2: General description of a scattering problem
20

Since the domain D− is a perfect electric conductor, further simplifications can
be effected to the equations (2.37a) and (2.37b) above. We can specialize equations
(2.6a) and (2.6d) and the equation of continuity to the case of H− = E− = 0 and
H+ = H , E+ = E where the subscript − and + denotes quantities in D− and D+ ,
respectively. We then, obtain for any r on the surface Ω,

ˆ
n(r) × E(r) = 0,

(2.39a)

ˆ
n(r) · H(r) = 0,

(2.39b)

ˆ
n(r) × H(r) = J(r),

(2.39c)

ˆ
n(r) · E(r) =

j
ωε

· J(r),

(2.39d)

where J is the surface current on Ω. Substituting into (2.37a) and (2.37b) and applying the boundary conditions on the tangential components of the fields, on the
surface (Ω) to arrive at
ˆ
n(r) × (Ei (r) + Es (r)) = 0

(2.40a)
r∈Ω

and
ˆ
n(r) × (Hi (r) + Hs (r)) = 0

,

(2.40b)

r∈Ω
gives us
ˆ
ˆ
n(r) × Ei (r) = n(r) ×

dr

jωµJ(r )g(R) + [

Ω

j
ωε

· J(r )]

g(R) , (2.41a)

and
ˆ
J(r) = T n(r) × Hi (r) + T − dr J(r ) ×
Ω

g(R).

Further simplification using vector identities and the observations that
21

(2.41b)
g(R) =

ˆ
− g(R) and n · J = 0, along with setting T = 2 (since r ∈ Ω) will lead to
ˆ
ˆ
n(r) × Ei (r)) = n(r) ×

dr

jωµg(R)J(r ) +

Ω

j
ωε

g(R) · J(r ) ,

(2.42a)

and
ˆ
n(r) × Hi (r) =

J(r)
ˆ
− n(r) × − dr J(r ) ×
2
Ω

which are the surface EFIE and MFIE for PEC bodies.

22

g(R).

(2.42b)

Chapter 3
Numerical Solution of Integral
Equations
The integral equations that were developed in chapter 2 are the so-called “surface”
integral equations since the field quantities in all space are expressed in terms of their
values (or the values of some equivalent quantities) on the surface of the scatterer.
These integral equations cannot be analytically solved except for canonical scatterer
geometries. This motivates the development of numerical solution techniques to these
equations.

3.1

Method of moments solution for integral equations

The method of moments (MoM), developed for electromagnetic problems in the seminal work by Harrington [24] proposes an elegant technique for the solutions to integral
equations. Consider an equation of the form

E(r) = L ◦ {J(r)},
23

(3.1)

on a domain Ω, where L represents either of the integral equations in (2.42a) or
(2.42b). The MoM solution to this equation is sought in two steps. The first is to
obtain a finite dimension approximation to the unknown J(r) using a set of basis
functions ji (r), namely
.
J(r) ≈ JN (r) =

N
ai ji (r).

(3.2)

i=1
Using this, and recognizing the linearity of the integral operator, (3.1) can now be
re-written as

N
ai L ◦ {ji (r)}.

E(r) =

(3.3)

i=1
To complete the process of discretization of this equation, the equation is enforced
using a set of testing functions {tj (r)}N to obtain
j=1

Ω

dr E(r)tj (r) =

Ω

ai L ◦ {ji (r)}.

dr tj (r)

(3.4)

i

This procedure transforms the integral equation into an equivalent matrix equation
of the form
e=Z a

(3.5a)

where the elements of the RHS vector e and “impedance matrix” Z are given by

ej =
Z j,i =

Ω

Ω

dr E(r)tj (r)

dr tj (r)L ◦ {ji (r)}
(3.5b)

and a = {a1 , a2 , · · · } is the vector of basis function coefficients. Significant work has
gone into the proof of the validity of the method of moments, especially in the special
case where the function spaces that {ji } and {tj } are identical (referred to as the
24

Ritz-Galerkin procedure) [12, 13]. While the issue of the choice of functions spaces
for {ji } and {tj } are important, and indeed critical; we will defer that discussion to
a later stage. At this point, we will assume that Ω is smooth enough and therefore,
that {ji } and {tj } can be picked from the standard L2 spaces.

3.2

Approximation errors

The method of moments procedure described above introduces errors in the solution
of the integral equation. The primary source of error is the approximation of the
function. This approximation error, εa is defined as follows.
Definition 1.
.
ε a = J − JN = J −

N
ai ji (r)

(3.6)

i=1
Theorem 3.2.1. For any linear and bounded integral operator L, ∃ a constant C
such that

Ω

dr tj (r)L ◦ {JN (r)} −

Ω

dr tj (r)L ◦ {J(r)} ≤ Cεa

(3.7)

Proof. By linearity of the integral operator, we have

Ω

dr tj (r)L ◦ {JN (r)} −

Ω

dr tj (r)L ◦ {J(r)} =

Ω

dr tj (r)L ◦ {JN (r) − J(r)}
(3.8)

Using Cauchy Schwartz inequality we have:

Ω

dr tj (r)L ◦ {JN (r) − J(r)} = tj (r) 2 L ◦ {JN (r) − J(r)} 2 .

(3.9a)

Since tj (r) ∈ L2 (Ω), we have
tj (r) 2 ≤ C2 ,
25

(3.9b)

and finally by the boundedness of the integral operator, we have

L ◦ {JN (r) − J(r)} 2 ≤ C1 JN (r) − J(r) 2

(3.9c)

Combining the two results, we have,

Ω

dr tj (r)L ◦ {JN (r) − J(r)} ≤ C

JN (r) − J(r)

2 = Cεa .

(3.9d)

The above result motivates the need for the accurate representation of the currents
in terms of the basis functions. In the MoM community, there have been a wide variety
of basis functions that were proposed for the solution of integral equations [59, 62,
63, 66, 69–71, 74]. While global basis functions can provide excellent approximation
qualities on canonical geometries; for most practical problems, the best that can
be hoped for, is a good local approximation to the currents. To this end, all basis
function schemes typically start out by subdividing the geometry using a simplicial
tessellation.
The ease and generality of triangular tessellations led to a basis function scheme
called the Rao-Wilton-Glisson (RWG) basis functions [71] that have become the defacto standard for the MoM solution to integral equations in electromagnetics. The
number of applications that use these functions today are innumerable. These functions have since been extended to the analysis of dielectric objects [89] , layered
media [90,91], transient scattering [92] and many other applications. Though the initial work was specific to first order basis functions, since then both higher order [74]
and singular [81] variants that have been developed.
The RWG basis functions are designed on each interior edge of the underlying
mesh. Fig. 3.1 shows two triangles sharing a common edge. On this edge, a basis
26

function is defined such that it vanishes everywhere on the surface except within the
two triangles attached to the edge. Points inside each triangle is designated by a
±
position vector ρ± associated with the free vertex Tn . The vector basis function on
these triangles can be defined as

 ln +

+

 + ρn r ∈ Tn
 2A
 n


ln −
−
fn (r) =
 2A− ρn r ∈ Tn
 n





0 otherwise

(3.10)

±
where ln is the length of the edge and A± is the area of the triangle Tn . The basis
n
function fn satisfies the following useful properties.
1. The normal component of the basis function vanishes on all boundaries except
the common boundary
2. The normal component of the basis function on the common edge is continuous
and, for each interior basis function, cancels the corresponding contribution
from the second triangle supported by the edge.
3. The surface divergence of the basis function is given by

 ln

+

 + r ∈ Tn
 2A

 n

−
·s ·fn (r) = − ln
 2A− r ∈ Tn


n




0 otherwise

(3.11)

The RWG basis functions have become the cornerstone for the method of moment
analysis over the past three decades. However, their intimate relationship with the
underlying tessellation, while providing some unique advantages, also has some im27

T−

T+
ln

¯
ρ−

¯
ρ+

2A−
n
ln

2A+
n
ln
Figure 3.1: A schematic showing the definition of the RWG basis functions on an
edge shared by two triangles

28

portant drawbacks, the amelioration of which will be the central focus of this thesis.

29

Chapter 4
The Generalized Method of
Moments
The standard RWG-MoM scheme as described in chapter 3 utilizes a set of basis
functions defined on a tessellation of the geometry. The RWG functions, as described
in [71] were designed for a triangular tessellation of the geometry and the function
spaces offer several advantages that are very useful in their implementation; (i) they
are divergence conforming, (ii) the normal component vanishes on the interior edges
of the geometry and (iii) the normal component is continuous across the edge. These
basis functions have become the cornerstone of most MoM methods for electromagnetic integral equations. However, notwithstanding the significant advantages of both
triangular tessellations themselves and RWG functions, the intimate relationship between the basis function and tessellation comes with some significant disadvantages.
The chief problem with having a basis function specification closely tied to the
discretization of the geometry is that it severely limits the size of the function space
from which basis functions can be chosen. It is worth mention that more than two
decades after the introduction of the RWG functions there are very few higher order
variants in active use today [19]. In fact all higher order and singular variations of
30

the RWG [74, 81] are constructed by multiplying suitably chosen polynomials with
first order RWG functions. This is a direct consequence of the construction of the
RWG functions and implies that the close relationship between tessellation and basis
functions is carried over to higher orders as well. Not only is the inclusion of arbitrary
functions in the basis space is almost completely precluded by the construction of a
basis function so tailored to the tessellation; even the use of mixed orders of basis
functions requires careful consideration to maintain order continuity across interior
boundaries [81]. This is especially critical in situations where some knowledge is
available about the local behavior of the current.
Consider for illustration purposes, scattering from a PEC cube. Fig. 4.1 shows
the distribution of currents when illuminated by an electric field, polarized along x
ˆ
and incident along −ˆ. What is clear from the figure is that the currents are smooth
z
on each face of the cube and the variations are limited to the edges. This would
suggest the use of lower order basis functions or larger size patches on the faces of
the cube which would serve to both reduce computational complexity and increase
accuracy. However, this procedure of mixing basis functions is extremely non-trivial
in the case of RWG basis functions.
While the simple case of a cube serves for illustration, this need for mixing basis
functions becomes magnified in the case of realistic, multi-scale, problems. In existing
basis function frameworks for MoM, it is difficult to mix different orders and impossible to mix different types of functions. Both these features deal with enlarging the
possible functions that can be used for approximating the currents; a feature that
would be beneficial in reducing the computational cost and increasing accuracy. For
instance, in a geometrically complex object that is electrically large, the current behaves differently in different regions. Near regions of geometric complexity, it may be
better to use either high density meshes with low order polynomials or vice-versa in
areas of low geometric complexity. Near an edge singularity, one needs to capture the
31

Figure 4.1: Currents on a PEC cube showing the singular nature of the currents
near the edges and the smooth currents towards the center. For interpretation of the
references to color in this and all other figures, the reader is referred to the electronic
version of this dissertation.

singular behavior of the fields. In smooth regions, the phase of the current may be
better captured using plane-wave type functions. The ability to mix approximation
spaces become important when it is necessary to capture fields with high fidelity.
Next, given the complexity of creating a tessellation for complex objects, this task
is typically partitioned to multiple sub-domains. More often than not, this results
in meshes that are non-conformal. The ability to include these class of meshes with
minimal effort is valuable to the accurate analysis of scattering from such geometries
and is another feature missing in the state-of-art of the basis functions today.
Finally, in many realistic geometries, often the smallest geometrical feature dictates the tessellation and thereby the basis function space that is used to discretize
the equation. This in turn leads to the “low-frequency problem” where the resulting
matrix equation is very badly conditioned. There have been several attempts to solve
this problem in literature where either special basis functions, called loop and star
functions, are constructed specifically for low frequency problems using linear combinations of RWG functions [93, 94], or a different, well conditioned, integral equation
is derived using mathematical properties of the Calderon operator [31, 95]. However
32

both the above approaches are very involved. The former approach of designing loop
and star functions require scaling the basis functions depending on the frequency and
will not work for truly multi-scale problems, where a large range of discretization
scales (ranging from 10−5 λ to 10−1 λ) are involved. Further the loop-star scheme
cannot be used for multiply connected domains [96]. The latter approach is applicable
to multi-scale problems but involve the careful construction of intermediary function
spaces to handle the Calder´n projector [97]. .
o
These shortcomings of the traditional RWG based discretization technique is the
chief motivation for the development of the Generalized Method of Moments. In the
remainder of this chapter, we will describe a mechanism for the construction of an
“umbrella” framework that allows for the inclusion of a wide variety of basis functions.
This framework will allow us to decouple the discretization of the geometry from the
definition of the basis function. Starting from this decoupling we will describe a
scheme that allows the easy mixing of various classes of functions in the basis space.
The scheme described herein will allow for the trivial transition from one order of
basis function to another. We will further describe a choice of basis functions that
alleviates the ill-conditioning problem at low frequencies.
In order to remove the usually close knit relationship between the basis functions
and the tessellation of the geometry, the GMM methodology follows a sequence of
steps in the design of basis functions for the solution of (2.42a) and (2.42b). In the
rest of this chapter, we will discuss each of these steps in detail.

4.1

Partition of Unity functions and Overlapping
Patches

The GMM procedure starts by replacing the regular tessellation of the geometry
with a set of overlapping patches as shown in Fig. 4.2. These patches provide for a
33

Ωj

Ωk

Ωi
Figure 4.2: Construction of GMM patches (Ωi ) shown as shaded region.
mechanism of modelling the surface of the geometry but at the same time, via the
nature of their overlap, allow for the mixing of different types of functions.
Definition 2. Open Cover Let Ω ∈ R3 be an open set. Define an open cover of the
domain Ω, {Ωi } such that

i Ωi = Ω that satisfies a point-wise overlap condition

∃M ∈ N

∀r ∈ Ω card i|r ∈ Ωi ≤ M

Once these patches are defined, we require a scheme to provide for the smooth
transition from one basis function to the next. This is effected through the use of a
partition of unity function. Mathematically, we provide the following definition.
Definition 3. Partition of Unity
Let

ψi

be a Lipschitz partition of unity subordinate to the cover

the following conditions:
(i) supp ψi ⊂ Ωi for all i;
(ii)

i ψi ≡ 1 on Ω;

(iii) ||ψi || ∞ 3 ≤ C∞ ; and
L (R )
C
.
(iv) || ψi || ∞ 3 ≤ Cd =
diam(Ωi )
L (R )
34

Ωi

satisfying

Here C∞ and C

are two constants. Then ψi is called an M, C∞ , C

of unity subordinate to

partition

Ωi . The partition of unity is said to be of degree k if

ψi ∈ C k (Rd ). The covering sets Ωi are called patches and they take the place of the
tessellation typically used in standard MoM.
Thus in the GMM framework, the standard tessellation is replaced by (i) a set of
overlapping patches that cover the domain Ω and (ii) a partition of unity subordinate
to this cover. Once this is achieved, we can define basis functions on these partitions
of unity and corresponding Galerkin testing yields a linear system of equations.
This concept of overlapping domains and corresponding PU functions is shown in
Fig. 4.3a. Generalization to higher dimensions can be effected using tensor products
as shown in Fig. 4.3b. Higher order extensions of these functions can be trivially
defined [98]. Thus, by definition, the PU functions and overlapping domains ensures
continuity across patch boundaries.
Once these patches and the partition of unity is defined, basis functions are defined
on these patches as a product of the partition of unity and an approximation function
in each patch. Due to the fact that the continuity is enforced using the partition of
unity, the approximation function can be chosen completely arbitrarily.

4.2

Definition of Basis Functions

Vector basis functions in each patch Ωi are defined as a product of the PU function
ψi (r) and an appropriate approximation function of order m, ν m (r). For a GMM
i
discretization, the approximation JN to the current J can then be written as
Definition 4. GMM Basis Functions and Approximation Spaces
N
J ≈ JN =
i=1

.
ji =

N
i=1 m

35

ai,m ψi ν m .
i

(4.1)

(a) 1-D
Figure 4.3: Examples of linear PU functions in (a) one and (b) two dimensions

36

(b) 2-D
Figure 4.3 Continued: Examples of linear PU functions in (a) one and (b) two dimensions
Then the local approximating space in each patch Ωi is given by
Vi = spanm ψi ν m .
i

(4.2)

Note that within each patch Ωi , there is no restriction on the set of functions ν m .
i
This set can be any mix of functions such as simple higher order complete polynomials,
functions that best approximate the local physics, etc. The construction of a scalar
basis function using a product of a two-dimensional linear PU function (shown in Fig.
4.3b) and a tensor product of higher order polynomials (in Fig. 4.4a) is shown in Fig.
4.4. For ease of presentation, the summation over m will be implicitly assumed and
the corresponding indices will be dropped in the rest of this work.
With this definition of the basis function, we can now state an important theorem
about the approximation errors introduced by the GMM scheme.

37

(a) ψi (x)
Figure 4.4: Construction of scalar GMM basis functions. (a) the GMM approximation
function and (b) the product of the approximation function and the partition of unity
function.

38

(b) ji (x) = φi (x)ψi (x)
Figure 4.4 Continued: Construction of scalar GMM basis functions. (a) the GMM
approximation function and (b) the product of the approximation function and the
partition of unity function.

39

Theorem 4.2.1. For the quantities, M , C∞ and C as defined above, let J ∈
H 1/2 (Ω) [37] be a function to be approximated. Assume that on each patch Ωi ∩ Ω,
J can be approximated by a function ji ∈ Vi such that
J − ji

L2 (Ωi ∩Ω)

≤ ε1

· {J − ji }

L2 (Ωi ∩Ω)

≤ ε2

× {J − ji }

L2 (Ωi ∩Ω)

≤ ε3

max
i
max
i
max
i

then the following holds for JN as defined in (4.1).
J − JN

≤
L2 (Ω)

· (J − JN )

≤
L2 (Ω)

× (J − JN )

≤
L2 (Ω)

√
M C ∞ ε1
√
√

(4.3)

2M (Cd ε1 )2 + (C∞ ε2 )2

1/2

2M (Cd ε1 )2 + (C∞ ε3 )2

1/2

(4.4)
(4.5)

Proof. The proof is presented here for (4.4) and (4.5) and that for (4.3) follows trivially. Using the basis function definition in equation (4.1), it can be shown that

× {J − JN } =

× {J −
×{

=
=

(

ψi ν i }
i
ψi (J − ν i )}

i
× {ψi (J − ν i )})

i
( ψi × {J − ν i } + ψi

=

× {J − ν i })

(4.6)

i
where the summation over the order m is assumed. To prove (4.4), note that the
identity in equation (4.6) holds identically for the divergence operator as well. Given
40

equation (4.6) above, it can be easily verified that
=
× {J − JN } 2 2
L (Ω)

( ψi × {J − ν i }
i
× {J − ν i }) 2 2
L (Ω)

+ ψi

ψi × {J − ν i }

≤ 2

2
L2 (Ω)

× {J − ν i }

2
L2 (Ω)

i
+ 2

ψi
i

ψi × {J − ν i }

≤ 2

2
L2 (Ωi ∩Ω)

× {J − ν i }

2
L2 (Ωi ∩Ω)

i
+ 2

ψi

i
≤ 2M 2 Cd 2 {J − ν i } 2 2
L (Ωi ∩Ω)
2
+ 2M 2 C∞
× {J − ν i } 2 2
L (Ωi ∩Ω)
≤ 2M 2 (Cd ε1 )2 + (C∞ ε3 )2

(4.7)

Theorem 4.2.1 essentially states that the total error in the approximation to the
current is bounded above by the l ocal error in the approximation. This result is
powerful in the sense that it allows for a direct way to control the approximation
error, by careful choice of approximation function.
As is evident from the above discussion, the GMM provides a framework to include
different approximation functions without explicitly requiring continuity constraints.
The order of the PU function determines the order of continuity between patches [99].
Thus, the global approximation inherits the smoothness of the PU function. An advantage offered by the GMM framework is the freedom of choice in the approximation
function. By disassociating the definition of the approximation function from the underlying tessellation, GMM allows for the mixing of multiple types and orders of
41

functions. While the approximation function can be chosen practically arbitrarily, we
describe a mechanism to choose functions that are suited to approximating surface
vector functions.

4.3

Design of the Vector Approximating Function

To define appropriate approximation functions ν m for the vector bases, we start with
i
a surface Helmholtz decomposition of the current. For any vector field J(r) defined
on the surface Ω, there exists two scalars fields φ(r) and ζ(r) such that

J(r) =

ˆ
s × n(r)ζ(r).

s φ(r) +

(4.8)

Assume that we have an approximation ν m to the local behavior of the current J(r)
i
on patch i.

Definition 5. Assume the knowledge of a local approximation, ν i to J such that on
patch i

• maxi

νi − J

L2 (Ω∪Ωi )

• maxi

s · νi −

• maxi

s × νi −

≤ ε1

s · J L2 (Ω∪Ω )
i

≤ ε2

s × J L2 (Ω∪Ω )
i

≤ ε3

Then, define two scalar functions φi (r) and ζi (r) such that
m,c
νi
=

m
s φi (r)

m,d
ˆ
νi
= ni (r) ×

m
s ζi (r)

42

(4.9)

and consistent with (4.1) define the approximation to the current as
N
JN =

i=1 q={c,d} m

m,q
ai,q,m ψi ν i .

(4.10)

Again, we allow for summation over order for each approximation function m.
Given the Definition 5 above and the definition of M , C∞ and Cd in Definitions
(2) and (3), the following result is a direct consequence of Theorem (4.2.1)
Result 4.3.1. For the local errors defined in (5) above the global errors in the approximation are bounded as follows
• J − JN

≤
L2 (Ω)

•

· (J − JN )

•

× (J − JN )

√

M C ∞ ε1

L2 (Ω)

≤

√

≤
L2 (Ω)

1/2
2M (Cd ε1 )2 + (C∞ ε2 )2

√

1/2
2M (Cd ε1 )2 + (C∞ ε3 )2

Fig. 4.5a and 4.5b show the basis functions constructed using third order Legendre
m
polynomials for φm (r) and ζi (r). These are obtained by multiplying the curl-free
i
and divergence-free components of the approximation functions with the corresponding PU function for the patch. Note, however, that the local approximation is only
quasi-curl and divergence free.
It is evident from the above exposition that there is no restriction on the functions
m
φm and ζi . Indeed, if the local asymptotic behavior is known it can be readily
i
incorporated into the approximation function space. Specifically, assuming that the
local behavior of current (j(r) ≈ ν i (r)) supports at least one derivative, a possible
way to compute the functions φi (r) and ζi (r) is given by the following lemma.
Lemma 4.3.1. Consider a vector function ν i (r) defined on Ωi . Let φi and ζi be
43

(a) ψi ν m=3,q=c
i
Figure 4.5: Vector GMM basis functions based on the quasi-surface Helmhotlz decomposition basis functions. (a) The quasi curl-free basis functions and (b) the quasidivergence free basis functions. Both functions are of order p = 3.

44

(b) ψi ν m=3,q=d
i
Figure 4.5 Continued: Vector GMM basis functions based on the quasi-surface
Helmhotlz decomposition basis functions. (a) The quasi curl-free basis functions
and (b) the quasi-divergence free basis functions. Both functions are of order p = 3.

45

solutions to
2 φ (r) =
i

· ν i (r)

2 ζ (r) = n(r) ·
ˆ
i

× ν i (r)

(4.11)
(4.12)

then
νi =

s φi (r) +

ˆ
s × ni (r)ζi (r) .

(4.13)

Scharstein [100] provides the proof of the above result along with some detailed
insight into the physical implications of equations (4.11) and (4.12). We have previously shown for a class of surfaces, that once appropriate approximation and basis
functions are defined the error propagated through the discretized integral operator
is bounded above by the local error in approximation [101]. This section provides a
mechanism to choose an approximation that minimizes this error using all available
knowledge of the local behavior of the function. In the next section, we will examine
the effect of this approximation error on propagation through the integral equation.

4.4

Error Bounds for the Discrete Integral Operators

For completeness, we state the following standard theorem from integral operator
theory. The proof can be obtained from [28].
Theorem 4.4.1. Any integral operator with a weakly singular kernel is compact and
therefore bounded.
From this theorem and the definition of the EFIE and MFIE operators in (2.42a)
and (2.42b), the following lemma immediately follows.

46

Lemma 4.4.2. The MFIE operator defined by
. 1
ˆ
ˆ
dr
n(r) ×
n(r) × Hs (r) = K ◦ {J(r)} =
4π
Ω
where g(R) =

g(R) × J(r )

(4.14)

.
exp −jκR
and R = |r − r |, is bounded. The EFIE operator, as defined
R

by
jη
. jkη0
Es (r) = T ◦ {J(r)} = −
dr g(R)J(r ) − 0
dr
4π Ω
4πk Ω

g(R) · J(r ) (4.15)

is not necessarily bounded.

In what follows, we will consider the MFIE and EFIE in sequence. The following
theorem holds directly for the MFIE.
1/2
1/2
Theorem 4.4.3. For any linear and bounded operator K : H
(Ω) → H
(Ω)
div
div
1/2
ˆ
with kernel GM (r, r ) := n(r) × g(R) and functions J, JN ∈ H
(Ω) and t ∈
div
1/2
H
(Ω) ∃ a constant CM such that
div

Ω

K ◦ {JN (r )} ti (r)dr −

Ω

K ◦ {J(r )} t(r)dr ≤ CM J − JN 2 (4.16)

Proof. The proof for the above follows from the Cauchy Schwartz inequality and the
1/2
observation that u (1/2),div ≤ u 2 ∀u ∈ H
.
div

Ω

Ω

GM (r, r ) × JN (r )dr

t(r)dr −

=
Ω

GM (r, r ) × J(r )dr

t(r)dr

GM (r, r ) × JN (r ) − J(r ) dr

t(r)dr

Ω
Ω

Ω

(4.17)
47

Applying Cauchy Schwartz inequality we obtain

Ω
≤
Ω

Ω

GM (r, r ) × JN (r ) − J(r ) dr

GM (r, r ) × JN (r ) − J(r ) dr

2 Ω

t(r)dr
t(r) 2 dr

(4.18)

1/2
Since, t(r) ∈ H
(Ω) we have that
div

Ω

t(r) 2 dx ≤ C < ∞

(4.19)

Further we have from the definition of boundedness, that

Ω

GM (r, r ) × JN (r ) − J(r ) dr

2

≤ C1

JN (r ) − J(r )

2

(4.20)

Combining equations (4.19) and (4.20) we have that

Ω

Ω

K ◦ Jn (r ) − J(r ) dr

t(r)dr ≤ CM

JN − J

2

(4.21)

where CM := CC1 .

These theorems demonstrate that, in each patch, the error in approximation
bounds the error propagated through the integral operator. Note that the theorems
hold only for compact operators. It is not directly applicable to unbounded operators
such as the EFIE. However, this can be overcome by casting the EFIE solution in a
Galerkin testing framework. First, the bilinear form of the EFIE as follows

Definition 6. For the EFIE operator, T given by
jη
. jkη0
T ◦ {J(r)} = −
dr g(R)J(r ) + 0
dr
4π Ω
4πk Ω
48

g(R) s · J(r )

Let BE be the bilinear form defined by
BE (J, t) =< t, T ◦ {J} >

(4.22)

Assuming that the domain Ω is smooth enough to permit the local basis functions
to support at least one derivative, it can be shown that

BE (J, t) =

T ◦ {J(r)} · t(r)d(r)
Ω

jkη0
dr t(r) ·
dr g(R)J(r ) +
4π Ω
Ω
jkη0
t(r) · g(R)J(r )dr dr −
=−
4π Ω Ω
=−

jη0
dr t(r) ·
dr g(R) s · J(r )
4πk Ω
Ω
jη0
· t(r)g(R) s · J(r )dr dr
4πk Ω Ω
(4.23)

Definition 7. The above equations give the following splitting of the Galerkin tested
EFIE operator
.
< t, T ◦ {J} > = < t, Tc ◦ {J} > + < tc , Tq ◦ {Jc } >
ikη
iη
.
= − 0 < t, Gk (r, r ) ◦ {J} > − 0 <
4π
4πk

· t, Gk (r, r ) ◦

(4.24)
·J>

It can be shown, using Theorem 4.4.1, that the operators Tc (r, r ) and Tq (r, r )
.
with kernel Gk (r, r ) = g(R) are compact and therefore, bounded. Therefore, Theorem 4.4.3 is directly applicable to the first term in equation (4.24). Further a scalar
specialization of 4.4.3 can be applied to the second term in equation (4.24)

Theorem 4.4.4. For any linear and bounded operator Tq : H 1/2 (Ω) → H 1/2 (Ω)
with kernel g(R) and (scalar) functions J, JN ∈ H 1/2 (Ω) and t ∈ H 1/2 (Ω)
49

∃

a

constant CE such that

Ω

Tq ◦ {JN (r )} ti (r)dr −

Ω

Tq ◦ {J(r )} t(r)dr ≤ CE j − jN 2 (4.25)

Combining Theorems 4.2.1 and 4.4.3 (or equivalently, 4.4.4), we have the following
results for both the EFIE and the MFIE.

Result 4.4.1. For local approximating functions JN that are at least C 1 (locally),
the approximation error on each patch bounds the error propagated through the integral equation. If we assume the existence of a convergent inverse for the discretized
equation, then the error in the solution will also be bounded by this error.

The local continuity assumption used in the result above is not restrictive since
all it implies is that the domain needs to be partitioned such that each patch does
not contain local singularities. In practice, this can be easily achieved by matching
patch boundaries with any geometric singularities. One can then define functions
1/2
JN ∈ H
as a set of functions that are piece-wise C 1 . Thus, finally the above
div
1/2
result along with the observation that if J, JN ∈ H
, then J − JN 2 ≤ C < ∞
div
makes it possible to obtain L2 bounded errors by choosing basis functions from the
trace space. The above theorems provide the foundation for the development of the
GMM and can be extended to any space of functions as long as the partition of unity
framework is maintained. This freedom of choice is the inherent advantage of the
GMM scheme as opposed to the standard MoM scheme where the basis functions are
tightly coupled to the discretization. For example, it can be shown that the scalar
higher order polynomial functions can be easily cast in this framework and error
bounds can be derived.
50

4.5

Benefits of the GMM framework

The GMM method described in the previous subsections comprises of two components; (i) defining overlapping domains and partitions of unity on these domains
and (ii) defining approximations functions in the partitions of unity. The definition
of basis functions provides a mathematically rigorous framework for designing local
approximation spaces. From an intuitive perspective, we have always used the PU
functions whenever the basis function space has comprised of pulses or hat functions for the analysis of scalar problems. Likewise, basis functions that are commonly
used in the computational electromagnetics community (eg. Rao-Wilton-Glisson [71],
Buffa-Christiansen [102], characteristic basis [103] or Graglia-Wilton-Peterson [74])
can all be interpreted as a product of a partition of unity function and higher order
approximating functions (with the correct definition of the PU domain). In most
cases, the PU function is implicitly assumed and is typically a pulse function. When
it is a pulse function, one has to worry about means to achieve conformal behavior
between neighboring elements. As in GFEM [99], the use of basis functions that
explicitly comprise of a product of an overlapping partition of unity and approximating functions removes the need for such boundary treatment between different sub
domains, and dramatically increases the number of functions that can be used for approximation. For instance, one can use a third order approximation in a domain next
to one where a first order is being used, without explicitly enforcing any condition.
Or for that matter, mixture of polynomial and non-polynomial basis functions. Use
of the these definitions of basis will be demonstrated in the Section 4.9. In the next
Section, we present specifics of the implementation of the general GMM scheme for
arbitrary 3D scatterers.
51

Ωj

Ωk

Ωi
Figure 4.6: Construction of GMM patches (Ωi ) shown as shaded region. The patches
are constructed as a set of triangles joined at a node.

4.6

Implementation of the GMM

The first need is to define the domain of support of each of these basis functions. This
is done by partitioning the domain Ω into subdomains/patches Ωi that overlap and
completely cover the domain. This construction of overlapping patches is illustrated in
Npat
Fig. 4.6. Formally, ∪i
Ωi = Ω, where Npat is the number of patches/subdomains.
Note, that at this point, a patch can be defined either within a meshless or a meshed
framework. The key requirement is that one have a surface description for each patch
Ωi . At this time, for a wide range of problems such a description is most easily
available as a surface tessellation. As a result, the rest of the development will focus
on developing the GMM scheme and basis functions within a meshed framework.
Parenthetically, we do observe that extending this to a meshless framework is not
very difficult. As, objects are typically represented using a triangulation of their
surfaces, that shall be our starting point.

4.7

Construction of Overlapping Patches

To define a set of overlapping patches starting from a triangular mesh, we designate
each node as the “center” of a patch, defined as the union of all triangles that share
52

the node.
Definition 8. Overlapping Patches
Consider a set of nodes NL = L {Ni } and a corresponding triangulation ∆N =
i=1
N {∆ } of a domain Ω. Let {N
n
nj=1−3 } be the set of nodes associated with
n=1
triangle ∆n . For each node Ni , denote a corresponding patch Ωi defined by a set of
triangles; Ωi = {∆n : ∆n ∈ ∆N , Ni ∈ Nnj }. The GMM covering is defined as the
set of these (overlapping) patches {Ωi }. The overlap between any two patches Ωp and
.
Ωq is another set of triangles ∆p,q = Ωp ∩ Ωq
This definition of GMM patches from a union of polygonal shapes can also be
seamlessly translated to higher order geometries defined on triangulations described
in [74]. In the remainder of this section, ∆n will be used to represent both flat
and equivalent higher order triangulations of the surface. The above patch design
also allows for integration on the triangles described by the original tessellations and
therefore, allows the use of available integration rules.
The next step is a consistent definition for basis functions on these patches. Note
that these patches are no longer coplanar and do not correspond to any canonical
geometry (unlike triangles, quads, etc.). This necessitates the careful definition of
functions on Ωi so as to maintain continuity. To this end, we first state the following
lemma for continuous portions of the surface.
Lemma 4.7.1. Consider a patch Ωi . Let the surface Ω ∩ Ωi be at least C 1 (R3 )
smooth. Let the location of Ni be denoted by ri . Let the surface normal at any point
ˆ
r ⊂ {Ω ∩ Ωi } be denoted by n(r)(consistent with the previous definitions). Then
by the definition of continuous surfaces, we have ∀ δi ≥ 0 ∃ Bi,ε ∈ R3 such that
∀ r ∈ {Bi,ε ∩ Ωi }

n(r) − n(ri ) 2 ≤ δi where Bi,ε is a ball of radius ε centered at

Ni
With this preamble, we can now define a patch normal and projection plane for
53

each patch.
Definition 9. Patch Normal and Projection Plane
Consider a patch Ωi composed of m triangles Ωi = m {∆k }. Define an average
k=1
m n(r ), where r are points chosen on
ˆ
ˆ
normal for the patch, ni as ni,ε = 1/m k=1 ˆ k
m
each of the triangles that make the patch such that rm ⊂ Ωi ∩∆k and rk − ri 2 ≤ ε.
The projection plane for patch Ωi , Pi is defined as the plane passing through ri and
ˆ
normal to ni .
For purposes of illustration, Fig. 4.7 shows the construction of one patch and
ˆ
its corresponding projection plane.Observe that the patch normal ni coincides with
the true surface normal at ri if the triangulations are smooth (and necessarily higher
ˆ
order), if not, the differences between the patch normal ni and the normals of the
individual member triangles provide an intuitive measure of the “sharp-ness” of the
surface at ri . This information can be used as a criteria for the design of approximation functions. If the surface is discontinuous, it is partitioned into areas that are
piecewise continuous and then the above procedure is applied to each portion.
Once the normal and the projection planes are defined for each patch, the basis
functions are defined on the projection of the entire patch on the plane using a local
coordinate system on the projection plane.
Definition 10. Local Coordinate System
Let Γi be the projection of a patch Ωi on Pi . For the rest of the discussion, primed
coordinates denote the projection of the corresponding global coordinate on Pi (note
that, by construction, ri = ri ). Two local projection vectors are generated for each
.
ˆ
ˆ
ˆ
patch up,i (r ) and vp,i (r ) as up = up / up 2 where up = r −ri such that up 2 ≥
ˆ
ˆ
ˆ
r − ri ∀r ∈ Γi and vp = ni × up .
2
With this preamble, it is now possible to define a mapping from regular x, y, z
coordinate space into the u, v coordinate space defined on the projection Γi as follows:
54

(a) Definition of GMM Patch Ωi

ˆ
n

(b) Construction of Projection Plane
Γi
Figure 4.7: Construction of (a) GMM Patches (Ωi ) shown as shaded region and (b)
Projection Planes (Γi ) from standard triangulations. The patches are constructed as
a set of triangles joined at a node.

55

• ∀r ∈ Ωi ∃ r ∈ Γi ,
ˆ
ˆ
• define u = r · up and v = r · vp , and
• then the local coordinates u and v are given by uˆ p and vˆ p .
u
v
Approximation and basis functions are now constructed by defining functions on
the corresponding planar coordinates. One consequence of the choice of the coordinate
system is that the projection coordinates span an area that contains, but is not limited
to, Γi . This is possible as functions defined on these patches are eventually multiplied
by a PU function that goes to zero at patch boundaries.
Given the definitions above, it is important to ensure that functions that are
continuous in the projection domain are so in the real space as well. The following
theorem (stated without proof) prescribes the necessary conditions.
Theorem 4.7.2. Given a surface description Ωi (r) that is at least C 1 (R3 ) smooth,
consider any function F (u(r), v(r)) defined on the projection Γi . Then, if F ∈
C m (Γi ) then F ∈ C m (Ωi ) for any m.
The proof follows from the definition of partial derivatives. While a corresponding
result does not hold for piecewise continuous triangulations, we have observed, in
practice, that it is possible to construct functions that are continuous across the
interior edges of the patch. This fact is shown in Fig. 4.8 where a set of polynomial
functions are defined on a flat projection plane, in Fig 4.8a, and then transformed to
a non-flat, piecewise linear, triangulation in Fig. 4.8b. It is apparent from the figures
that there are no discontinuities introduced due to the projection.

4.8

Computation of Matrix Elements

Next, the discrete versions of the EFIE and the MFIE are constructed using Galerkin
testing. From the definition of the discretized EFIE and the MFIE in the previous
56

(a)

(b)
Figure 4.8: Basis functions defined on (a) the projection plane and (b) the true
geometry (surface of a tessellated sphere). The figure (b) shows that there is no
discontinuities introduced by the projection.

57

chapter, it can be verified that the computation of the matrix elements involves the
evaluation of the following three types of integrals; viz.

I1 =
I2 =
I3 =

Ωi
Ωi
Ωi

dr˜i (r) ·
fn
dr˜i (r) ·
fn
dr˜i (r) ·
fn

k
dr gq (R)fj (r )

Ωj

k
gq (R) × fj (r )

dr
Ωj

k
gq (R) · fj (r )

dr
Ωj

(4.26a)
(4.26b)
(4.26c)

n
n
ˆ
where ˜i (br) can stand for either fi (r) or n(r) × fi (r). There is extensive literature
fn
on evaluating both the singular and hypersingular integral that occur in Eqns. (4.26a)
and (4.26b) [104–106]. The methods that we have used follow those presented in these
papers. In what follows, we shall deal with terms of the form (4.26c). The evaluation
of (4.26c) is more involved. Assume that Ωi(j) is the union of a set of triangles
denoted by ∪∆i(j) . Then (4.26c) can be rewritten in either of the following two
forms:

dr t(r) ·
Ωi
∪∆i

=−
+

dr t(r) ·

∪∆i

∪∂∆i

gq (R) · j(r ) =

dr
Ωj

dr

∪∆j

s · t(r)

∪∆j

ˆ
dr mi (r) · t(r)

gq (R) · j(r )

dr

dr

∪∆j

58

dr

s gq (R) · j(r )
s gq (R) · j(r )

(4.27a)

or

=−
+
−

∪∂∆i

∪∂∆i
+

∪∆i

dr

s · t(r)

ˆ
dr mi (r) · t(r)

ˆ
dr mi (r) · t(r)

∪∆i

∪∆j

dr

s · t(r)

dr gq (R) s · j(r )

∪∆j

∪∂ ∆j

∪∂ ∆j

dr gq (R) s · j(r )

(4.27b)

ˆ
dr gq (R)mj (r ) · j(r )

ˆ
dr gq (R)mj (r ) · j(r )

ˆ
ˆ
where mi and mj are the unit normals to each boundary ∂∆i and ∂ ∆j of ∆i and
∆j , respectively. The divergence operator ( ( ) ) can be reduced to the corresponding
()
surface divergence operator ( s ) because j(r)(t(r)) are surface functions. In typical
ˆ
ˆ
RWG settings, the normals mi = −mj and so, very conveniently, the last three terms
in (4.27b) drop out. This leaves an integrable singularity using the above mentioned
techniques for the magnetic vector potential.

However, in GMM the set of triangles forming a patch are not co-planar. This
means that the normal mi(j) is not continuous within a patch and so, we need to
compute all the terms on the right-hand side of the expansion in (4.27b). On the
outer boundary of the patch though, the integral will vanish as a result of the PU
function and thus only interior triangle edges need to be considered. (4.27a) involves
integrals that are an order of singularity higher than (4.27b). However, as r → r
(gq is singular), (4.27b) can no longer be used since (1/R) is not integrable on a line
when r → r . As a result, in the GMM implementation, this is handled depending
on the location of r and r .
59

r and r lie on the same triangle
In this case, we employ a singularity subtraction technique which takes into account
the measure provided by both the source and testing integral, in order to mollify the
second order singularity. This technique was first proposed by [107] and computes
the double integrals in (4.27a) together rather than first computing the source and
then the testing integral as is usually done in practice.

r and r lie on triangles that share an edge
In this case, the integral is not singular but the 1/R3 term needs to be handled
carefully. This case is evaluated using (4.27a). The first term in the expansion can
be computed for near singular (but non-singular) terms using a slight modification of
the method described in [104] to integrate terms of the type 12 .
R

r and r are not in neighbouring triangles but are close (|r − r | ≤ 0.15λ)
Here we use (4.27b) to evaluate the integral since the line-line singularity is no longer
a concern.

r and r are sufficiently separated from each other
In this case, the (now) regular integral can be evaluated using either of the two
expressions in (4.27a) or (4.27b), using standard numerical quadrature.

A transformation for the near singular integration of

g(R)j(r ) for r → r

Consider the evaluation of the following integral for r → r but r = r :

dr

g(R)j(r )

∆
60

(4.28)

We modify the steps outlined in [104] for the evaluation of

dr g(R)j(r )

(4.29)

∆
In [104], the authors split the triangle into three parts around the projection of r on
to the tangent plane and use a geometric transformation to transform the integral to
the following form
3

hi

i=1 0

dyi

xi (yi )

dxi j(r (xi , yi ))eikR

1
2
(x2 + yi + z 2 )
i

(4.30)

where z is the height of the r above the equivalent tangent plane of the triangle ∆ ,
hi and xi are the height and the base length of each sub triangle. [104] further uses
an sinh−1 transformation to transfer the singularity of 1 x2 + y 2 + z 2 from the
integrand to the limits of the integral. It can be shown that a similar transfer of singu2
larity can be achieved for 1 |r − r |2 = 1 x2 + yi + z 2 using a tan−1 transformation.
i
The rest of the steps remain identical.

4.8.1

r and r are not in neighbouring triangles but are close
(|r − r | ≤ 0.15λ)

Here we use (4.27b) to evaluate the integral due to the (now) better behaved nature
of the singularity.

4.8.2

r and r are sufficiently separated from each other

In this case, the singularity is no longer a concern and the integral can be evaluated
using either of the two expressions above.
61

4.8.3

Integration with acceleration techniques

At this point, we simply note that once the basis functions are defined and the
appropriate integration rules chosen, integration with the wideband fast multipole
method (FMM) is straightforward [38, 48, 108].

4.9

Numerical Experiments

Next, via a series of numerical experiments, we will demonstrate that the GMM
discretization framework proposed here (i) provides excellent approximation of the
unknown currents and its derivatives and (ii) accurate results for scattering from
several topologically different objects, (iii) demonstrate the viability of this method
to accurately analyze scattering from objects using a mixture of basis functions, (iv)
demonstrate analysis from objects meshed with non-conformal meshes and (v) low
frequency stability the choice of basis. Where necessary, results have been obtained
using a GMM code augmented by FMM.

4.10

Quality of Approximation

The first goal is to test whether multiple functions can be seamlessly mixed together
within the GMM approximation space. This will be verified by a series of numerical
experiments described next. It is assumed that the current in a domain Ω satisfies
some prescribed functional dependence, say J(r) = f (r). The coefficients of these
basis functions
N
JN (r) =

ai ji (r) =
i

i=1 q=d,c m
62

ai,m ψi ν m (r)
i,q

(4.31)

are computed by solving a discrete system
¯a ¯
Z¯ = b

(4.32)

where Zi,j =< ji , jj > and bi =< ji , J > The approximation error defined as
N aj
i=1 i i
|f (r)|

f (r) −
ε(r) =

(4.33)

is computed by oversampling the currents on the geometry. Similar errors are computed for both the curl and the divergence of the current using the coefficients {ai }.
For each test described next, the approximation was computed on a square domain
of side 1.0m using overlapping patches, each of size h = 0.2m. Fig. 4.9a-4.9b show
the convergence of the maximum approximation error with approximation order p for
each of the tests.

4.10.1

Approximation of Polynomial Functions

It was earlier argued that the GMM scheme permits the mixing basis functions of
various types. For this to be effective, it is important that the higher order functions
do not corrupt the approximation of smooth functions. Fig. 4.9a shows the converˆ
gence of the approximation error for two functions (i) a constant f (r) = x and (ii) a
ˆ
ˆ
polynomial f (r) = xy 3 + yx3 . As expected, the error is consistently less than 10−12
for the constant function. Similarly, it is very interesting to note here that once the
order of the approximation reaches the order of the polynomial function (p = 3), the
error in the GMM scheme drops to machine precision.
63

ˆ
ˆ
ˆ
(a) f (r) = x (m = 0) and f (r) = xy 3 + yx3 (m = 3)
Figure 4.9: Quality of approximation of (a) polynomial function and (b) sinusoidal
functions.

4.10.2

Approximation of sinusoidal functions

ˆ
Fig. 4.9b shows the convergence graphs for f (r) = x sin(x) sin(y) and f (r) = x exp(−jk·
ˆ
r). The figures show the convergence curve for both the approximation in the current
and its curl and divergence. As seen in the figures, the error in both

· J and

×J

is the same and of order (p − 1), as expected.

4.10.3

Approximation of non-smooth and near-singular functions

The true power of the partition of unity comes into play when the functions are
discontinuous across the patch boundaries. To demonstrate this, consider a function
that is first order discontinuous on a line that does not coincide with the patch
64

ˆ
(b) f (r) = x sin(x) sin(y) and f (r) = x exp(−jk · r)
ˆ
Figure 4.9 Continued: Quality of approximation of (a) polynomial function and (b)
sinusoidal functions.

65

(a) Error without PU
Figure 4.10: Quality of approximation of discontinuous function: (a) without using a
partition of unity; and (b) with a partition of unity.

boundary described by
∀ x ≤ x0



xx/x0
ˆ

f (r) =



x
ˆ

∀ x ≥ x0

.

(4.34)

Fig. 4.10a shows the approximation error obtained by using two separate functions
in each patch without overlapping patches or a partition of unity. As can be clearly
seen, the error in the function propagates through the patch because the discontinuity
does not coincide with the patch boundary. Fig. 4.10b on the other hand, shows the
same function approximated by overlapping patches and a partition of unity. As
the figures clearly show, the error is significantly better in the GMM basis function
framework.
Next, we examine the efficacy of using a GMM framework in approximating cur66

(b) Error with PU
Figure 4.10 Continued: Quality of approximation of discontinuous function: (a) without using a PU and (b) with a PU.

67

rents near a one dimensional edge upto a distance 1m from the tip; see [79, 80] and
references therein for details of some of the basis functions in current use. Fig. 4.11
demonstrates the error convergence as a function of the number of unknowns for
three different types of basis functions that are detailed next. Typically, the current
perpendicular to the edge is modeled using leading order terms from a Meixner series
√
expansion (as 1/ d where d is the distance from the edge). In [79], this dependence
is incorporated by changing the basis function in the patch near the edge. Specifically, in [79] the domain is partitioned into segments of size h; in the patch nearest
the edge, the basis function is chosen to be the leading term of the Meixner series,
and in all other patches it is either a pulse or a hat function. In figure 4.11, we
examine the error convergence, as a function of the number of unknowns, for the
approximation to the current near a one dimensional edge. Fig. 4.11 shows three
error curves, the first is obtained using the basis functions in [79]. In this case, the
√
basis function in the first patch is 1/ d where d is the distance from the edge/tip.
Starting from the second patch onwards hat functions are used as basis. This result
establishes a baseline for comparison of convergence. The second curve is obtained
using basis functions that are a combination both singular functions and higher order
polynomials. Specifically, patches comprise of two adjacent segments. On patches
that lie within a perpendicular distance of 0.1m of the edge, it is assumed that basis
√
functions are a superposition of singular functions of the form 1/ d where d is the
distance from the origin and polynomials that are complete to order p = 2. As is
evident from fig. 4.11, this class of basis functions provides better approximating
properties. The third curve in fig. 4.11 demonstrates the convergence obtained by
using a partition of unity function. As before, patches are defined to comprise of two
adjacent segments. For all patches whose perpendicular distance is less than 0.1m
from the edge, basis functions comprise of a linear combination of singular functions
√
(1/ d) and polynomials complete to p = 1, both of which multiplied by a linear
68

PU function. For all patches that are further from this point, the basis function is a
product of a linear PU function and a polynomial of order p = 1. The use of a lower
order polynomial ensures that the order of the basis function is p = 2. As is evident,
the approximation provided by this scheme is considerably better than the previous
case. In large part, this result is due to better continuity properties offered by the
PU function; proofs for convergence of approximation to scalar and vector functions
have been provided in [76, 77].
As can be seen from all the examples presented in this Section, for a series of
discretizations, the PU-based approximation consistently performs well. To a large
extent, this is a consequence of the GMM framework which decouples the traditional
link between tessellation and the basis function definition. In other words, it is possible to define dramatically different basis functions in the same patch (or adjacent
patches) without worrying about auxiliary conditions to maintain continuity. Continuity of these functions is provided by the partition of unity function. As a result,
the quality of approximation is better. In the next Section, we further explore the
capabilities of this representation by applying it to compute electromagnetic scattering

4.11

Scattering Results

In this section, we compare the radar cross-section from a variety of PEC objects.
Unless otherwise noted, in all examples that follow, the incident field is assumed to
ˆ
be x polarized and propagating along the −ˆ direction. The object is discretized
z
using a triangulation such that the average edge length is 0.1λ and GMM patches
are constructed using the scheme outlined in this paper. Since the GMM patches
are of arbitrary shape, we define the side of the smallest enclosing cube for each
GMM patch as an equivalent “edge length” (h). In the following examples, h varies
69

Figure 4.11: Comparison of approximation error with and without PU.
from 0.1λ to 0.4λ, where λ is the wavelength of the incident field. Unless specified
m
otherwise, the GMM approximation functions, φm and ζi , were constructed using
i
Legendre polynomial functions complete to order p. Finally, note that the underlying
triangulation used is the same for both GMM and RWG or its higher order variant.
In what follows, we have used the notation “higher order RWG” to denote basis
functions develop by Graglia, Wilton and Peterson [74].

4.11.1

Analytical Validation and p-convergence

Fig. 4.12a compares the RCS of a PEC sphere obtained using GMM against analytical
solutions evaluated using a Mie-series representation. The number of GMM unknowns
(NGM M ) used is 200, 776 and 3080 for spheres of radii 0.2λ, 0.6λ and 1.0λ respectively. As is evident, the agreement between the two is excellent. Fig. 4.12b compares
the convergence of the error as a function of the polynomial order of the basis function
70

used, for both the GMM and higher order RWG [74]. The comparison is performed
for backscattered data obtained from two scatterers (i) a sphere of diameter 1.0λ and
(ii) a NASA almond that fits in a box of dimensions 0.4λ×1.15λ×3λ. The number of
unknowns for the sphere is as follows: NGM M = 964, 2604 and 5208 for p = 1, 2, 3,
respectively. The number of higher order RWG unknowns are NRW G = 4230, 9072,
15552 at p = 1, 2 and 3, respectively. For the almond, NGM M = 948, 4656 and
7560 at p = 2, 3 and 4, respectively and NRW G = 3800, 7980 and 13680 at p = 2, 3
ˆ
and 4, respectively. In all cases, the incident field was x polarized and incident along
ˆ
ˆ
z. The axis of the almond was oriented along x. The convergence is evaluated for
the error between the numerical and analytical solution for the sphere and for the
relative error between successive steps in the case of the almond. In each case, the
approximation functions used for GMM were p-order complete polynomials in each
of the local projection co-ordinate spaces. Correspondingly the basis functions used
for higher order RWG were of the type suggested in [74] and are complete to order
p. The exponential nature of the convergence is clearly evident from the figure. It is
interesting to note that the convergence rate for the GMM scheme seems consistently
better than that for RWG for both the geometries.

4.11.2

Comparison with standard RWG

In this section, the RCS obtained due to scattering from a variety of geometries is
computed using the GMM basis functions and compared against that obtained using
RWG basis functions. In all cases, first order complete polynomials were used to
construct basis functions for the GMM. Similarly first order RWG functions were
used for the RWG scheme.
Fig. 4.13a shows the RCS from a square PEC plate of side λ discretized using 248
effective GMM unknowns (NGM M ) and 228 RWG unknowns (NRW G ); and an open
cavity of dimensions λ × λ × λ (NGM M = 402, NRW G = 1746). The incident field
71

(a) Comparison with analytical solution for PEC sphere of various
radii
Figure 4.12: Validation studies: Scattering from PEC spheres of different radii. Fig.
(a) shows the RCS comparison and figure (b) shows the convergence of solutions with
basis order p.

72

(b) Convergence with polynomial order (Sphere and NASA Almond)
Figure 4.12 Continued: Validation studies: Scattering from PEC spheres of different
radii. Fig. (a) shows the RCS comparison and figure (b) shows the convergence of
solutions with basis order p.

73

ˆ
ˆ
is x polarized and propagates along z. The RCS is computed by an EFIE formulation
discretized using both the GMM and a standard RWG. It is evident that the GMM
result agrees very well with the standard RWG formulation.
Fig. 4.13b compares the RCS obtained due to scattering from a sphere of radius
1.0λ (NGM M = 1736, NRW G = 4320), an EMC benchmark cone-sphere of height
2.6λ and radius 0.5λ (NGM M = 1792, NRW G = 1656) and a NASA almond that fits
in a box of dimensions 0.4λ × 1.15λ × 3λ (1:3:9) (NGM M = 1528, NRW G = 2800).
ˆ
The cone-sphere and the almond are both oriented with their long axis along x. In
ˆ
ˆ
each case, the incident field is x polarized and propagates along the z direction. Both
the GMM and the MoM results are computed using a CFIE. While the sphere and
cone-sphere results are obtained using polynomials of order p = 1, the almond results
are obtained for both p = 1 and p = 2. As is evident from the figures, the agreement
between the GMM and MoM is excellent.

4.11.3

Comparison with measured data

Finally, we consider scattering from a PEC cube of side 0.7λ. The cube was discretized
using NGM M = 320 and NRW G = 1500 unknowns. Measured data has been
reported for scattering from this structure in [109]. Again, electric field polarized
ˆ
ˆ
along x is assumed in the z direction. Fig. 4.14 shows the comparison of RCS
obtained by GMM and RWG against that obtained from experimental measurement.
Again, there is excellent agreement between all three RCS data sets. These numerical
experiments demonstrate the ability of the GMM technique to accurately use first and
higher order polynomial basis functions to compute scattering cross-sections from a
wide variety of geometries.
74

(a) Plate (λ × λ) and Open cavity (λ × λ × λ) (using EFIE)
Figure 4.13: Validation studies: scattering from (a) a PEC plate and an open cavity.
Plot of bistatic RCS (φ = 0). Fig. (b) bistatic scattering RCS from a NASA thin
almond and a PEC conesphere. The RCS is compared at φ = 0 for the almond and
θ = π/2 for the conesphere. Inset figures show the surface current on the almond and
the conesphere.

75

(b) Sphere (λ), Conesphere (2.6λ × 0.5λ) and NASA Almond
(1:3:9) (using CFIE)
Figure 4.13 Continued: Validation studies: scattering from (a) a PEC plate and an
open cavity. Plot of bistatic RCS (φ = 0). Fig. (b) bistatic scattering RCS from a
NASA thin almond and a PEC conesphere. The RCS is compared at φ = 0 for the
almond and θ = π/2 for the conesphere. Inset figures show the surface current on the
almond and the conesphere.

76

Figure 4.14: RCS comparison for a PEC box (0.7λ). The figure compares RCS due
to scattering of a −ˆ directed plane wave polarized along −ˆ . The bistatic RCS is
x
y
compared at φ = 0 against that obtained using an RWG discretization and against
measured data.

77

4.12

Flexibility of the GMM scheme

The results presented thus far demonstrate the approximation quality of GMM basis functions, the convergence characteristics and applicability to calculating scattered fields from a variety of topologically different objects. In addition, the GMM
scheme also provides the flexibility to mix different types of orders and types of
basis functions, and use non-conformal meshes with minimal change to the overall
computational structure. The examples presented in this section demonstrate these
advantages. In all cases analyzed in this section, the incident field is polarized along
the x direction and the field is incident along the −ˆ direction.
ˆ
z

4.12.1

Mixing orders of basis functions

First, consider a PEC box of dimensions λ × 4λ × λ that is oriented along the hatx
direction. This box is meshed at 0.1λ everywhere to obtain an initial triangulation.
GMM patches and basis functions complete to order p = 1 are constructed to obtain the RCS due to illumination by a z directed incident electric field polarized
ˆ
ˆ
along Ex . Next, the same geometry is discretized using an edge length 0.25λ for
0.75λ ≤ x ≤ 2.25λ and of edge length 0.1λ elsewhere. Basis functions of order p = 2
are used in the larger patches while basis functions of order p = 1 are used in the
small patches. RCS computed using this configuration is then compared against that
obtained earlier using a denser (almost) uniform discretization in figure 4.15a. As
is evident, the agreement is excellent. While the agreement between the two RCS
data sets demonstrates the ability of the technique to mix different orders of basis
functions, the ease with which this mixing can be effected within the GMM scheme
is equally important. In our implementation, types and orders of basis functions to
be used in a patch are specified by an input flag. This level of flexibility in choice of
basis functions is unique to the GMM scheme.
78

Fig. 4.15b demonstrates two RCS obtained due to scattering from a PEC almond
that fits in a box of size 0.3λ × λ × 3.2λ. The first RCS is obtained using basis
functions of order p = 1 everywhere and then second is obtained using basis functions
of order p = 2 in all patches within 0.2λ of the tip and p = 1 everywhere else. Again,
the agreement of between the two RCS data is excellent and demonstrates the ability
of the GMM technique to mix different basis functions.
Finally, figure 4.15c shows RCS due to scattering from a PEC box of size 3λ ×
6λ × 3λ; meshed at three levels of discretization as shown in the inset figure. The
discretizations are at 0.1λ, 0.2λand0.25λ. The RCS is computed using polynomials
of orders p = 1, 2, 3 in each discretization range. This RCS is compared against one
obtained using RWG basis functions on a uniform discretization of 0.1λ; and as can
be seen from the figure, both results agree very well with each other.

4.12.2

Mixing types of basis functions

The GMM scheme not only permits mixture of different orders of polynomial basis
functions, but enables mixtures of different types of basis functions. As was mentioned
earlier, the surface Helmholtz decomposition based system using scalar polynomial
functions is only one possible way to define basis functions. In the following example,
we demonstrate that these functions can be effortlessly combined with functions of
the form ceik·r where c is a pilot vector and k is the propagation vector.
ˆ
ˆ
Consider scattering from a PEC box of size λ × 0.5λ × 0.5λ. Three sets of RCS
data are compared in figure 4.16. One is obtained using GMM based polynomial basis
functions of order p = 1 everywhere. The other two use a combination of polynomials
of order p = 1 and plane-wave basis as described next. The plane wave basis functions
m,1
used are defined as follows: ν m for patch Ωi is defined as ν i
= ueikm ·(r−ri ) and
ˆ
i
m,2
νi
= v eikm ·(r−ri ) where u and v are the local coordinate vectors, ri is the center
ˆ
ˆ
ˆ
of patch i, and km is defined as (k · u) cos(mπ/(M + 1))ˆ + (k · v ) cos(mπ/(M +
ˆ
u
ˆ
79

(a) PEC box (λ × 4λ × λ): GMM mixed order (p =
{1, 2}) vs. standard GMM p = 1.

(b) NASA Almond (1:3:9): GMM mixed order (p =
{1, 2}) vs. standard GMM p = 1.
Figure 4.15: RCS comparison for regular GMM and mixed order GMM: (a) comparison of bistatic RCS (φ = 0) for a PEC box. Fig. (b) shows the comparison of bistatic
RCS (φ = 0) for a PEC thin almond. In each case, the incident field is directed along
ˆ
−ˆ and polarized along x and bistatic RCS is evaluated at φ = 0.
z

80

(c) PEC Box 3λ × 6λ × 3λ GMM mixed order (p = {1, 2, 3}) vs.
RWG p = 1.
Figure 4.15 Continued: RCS comparison for regular GMM and mixed order GMM:
Fig. (c)comparison of bistatic RCS (φ = 0) for a PEC box using three different orders
ˆ
of basis functions. The incident field is directed along −ˆ and polarized along x and
z
bistatic RCS is evaluated at φ = 0.

81

ˆ
Figure 4.16: RCS obtained due to scattering of a x polarized plane wave directed
ˆ
along −ˆ from a PEC box (λ × 0.5λ × 0.5λ) directed along x obtained using regular
z
GMM and mixed basis GMM with polynomials and plane waves.

1))ˆ ∀ m ∈ {0 . . . M }. The three cases that are compared are as follows: (i) only
v
ˆ
polynomials, (ii) plane waves (M = 1) in all patches for 0.4λ ≤ x ≤ 0.6λ and
polynomials elsewhere and finally (iii) using plane waves (M = 2) in patches 0.4λ ≤
ˆ
ˆ ˆ
x ≤ 0.6λ and polynomials elsewhere. A plane wave electric field polarized along x − y
ˆ ˆ
is assumed incident along x + y and the RCS is observed for all θ at φ = 0. It is
observed that as the number of plane wave is increased, the RCS obtained rapidly
approaches that obtained using purely polynomials. This example illustrates that
idea that not only can one used non-polynomial basis functions with ease within the
GMM scheme but also use a combination of polynomial and non-polynomial basis
functions.
82

4.12.3

Use of non-conformal meshes

Finally, another advantage of the GMM arises from the initial motivation to separate
the basis function definition from the tessellation of the geometry. Since, unlike the
RWG, the GMM is not reliant on the tessellation; the GMM can be used with nonconformal meshes. It is common in practice to have different vendors mesh different
areas of a geometry and this results in sets of meshes that need to be explicitly
“stitched together” at the boundaries. There has been a considerable interest in
the design of techniques that can handle non-conformal meshes [110]. However, a
consequence of the design of the GMM scheme is that it can implicitly handle these
kinds of meshes without any change to the framework or the implementation. Here
we present one example illustrating this ability. Fig. 4.17a shows a λ × λ plate that is
discretized at two slightly different scales; one at 0.1λ and the other at 0.15λ. These
two halves are put together and a GMM covering is defined on this “non-conformal”
tessellation. Fig. 4.17b shows the RCS obtained using the GMM scheme compared
with the RCS obtained using a regular RWG discretization using a conformal meshing
at 0.1λ of this plate.

4.13

Condition Number of the System

One of the consequences of using GMM is that the resultant discrete system appears
to have a constant condition number over a wide range of frequencies. We have
examined the condition number of the system over a wide range of frequencies for a
variety of geometries. In each case, we use a constant mesh size that corresponds to
an average edge length of 0.1λ at the highest frequency. This is used to construct
the GMM patches (h varies from about 0.1λ to 0.4λ). The frequency is continuously
reduced from 300MHz to 3Hz. The condition numbers of the impedance matrix
obtained using regular RWG, GMM and loop-star basis functions are compared.
83

(a) PEC plate (λ × λ) discretized using two nonconformal meshes

(b) RCS comparison with standard RWG (on conformal mesh)
Figure 4.17: Backscatter RCS comparison for scattering from a non conformally
meshed PEC plate (shown in figure (a)) RCS obtained by using the GMM discretization on a non conformal mesh compared against that obtained by discretizing on a
(standard) conformal mesh for a PEC plate.

84

Figure 4.18: EFIE Condition number comparisons for PEC sphere (0.2), square plate
(1m) and cube(1m) using GMM, RWG and loop star basis.

Fig. 4.18 shows the condition numbers for three geometries (i) PEC plate (size
1m × 1m, NGM M = 248, NRW G = 228) (ii) PEC cube (side 1m, NGM M = 320,
NRW G = 1500) and a sphere (0.2m diameter, NGM M = 384, NM oM = 480). As
expected, the regular MoM discretization is ill-conditioned at low frequencies. The
loop-star scheme and the GMM scheme results in discrete systems with a constant
condition number across a range of frequencies. It is evident from the figures that
the condition number of the GMM scheme stays constant and is always less than or
equal to the condition number of the loop-star scheme.
The above results demonstrated condition number for discretization using linear
basis functions. Fig. 4.19 shows the condition numbers obtained on discretizing
using higher order basis functions for a 0.2m diameter PEC sphere and PEC strip of
dimensions 1m × 0.2m. The approximation functions used are third order complete
85

Figure 4.19: EFIE Condition number comparisons for PEC sphere (0.2m) and strip
(1m × 0.2m) using higher order (p = 3) complete GMM and RWG basis.

and the condition numbers are compared against those obtained using a higher order
RWG scheme using higher order RWG basis functions of the same order. Again, it
is evident that the RWG condition number increases exponentially with frequency,
whereas the GMM condition number remains constant.
The purpose of a low frequency stable system is the ability to handle meshes
involving a large range of discretization sizes. To test the efficacy of the GMM scheme
on such structures, we compare the GMM scheme with a standard RWG scheme on
two geometries using multi-scale discretizations. Fig. 4.20 shows the variation of the
condition number with frequency of the impedance matrix generated for a square
PEC plate (side 1m) and a thin almond of aspect ratio 1:3:9, both discretized nonuniformly. The ratio of the largest to smallest edge length used is 1:4 for the plate
and 1:10 for the almond. The condition number is compared with both an RWG and
86

Figure 4.20: EFIE Condition number comparisons for non-uniformly discretized
square plate (1m) and thin almond(1 : 3 : 9) using GMM and RWG and loop star
basis. The plate is discretized at average edge length ranging from 0.1λ to 0.025λ
and the almond is discretized from 0.1λ to 0.02λ.

a loop star scheme as the frequency is reduced. As can be clearly seen the condition
number of the MoM scheme grows exponentially whereas the GMM condition number
varies very slowly with frequency.
Finally, Fig. 4.21 shows the RCS for a square PEC plate (of side 1m) obtained
using the GMM and loop star basis functions. The figure shows the RCS at three
different frequencies (30M Hz, 30kHz and 3kHz). The excellent agreement at all
three frequencies demonstrate the low frequency accuracy of the proposed GMM
scheme.
The variety of examples presented above demonstrate the ability of the GMM
scheme coupled with the quasi-helmholtz basis functions to construct a well conditioned discretization scheme for the EFIE at low frequencies. While the full rationale
87

Figure 4.21: Comparison of RCS (GMM vs. loop star) at 30M Hz, 30kHz and 3kHz
for constant h for a Plate: side 1m. The RCS is compared against that obtained
using Loop Star basis functions.

88

for this behavior is still under investigation; the quasi-Helmtholtz nature of the basis
functions is expected to be the major contributing factor to the well conditioned nature of the scheme. As an illustration of this idea; in figure 4.22 we study condition
numbers for a PEC box of sides 0.5m × 0.5m × 0.1m where we use just one of the
two types of basis functions in each patch (we use the quasi-divergence free basis
function). The RCS obtained by using just the quasi-divergence free basis function
is compared against the RCS obtained using the complete GMM basis set, to ensure
that the currents are represented adequately. The comparison of RCS is shown inset
in the figure.
The condition numbers don’t seem to maintain the near-constant nature that is
seen when the complete basis set is used. This points, indirectly, to the importance
of the role played by the quasi-Helmholtz nature of the basis functions. However,
it should be pointed out that the GMM framework is crucial to the implementation
of a locally quasi-Helmholtz basis function set and therefore; is critical to the well
conditioned nature of the discretized integral equation.
Another key feature of the quasi-surface Helmholtz decomposition implemented
in the GMM framework is that the decomposition is local to each patch. As a result,
this technique is directly applicable to multiply connected domains; as opposed to the
loop star system [93] which uses a global Helmholtz decomposition. A problem with
a global Helmholtz decomposition can be viewed as a result of the Hodge decompositions of the surface fields. This means that the loop start technique is not applicable
to multiply connected domains. In contrast, the GMM scheme is directly applicable
to multiply connected domains such as a toroid.
To illustrate this, figure 4.23 below compares the scattering from a toroidal scatterer computed using the GMM discretization compared to that obtained using a
standard RWG discretization. The incident field direction and polarization is shown
inset in the figure. The excellent comparison of the RCS demonstrates the ability
89

Figure 4.22: Condition number vs. frequency for a PEC box: The figures show
the condition numbers using the regular GMM basis set and using just the quasidivergence free basis functions from the GMM basis set. In each case; the RCS is
computed to good accuracy as shown inset.

90

Figure 4.23: Bistatic RCS (φ = 0) comparison for a dielectric toroid (outer radius
0.75λ0 and inner radius 0.5λ0 , εr = 4.0). This result demonstrates the advantage of
using an local quasi-surface Helmholtz decomposition basis function set. The inset
shows direction and polarization of the incident field along with an image of the
surface current density on the toroid.
of the GMM framework using the quasi-surface Helmholtz basis functions to handle
multiply connected domains.

4.14

Summary

In this work we have presented a framework for the discretization of electromagnetic
integral equations that (i) provides considerable freedom over the choice of basis
functions and (ii) eliminates the traditional coupling of the basis function and the
discretization. Consequently, this provides a general structure wherein different basis
functions can be seamlessly integrated to best approximate the currents locally. More
importantly, it implies that it is not necessary to have basis functions that satisfy
91

specific continuity requirements. As this is a departure from standard basis functions
defined on tessellations, we have devised a methodology to use existing meshes to
define basis functions and their domains of support. Results presented demonstrate
the excellent approximation properties of the method as well as far field RCS data
from three dimensional objects. Several results have been presented to showcase
several of the advantages of the GMM scheme. The results show that the scheme is
low frequency stable, can handle mixture of types and orders of basis functions and
that the mechanism is able to handle non conformal meshes.
In the next two chapters; we will extend the GMM discretization scheme to the
solution of scattering from dielectric bodies and to the analysis of transient scattering
respectively.

92

Chapter 5
Application of the Generalized
Method of Moments for Dielectric
Bodies
So far, we have restricted our study to scattering from bodies that are perfectly
conducting. We will now show how the GMM technique can be extended to the
analysis of penetrable scatterers. We will start from the general integral equations
developed in chapter 2. We restate the problem and equations slightly differently
here, in a form that allows us to specialize to penetrable bodies.

5.1

Problem Statement

Consider a penetrable body that occupies a domain D and is bounded by a surface Ω.
Let D+ and D− denote regions that are exterior and interior to Ω; henceforth, the
subscript q = ± will be used to denote all quantities associated with these regions.
If no subscript is used, then it is assumed that the quantity is associated with region
D+ . Let the surface Ω be sufficiently smooth so that it is possible to define a unique
93

ˆ
ˆ
normal n± (r) that points into regions D± , respectively. It follows that n+ (r) =
ˆ
−ˆ − (r) = n(r) ∀r ∈ Ω. Further, assume that the object is homogeneous, has
n
permittivity and permeability {ε− , µ− }, respectively, and resides in an unbounded
homogeneous background with constitutive parameters {ε+ , µ+ }. Sources residing
within Dq produce incident electromagnetic fields {Ei (r), Hi (r)}, and the interaction
q
q
of these field with Ω produces scattered field {Es (r), Hs (r)}. The total fields in any
q
q
region, {Etot (r), Htot (r)}, comprise of the incident and scattered fields
q
q
Etot (r) = Ei (r) + Es (r)
q
q
q
Htot (r) = Hi (r) + Hs (r)
q
q
q

(5.1)

The scattered fields within region Dq are produced by equivalent sources that reside
on the surface Ω and are defined as
Mq (r) = −ˆ q × Etot (r)
n
q
ˆ
Jq (r) = nq × Htot (r)
q

(5.2)

From these definition and the continuity of the tangential components of the fields, it
follows that J+ (r) = −J− (r) and M+ (r) = −M− (r). These currents can be used
to define magnetic and electric vector potentials

Aq (r) = µq
Fq (r) = εq

Ω
Ω

dr gq (R)Jq (r )

(5.3a)

dr gq (R)Mq (r )

(5.3b)

where κq is the wavenumber in region q, R = |R| = |r−r |, and gq (R) =

exp[−jκq R]
4πR

is homogeneous Green’s function with parameters equal to that of region q. These

94

potentials can in turn be related to the scattered fields via
1
Es (r) Jq , Mq = −jω I +
q
κ2
q
Hs (r) Jq , Mq =
q

1
µq

· Aq (r) −

1
εq

× Aq (r) − jω I +

1
κ2
q

× Fq (r)

· F(r)

(5.4a)

(5.4b)

In these equations, I is used to denote the idempotent. It is well known that knowledge of the equivalent currents is sufficient to express the total fields everywhere via
(5.1) and (5.4). Given four equations (two in each region) that can be formulated
using these currents, different combinations of these have been developed to solve for
the two currents [23]. Of these, the most often used sets are the PMCHWT [26] and
the M¨ller equations [23, 111]. The former is a Fredholm first kind integral equation,
u
whereas the second is a Fredholm second kind equation. In what follows, we will
present a succinct derivation of both these equations for completeness.
The PMCHWT set of equations is derived by imposing the continuity of each field
ˆ
ˆ
ˆ
across Ω. That is, imposing n(r) × Etot (br) = n(r) × Etot (r) and n(r) × Htot (br) =
+
+
−
ˆ
n(r) × Htot (r) results in
−
ˆ
ˆ
n(r) × Einc (r) − Einc (r) = n(r) × Es (r) J− , M− − Es (r) J+ , M+
+
+
−
−
= −ˆ (r) × Es (r) J+ , M+ + Es (r) J+ , M+
n
+
−
= −ˆ (r) × Es (r) + Es (r)
n
+
−

J+ , M +
(5.5a)

ˆ
ˆ
n(r) × Hinc (r) − Hinc (r) = n(r) × Hs (r) J− , M− − Hs (r) J+ , M+
+
+
−
−
= −ˆ (r) × Hs (r) J+ , M+ + Hs (r) J+ , M+
n
+
−
= −ˆ (r) × Hs (r) + Hs (r)
n
+
−

J+ , M+
(5.5b)

95

On the other hand, the M¨ller equations are obtained using a linear combination of the
u
currents in the two regions and the total fields. In other words, a linear combination
of the magnetic currents in the two regions together with the total fields yields
ˆ
− n(r) × Einc (r+ ) + αEinc (r− )
+
−
ˆ
= (1 + α)M+ + n(r) × Es (r+ ) J+ , M+ + αEs (r− ) J− , M−
+
−
ˆ
= (1 + α)M+ + n(r) × Es (r+ ) − αEs (r− )
+
−

(5.6a)

J+ , M+

Similar operations on the electric current yields
ˆ
− n(r) × Hinc (r+ ) + βHinc (r− )
+
−
ˆ
= −(1 + β)J+ + n(r) × Hs (r+ ) J+ , M+ + βHs (r− ) J− , M−
+
−
ˆ
= −(1 + β)J+ + n(r) × Hs (r+ ) − βHs (r− )
+
−

(5.6b)

J+ , M+

where r+ and r− denote points that are just above and just below r for all r ∈ Ω.
In these equation, α and β are constants. In our implementation, we use α = ε− /ε+
and β = µ− /µ+ as suggested by M¨ller. To solve these equations numerically, one
u
typically represents the currents in terms of spatial basis functions. That is, J+ (r) =
Ns J f (r) and M =
−
i=1 i i

Ns M f (r), where J , M are unknown coefficients, f (r)
i i
i
i=1 i i

are basis functions, and Ns is the total number of spatial degrees of freedom for a
current. Using Galerkin testing with the same space of basis functions yields a matrix
equation that may then be solved for the unknown coefficients. Traditionally, RWGbasis functions have been used to represent the spatial variation of the currents and
the fields [89, 112].

96

5.2

GMM Discretization of the PMCHWT and
M¨ ller Integral Equations
u

The PMCHWT formulation involves the discretization of the L± |t and the K± |t
operators (where |t denotes the tangential component of the functional). While the
former can be discretized in a manner identical to that done for conducting bodies in
chapter 4, the K± operator needs to be handled slightly differently. For self triangles,
the integral in the PEC MFIE can be computed for flat patches as as just the principal
ˆ
value component due to the fact that n × ji (r) × g± (R) = 0, when r = r . However
for the PMCHWT, the tangential component of ji (r) ×

5.2.1

g± (R) is 0.

M¨ ller Operator
u

The discretization of the M¨ller operator involves the computation of the following
u
bilinear forms.
.
ˆ
ˆ
ˆ
< tj , n × L± ◦ ji >= tj · n × L± · ji = n × tj · L± · ji

(5.7a)

ˆ
tj · n × K ± · j i

(5.7b)

and

As shown in equations (5.7a) and (5.7b), the bilinear form arising from the K operator
is identical to the case of the PEC in chapter 4 and that arising from the L operator
ˆ
can be effected by testing the operator using n × tj (as opposed to tj . Thus the
ˆ
only modification required to implement the GMM is the evaluation of n × ji and
ˆ
· n × ji . To perform these operations, we revisit our basis functions described in
chapter 4. The basis function set is formed of two functions given by
.
jm (r) = ψi s φm (r)
i
i,d
97

(5.8a)

and
.
ˆ
jm (r) = ψi n ×
i,c

m
s ζi (r).

(5.8b)

ˆ
To compute n× these basis functions, if we assume that φm (r) is chosen to be the
i
m
ˆ
same as ζi (r), we observe that n × ji,d = ji,c . From an algorithmic view point, this
means that if we have a subroutine that computes ji,c , we can recycle that to compute
ˆ
ˆ
ˆ
ˆ
n × ji,d . Also, n × ji,c requires the computation of n × ψi n ×

s · = −ψi s · which

is algorithmically identical to the computation of jm . Thus the changes required
i,d
to implement the GMM for dielectric bodies can be trivially effected using the same
code as used for the PEC case.

5.3

Implementation Results

Next, we will present a series of results that serve to demonstrate several features of
the GMM framework. We begin by presenting results that validate the application of
the GMM framework to the solve both the PMCHWT and M¨ller equations. Here, the
u
data obtained is compared either against analytical data or those obtained using RWG
basis functions. Following this set, we will present results that (i) demonstrate hpconvergence, (ii) show well conditioned behavior for the PMCHWT formulation, (iii)
demonstrate ability to handle non-conformal meshes with relative ease, (iv) ability
to easily mix different order/types of basis function within a single simulation, and
finally (v) ability to handle multiply connected domains. In all the cases considered
herein, we will compare bistatic radar scattering cross-section data. It will be assumed
that the scatterer geometry is discretized using triangular patches whose average edge
length is approximately a tenth of the wavelength in the medium (or 0.1λ− ). Further,
unless otherwise specified, the same mesh is used to obtain results for both GMM and
RWG basis functions. Also, λ0 will be used to denote the wavelength of the incident
field in free space.
98

5.3.1

Validation of GMM PMCHWT and M¨ ller Solvers
u

Next, we validate GMM based PMCHWT and M¨ller integral equation solvers. Unu
less otherwise mentioned, the permittivity and permeability of the dielectric object is
.
assumed to be ε− = εr ε0 = 4.0ε0 and µ− = µ0 , where ε0 and µ0 are the permittivity and permeability of free-space, respectively. Further, it is assumed that all objects
considered herein are immersed in free-space, i.e., ε+ = ε0 and µ+ = µ0 . In each
case, the object under study is discretized using triangular patches of an average edge
length 0.1λ− . GMM patches and basis functions are constructed starting from this
triangulation. The RCS obtained is compared against analytical data (for spheres)
or against a traditional RWG scheme constructed on the original tessellation.
First, Figure 5.1a compares the RCS data for three dielectric spheres obtained
using the GMM and RWG based codes, as well as that obtained using Mie series [26].
The radii of the spheres are 0.1λ0 , 0.3λ0 and 0.5λ0 respectively. The incident field is
ˆ
polarized along x and propagates along −ˆ. The GMM discretization for the spheres
z
result in NGM M = {400, 776, 9416} total (J and M) unknowns, respectively. The
order of basis function p is chosen to be one. The number of unknowns when using
RWG basis are NRW G = {596, 1152, 12004}, respectively. The bistatic RCS (at
φ = 0) obtained from the GMM, RWG based solvers for PMCHWT and M¨ller
u
equations are compared in Figures 5.1a against analytical data. As is evident, the
three sets of data agree very well with each other, demonstrating that the low order
GMM performs as well RWG.
ˆ
Next, we consider an x polarized plane wave, traveling in the −ˆ direction that
z
is incident on a dielectric cube of side length λ0 , where λ0 is the wavelength in free
space. The object was discretized using triangular patches, and the corresponding
number of the GMM and RWG unknowns for this mesh are NGM M = 2404 and
NRW G = 3600. The GMM discretization is applied to both the PMCHWT and the
99

M¨ller equations. As is evident from the RCS data at the φ = 0 plane shown in
u
Figure 5.1b, the two sets of results agree very well with each other.
Next, we consider scattering of a plane wave from two non-canonical geometries–
a NASA almond and a conesphere. The NASA almond fits in a box of dimension
ˆ
0.4λ0 ×1.15λ0 ×3λ0 . The incident field is polarized along x direction, and propagates
along the −ˆ direction. The number of GMM and RWG basis are NGM M = 3056
z
and NRGW = 5600, respectively. Figure 5.2a compares the RCS data obtained in
the φ = 0 plane using PMCHWT and M¨ller integral equations that are solved using
u
GMM and RWG basis functions. It is evident the agreement between the four RCS
data is excellent.
Finally, Figure 5.2b the compares the RCS data obtained for a dielectric cone
sphere of maximum radius 0.5λ0 and cone height 2.6λ0 . As before, data is obtained
using both PMCHWT and M¨ller integral equations, using both GMM and RWG
u
basis functions. In this test, the incident field is polarized along −ˆ , propagates
y
ˆ
along x and bistatic RCS is observed at θ = π/2. The number of total unknowns
is NGM M = 2900 and NRW G = 3312 respectively. As can be seen in each of the
above cases, the agreement between the RCS is excellent. All these results provide a
measure of confidence on the GMM

5.3.2

Salient features of the GMM scheme

Thus far, we have demonstrated that GMM based discretization of both the PMCHWT and M¨ller equation produces accurate results. In what follows, we shall
u
present data that highlight some of the features of is approach, and contrast it with
what can be done using RWG basis functions.
100

(a) Spheres
Figure 5.1: RCS obtained using GMM basis functions used to discretize the PMCHWT amd M¨ller equations for the solution of scattering from (a) dielectric spheres
u
(radii 0.1λ0 ,0.3λ0 and 0.5λ0 and εr = 4.0) and (b) a cube (side λ0 , εr = 4.0) due
to plane wave incidence. The bistatic RCS (at φ = 0) are compared against that
obtained using an RWG discretization and analytical results (for spheres).

101

(b) Cube
Figure 5.1 Continued: RCS obtained using GMM basis functions used to discretize
the PMCHWT amd M¨ller equations for the solution of scattering from (a) dielectric
u
spheres (radii 0.1λ0 ,0.3λ0 and 0.5λ0 and εr = 4.0) and (b) a cube (side λ0 , εr = 4.0)
due to plane wave incidence. The bistatic RCS (at φ = 0) are compared against that
obtained using an RWG discretization and analytical results (for spheres).

102

(a) Almond (0.4λ0 × 1.2λ0 × 3.0λ0 )
Figure 5.2: RCS obtained using GMM basis functions due to scattering from (a)
NASA thin almond of εr = 4.0 that fits in a box of size 0.4λ0 × 1.2λ0 × 3.0λ0 due to
ˆ
a plane wave incident along −ˆ and polarized along x and (b) a cone-sphere of radius
z
0.5λ0 , cone height 2.6λ0 and εr = 4.0, due to a plane wave polarized along −ˆ and
y
ˆ
incident along x. Bistatic RCS is computed at φ = 0 for the almond and at θ = π/2
for the conesphere and used to validate the GMM scheme against RWG.

103

(b) Cone-sphere (0.5λ0 , 2.6λ0 )
Figure 5.2 Continued: RCS obtained using GMM basis functions due to scattering
from (a) NASA thin almond of εr = 4.0 that fits in a box of size 0.4λ0 ×1.2λ0 ×3.0λ0
ˆ
due to a plane wave incident along −ˆ and polarized along x and (b) a cone-sphere
z
of radius 0.5λ0 , cone height 2.6λ0 and εr = 4.0, due to a plane wave polarized along
ˆ
−ˆ and incident along x. Bistatic RCS is computed at φ = 0 for the almond and at
y
θ = π/2 for the conesphere and used to validate the GMM scheme against RWG.

104

Mixtures of Basis functions
As was alluded to earlier, the GMM framework permits mixtures of basis functions
within the simulation domain. This implies that one can use a set of basis function
in a sub-domain and use an entirely different set in neighboring sub-domains without
invoking any special condition regarding continuity of different orders. This is possible
as patches overlap with each other and continuity is built into the basis functions.
This property is extremely useful when analyzing complex structures where h-, p- or
hp-adaptivity is required or if known asymptotic solutions are to be built into the
approximation space. In what follows, we shall explore these features of the GMM
framework. Note, it should be emphasized that the use of different basis functions
with different patches is as simple as choosing a flag associated with a node in the
input file that describes geometry. This is so as each node is associated with a patch
and this flag can describe the function that is to be used in the patch.

Mixed Order Bases: The first class of numerical tests will be carried out using
mixtures of basis functions of different orders. Figure 5.3 compares the RCS obtained
due to scattering from the almond and box using two different orders of GMM basis
ˆ
functions. In both cases, the incident wave is x polarized and travels in the −ˆ
z
direction. The relative permittivity of both scatterers is εr = 4, and the RCS data is
obtained in the φ = 0 plane.
Figure 5.3a compares the RCS data for an almond obtained using p = 1 everywhere on the scatterer with that obtained using p = 1 or p = 2 in different regions.
Specifically, basis functions on all patches that fall within 0.2λ− of the tip were of
order p = 2 and all others were of order p = 1. As is evident from Figure 5.3a, the
two sets of data agree very well with each other.
Next, Figure 5.3b compares the RCS data obtained for scattering from a thin
dielectric box of size 0.2λ0 × 0.2λ0 × 4λ0 with varying h and p with a reference
105

(a) Almond (0.4λ0 × 1.2λ0 × 3.0λ0 )
Figure 5.3: RCS comparison of standard GMM vs. Mixed-order GMM (a) NASA
almond discretized with GMM basis functions complete to order p = 4 near edges
and p = 1 elsewhere. (b) Dielectric box using p = 1 near edges and p = 2 in larger
patches near center. The bistatic RCS is computed at φ = 0 and the incident field is
ˆ
polarized along x and incident along −ˆ.
z
data, obtained using GMM basis functions in a dense mesh with NGM M = 1552
unknowns, and p = 1. As illustrated in the inset, the box was discretized using the
following rule: h = 0.1λ− within h = 0.25λ− of any edge and 0.2λ− elsewhere. For
patches with h = 0.1λ− the order of basis function was p = 1, and for h = 0.25λ−
the order was p = 2. This resulted in a system with NGM M = 440 unknowns. As is
evident from Figure 5.3b that the two sets of results agree very well with each other.

Mixed Types of Bases: The next example demonstrates the use of a mixture of
types of basis functions. Figure 5.4 compares the RCS obtained from a dielectric
box of dimensions 5.0λ0 × 2.0λ0 × 0.5λ0 and εr = 2.0. This box is first meshed
106

(b) Box 0.2λ0 × 0.2λ0 × 4λ0
Figure 5.3 Continued: RCS comparison of standard GMM vs. Mixed-order GMM (a)
NASA almond discretized with GMM basis functions complete to order p = 4 near
edges and p = 1 elsewhere. (b) Dielectric box using p = 1 near edges and p = 2 in
larger patches near center. The bistatic RCS is computed at φ = 0 and the incident
ˆ
field is polarized along x and incident along −ˆ.
z

107

uniformly using a triangulation of average edge length 0.1λ− . GMM basis functions
are constructed on this tessellation and the RCS is computed by solving the M¨ller
u
integral equation; using NGM M = 21680 unknowns. The box is then re-meshed
using two sets of elements, of average edge length 0.1λ− near the edges and 1.0λ−
at the center; as shown in the inset figure. GMM basis functions are constructed
using polynomials as approximation functions in the smaller patches, using the quasisurface Helmholtz technique described in Section 4.3 and plane wave approximation
.
√
ˆ
ˆ
functions of the form ν m (r) = u cos(θm ) exp −jω µ− ε− cos(θm )k · (r − ri ); where
i
ˆ
ˆ
u is a pilot vector along the surface of the patch, k is the incident field direction, ri is
the patch center and θm is an order parameter between 0 and 2π; in the larger patches.
These functions are used to discretize the M¨ller equations with NGM M = 12350
u
unknowns; and the RCS’s are compared.
The RCS obtained using the regular GMM and the mixed order GMM is demonstrated in Figure 5.4. It is evident from the figure that the RCS agrees very well.
This demonstrates the efficacy of the GMM scheme in mixing different types of basis
functions; a capability that is not available in any of the basis function frameworks
currently in use for the discretization of electromagnetic integral equations.

hp-adaptivity: In this section, we utilize the flexibility of the GMM framework to
study the hp− convergence of the RCS due to scattering from a dielectric arrow. The
arrow has dimensions 2.0λ0 × 3.0λ0 × λ0 and εr = 2.0. RCS data is obtained for the
arrow discretized at λ− /40 to provide a baseline for comparison. The RCS data is
then computed using the GMM basis functions constructed on this tessellation with
NGM M = 64, 048 unknowns. The incident field is polarized along −ˆ and incident
z
ˆ
along x. This RCS is then compared against the following different discretizations
in Figure 5.5. First, the arrow is meshed using elements with an average edge-length
of λ− /10 and approximation functions complete to order p = 1 are used everywhere
108

Figure 5.4: RCS comparison of standard GMM vs. Mixed-Basis GMM (5.0λ0 ×
ˆ
2.0λ0 × 0.5λ0 ). Electric field is incident along −ˆ and polarized along x. The bistatic
z
RCS (φ = 0) is compared for GMM using polynomial basis functions against GMM
using polynomials mixed with plane wave basis functions.

109

Figure 5.5: RCS comparison for scattering from a dielectric arrow (2.0λ0 ×3.0λ0 ×λ0 ).
ˆ
Electric field is incident along −ˆ and polarized along x. The bistatic RCS (φ = 0) is
z
compared for four different combinations of h and polynomial order of basis function
p.

(NGM M = 8960). In the second test, the mesh is refined to an average edge length
of λ− /20 near the tips and basis functions of order p = 1 are applied near the
tip and basis functions of order p = 2 everywhere. This case is shown as hp − 1
(NGM M = 12, 048) in Figure 5.5. Finally, the mesh is refined further, in two levels,
with λ− /20 at the tip, λ− /10 in the smooth areas and λ− /15 in the transition region.
In this case, basis functions of order p = 1 are used in the patches of size λ− /10 and
p = 2 everywhere else. This case is shown as hp − 2 (NGM M = 22, 860). As can be
seen the RCS demonstrates the hp-convergence of the method.
110

Non Conformal meshes
One advantage of disassociating the definition of the basis function from the underlying tessellation of the geometry is that it allows for the use of non conformal meshes.
It is common, in the case of practical geometries, to have different vendors mesh different parts of a geometry. This results in sets of meshes that need to be explicitly
“stitched together” at the boundaries. This motivates the design of techniques that
can handle non-conformal meshes [110]. However, a consequence of the design of
the GMM scheme is that it can implicitly handle these kinds of meshes without any
change to the framework or the implementation.
Here, we present one example illustrating this ability. We consider scattering from
a dielectric box (εr = 4.0) of dimensions 2.0λ0 × 0.2λ0 × 0.5λ0 . Two halves of the
box are meshed at slightly different discretization rates to obtain a non conformal
mesh. The coarser mesh is 0.125λ− while the finer is 0.1λ− . The RCS is computed
ˆ
ˆ
due to scattering of a x polarized plane wave incident along z and is observed along
φ = 0. The RCS obtained due to scattering from a discretization constructed on
this non conformal mesh is compared against that is obtained by constructing a fully
conformal mesh on the geometry. The number of unknowns is NGM M = 5600 for
the non conformal mesh and NGM M = 6004 for the conformal mesh. The agreement
of the RCS data, shown in Figure 5.6, demonstrates the ability of the GMM scheme
to seamlessly handle non conformal meshes.

5.3.3

Low frequency stability and multiply connected scatterers

It is well known that the condition number of the discrete integral equation is intimately related to the discretization size. A rule of thumb for linear basis functions
is that the characteristic dimension h vary approximately as λmin /10. This means
111

(a) Non conformal mesh for a box

(b) RCS comparison
Figure 5.6: RCS comparison for scattering from a non conformally meshed dielectric
box. Bistatic RCS (φ = 0) obtained by discretizing the M¨ller integral equation on a
u
non conformal mesh (shown inset) compared against that obtained by discretizing on
a (standard) conformal mesh for a thin dielectric (εr = 4.0) box (2.0λ0 ×0.2λ0 ×0.5λ0 )

112

that the size of discretization is governed by both the frequency and the dielectric
constant. Further, the size of discretization is often governed by the need to accurately capture geometric features and, as a result, may be quite small. Next, the
two equations that are used in this paper, the PMCHWT and the M¨ller fall under
u
the category of Fredholm first and second kind, respectively. It can be proven that
the eigenvalue of the former spread further apart as the frequency becomes higher,
with the smallest accumulating around zero and the largest tending to infinity. In
contrast, the eigenvalues of the M¨ller cluster around a non-zero value [113]. As a
u
result, discrete system that arise from the PMCHWT tend to be ill-conditioned when
the wavelength is large compared to the feature size. In other words, they suffer from
low-frequency breakdown. The common remedy to solve this problem has been the
development of loop-star or loop-tree basis functions, that separate the basis functions
into an solenoidal and non-solenoidal components [93]. This decomposition confers
low-frequency stability to the discrete PMCHWT equations, however, unless these
decompositions are augmented, it is not possible to use them for analysis of multiply
connected objects. As the GMM basis functions are neither purely solenoidal nor
non-solenoidal, and the basis functions are purely local to a patch. They have been
shown excellent approximation properties of the function, its divergence and curl. In
what follows, we will demonstrate that they yield well conditioned systems that can
be used for analysis of multiply connected objects.
We compare the condition number of systems resulting from the discretization of
PMCHWT equations using GMM and RWG basis for three different scatterers; (i) a
sphere of radius 0.1m, (ii) a cube of side 1m and (iii) a thin NASA almond that fits in
a box of size 0.4m × 1.15m × 3m. The relative permittivity of each body is εr = 4.0.
Figure 5.8 demonstrates the condition number of the impedance matrix obtained
using the GMM basis functions and using RWG functions over a range of frequencies,
starting from 30Hz to 30MHz. In each case, the object is discretized at 0.1λ− at
113

the highest operating frequency. The number of GMM and RWG unknowns are
NGM M = {400, 2404, 3056} and NRW G = {596, 7200, 5600} for the sphere, cube
and almond, respectively. As is immediately apparent from Figure 5.8, the condition
number of systems discretized using RWG basis grows exponentially as the frequency
decreases whereas those discretized using GMM are relatively constant in across a six
decade change in frequency.

(a) Sphere (0.1m)
Figure 5.7: Condition number comparisons. All objects are discretized at 0.1λm at
30M Hz and the operating frequency is swept from 30M Hz to 30Hz and condition
number of the impedance matrix obtained on GMM discretization is compared against
that obtained using a standard RWG discretizaion. (a) Sphere (radius 0.1m, εr = 4.0)
and (b) Cube (side 1m and εr = 4.0).
Next, we illustrate that these basis functions can be used to analyze scattering
from multiply connected objects unlike loop-star or loop-tree functions [96]. Figure
5.9 below compares the scattering from a toroidal scatterer of outer radius 0.75λ0 ,
114

(b) Cube (1m)
Figure 5.7 Continued: Condition number comparisons. All objects are discretized
at 0.1λm at 30M Hz and the operating frequency is swept from 30M Hz to 30Hz
and condition number of the impedance matrix obtained on GMM discretization is
compared against that obtained using a standard RWG discretizaion. (a) Sphere
(radius 0.1m, εr = 4.0) and (b) Cube (side 1m and εr = 4.0).

115

Figure 5.8: Condition number comparisons. All objects are discretized at 0.1λm at
30M Hz and the operating frequency is swept from 30M Hz to 30Hz and condition
number of the impedance matrix obtained on GMM discretization is compared against
that obtained using a standard RWG discretizaion. (c) NASA thin almond (side
0.4m × 1.2m × 3.0m and εr = 4.0)
inner radius 0.5λ0 and of relative permittivity 4.0 computed using the GMM discretization of the M¨ller equation (NGM M = 3406) compared to that obtained
u
using a standard RWG discretization (NRW G = 7800). The incident field direction
(−ˆ) and polarization (ˆ ) is shown inset in the figure. The excellent comparison of
z
x
the RCS demonstrates the ability of the GMM framework using the quasi-surface
Helmholtz basis functions to handle multiply connected domains.

116

(a) Surface current on dielectric toroid

(b) RCS comparison for dieletric toroid
Figure 5.9: Bistatic RCS (φ = 0) comparison for a dielectric toroid (outer radius
0.75λ0 and inner radius 0.5λ0 , εr = 4.0). Figure 5.9(a) shows direction and polarization of the incident field along with an image of the surface current density on the
toroid.

117

Chapter 6
Analysis of transient scattering
using the Generalized Method of
Moments
In this chapter, we extend the GMM technique to the solution of transient scattering
from perfectly conducting bodies. We will proceed by deriving the time domain EFIE
and MFIE. While these equations can be derived from first principles, in this work
we will restrict to a simpler derivation starting from the frequency domain equations
using the Fourier transforms.

6.1

The Fourier transform pair

The temporal Fourier transform pair that we will use for the derivation of the TDEFIE and TD-MFIE are given by
.
F(r, ω) = F{F(r, t)} =
118

∞
−∞

dt F(r, t)e−jωt

(6.1a)

and
. 1 ∞
F(r, t) = F−1 {F(r, ω)} =
dt F(r, ω)ejωt .
2π −∞

(6.1b)

We will freely use the following identities without proof:
∂
F(r, t)
∂t

F

= jωF(r, ω),

(6.2a)

F {δ(t − τ )} = e−jωτ ,

(6.2b)

which directly lead to
F−1

e−jkR
4πR

= F−1

e−jωR/c
4πR

=

δ(t − R )
c
4πR

(6.2c)

and
F−1 {F(r, ω)G(r, ω)} = F(r, t) ∗t G(r, t),

(6.2d)

where ∗t denotes convolution in time defined as
∞
F(r, t) ∗t G(r, t) =

6.2

−∞

dτ F(r, τ )G(r, t − τ ).

(6.2e)

Time Domain EFIE and MFIE

For the sake of completeness, we restate our scattering problem in a transient framework here. Consider a perfectly electric conducting body occupying a domain D− ⊂
R3 in free space (µ0 , ε0 ). Let Ω denote the bounding surface of D− and let the outˆ
ward pointing unit vector normal at any point r on Ω be denoted by n(r). Assume
that a field {E i (r, t), Hi (r, t)} is incident on this body. This field induces currents,
J (r, t), on the surface Ω, that radiate scattered fields {E s (r, t), Hs (r, t)}. Thus, the
total electric and magnetic fields in all of D+ = R3 /D− , denoted by E(r, t) and
119

H(r, t), can be decomposed into the incident and scattered fields by
E(r, t) : R3 × R+ → R3 = E i (r, t) + E s (r, t)

(6.3a)

H(r, t) : R3 × R+ → R3 = Hi (r, t) + Hs (r, t).

(6.3b)

and

{E i (r, t), Hi (r, t)}
ˆ
ns
D

{E s (r, t), Hs (r, t)}

−

Ω

Figure 6.1: General description of a scattering problem

Applying the Fourier transform to equations (2.11) and (2.12), we obtain their
time domain equivalents as
Hs (r, t) =

1
µ0

× A(r, t)

(6.4a)

and
E s (r, t) + ∂t A(r, t) = − φe (r, t)

(6.4b)

Coupling (6.4a) and (6.4b) with the Lorentz gauge given by

· A(r, t) + µ0 ε0
120

∂φe (r, t)
=0
∂t

(6.5)

, we obtain
E s (r, t) = −∂t A(r, t) +

t
1
du
µ0 ε0 −∞

· A(r, u).

(6.6)

This provides us a relationship between the magnetic vector potential and the fields.
The magnetic vector potential can be related to the induced current J (r, t) by the
Fourier transform of equation (2.27) as

A(r, t) = µ0

Ω

dr J (r , t) ∗t

δ(t − τ )
.
4πR

(6.7)

where τ = t − R/c and R = |R| = |r − r | and using (6.7) in (6.6) and (6.4a). We
assume that all the integrals contain finite energy and that we can switch the order
of integration to obtain the following equations for the scattered fields generated by
J (r, t):
E s (r, t) = −µ0 ∂t

dr
Ω

1
J (r , τ )
+
dr
4πR
ε0 Ω

τ
·

−∞

dt

J (r , t)
4πR

(6.8a)

and
Hs (r, t) =

×

dr
Ω

J (r , τ )
.
4πR

(6.8b)

Equation (6.8b) can be, after some manipulation, reduced to
Hs (r, t) =

1
dr
4π Ω

1
1
∂t J (r , τ ) +
J (r , τ ) × R.
cR
R2

(6.8c)

Applying the PEC boundary conditions from chapter 2 and proceeding very similar to
the derivation of the frequency domain EFIE and MFIE, we obtain the time domain
versions of the EFIE and the MFIE as
ˆ ˆ
ˆ ˆ
n × n ×E i (r, t) = n × n ×

µ0
1
∂t dr J (r , τ ) −
dr
4π
4πε0 Ω
Ω

τ
·

−∞

dt J (r , t)
(6.9a)

121

and
ˆ
n × Hi (r, t) = −ˆ ×
n

1
dr
4π Ω

1
1
∂t J (r , τ ) +
J (r , τ ) × R
cR
R2

(6.9b)

Finally, to avoid the integration involved in (6.9a), we use the time derivative of
these equations to obtain the time domain EFIE (TD-EFIE) the time domain MFIE
(TD-MFIE) as
ˆ ˆ
ˆ ˆ
n × n × ∂t E i (r, t) = n × n ×

µ0
1
2
dr ∂t J (r , τ ) −
dr
4π Ω
4πε0 Ω

· J (r , τ )
(6.10a)

and
ˆ
n × ∂t Hi (r, t) = −ˆ ×
n

1
dr
4π Ω

1 2
1
∂t J (r , τ ) +
∂t J (r , τ ) × R
cR
R2

(6.10b)

The standard marching on in time (MOT) [92, 114–116] solution to the above
system of equations expands the unknown current in a set of spatio-temporal basis
Nt
functions as J(r, t) = Ns
n=1 j=1 In,j jn (r)Tj (t). The spatial basis functions jn (r)
typically chosen are RWG functions. However, the temporal basis functions Tj (t)
need to be carefully chosen. In the next section, we will visit this topic in some
detail.

6.3

The Marching on in Time (MOT) scheme

To solve the TD-EFIE and TD-MFIE above, we will proceed similar to the frequency
domain case. The spatio-temporal current J (r, t) is discretized using a set of basis
functions that are assumed to be a product of Ns spatial vector basis functions
(eventually the GMM basis functions) and Nt scalar shifted temporal functions (of
122

.
the form Tj (t) = T (t − jδt), that partition the time into segments of size ∆t.
Ns Nt
J (r, t) ≈

Ij,n T (t − j∆t)jn (r).

(6.11)

i=1 j=1
Substituting (6.11) into the TD-EFIE/TD-MFIE and performing point testing (δ(t −
i∆t)) in time and Galerkin testing (jm (r)) in space results in a set of Nt , order-Ns
matrix equations of the form
Nt
Zi−j Ij = Vi ,

(6.12)

j=1
where
Ij = Ij,n ,






Vi =





Ω
Ω

(6.13a)

ˆ ˆ
dr n × n × ∂t E i (r, i∆t) · jm (r)
(6.13b)
ˆ
dr n × ∂t Hi (r, i∆t) · jm (r)

and


 dr jn (r) · n × n × µ0 dr jn (r )∂ 2 T ((i − j)∆t)
ˆ ˆ


t
 Ω
4π Ω




1


dr
· jn (r )T ((i − j)∆n)
−

4πε0 Ω
Zi−j =
1

2
 − dr jn (r) · n × 1
ˆ
dr
jn (r )∂t T ((i − j)∆t)



4π Ω
cR
Ω




1


+ jn (r )∂t T ((i − j)∆t) × R

R2

(6.13c)

The most critical observation about the (6.13c) is that, if the temporal basis functions
have finite support, that is ∃ k, l · −
⊃
−· T (t) = 0 ∀ t ≤ −(l + 1)∆t, t ≥ (k + 1)∆t,
then if Rmax is the maximum spatial extent of the scatterer (Rmax ≈ diam(Ω)), we

123

have that
Rmax
Zi−j = 0 ∀ (i − j) < −l and (i − j) > k +
.
c∆t

(6.14)

Further, if we assume that the temporal basis functions are causal, that is T (τ ) =
0∀τ ≤ −∆t, we have that Zi−j = 0 ∀ i < j − 1. This is called the “marching
criterion”, in that if this criterion is satisfied, the matrix equation (6.12) can be
re-written as a marching scheme, given by
i−1
Z0 Ii = Vi −

Zi−j I j

(6.15a)

j=i−km
where km is given by
km =

Rmax
+k
c∆t

(6.15b)

This naturally leads to a marching solution of the TDIE whereby the unknowns solved
for in (a finite number of) prior time steps leads to the solution at the current time
step.

6.4

Temporal basis functions

Given the marching scheme above, it is desirable to choose basis functions Tj (t) that
(i) have short temporal support, (ii) do not look forward in time, and (iii) provide
good interpolation properties.
The most commonly used function in the literature are piecewise continuous functions of order p given by

T (t) = fl (t)gp−l (t) for − ∆t + l∆t ≤ t ≤ l∆t ; l = 0 . . . p
124

(6.16a)

where

fl (t) =




1








l=0
k

i=1

and

t − i∆t
−i∆t
p−l

gp−l (t) =
i=1

(6.16b)
l=0

t + i∆t
i∆t

(6.16c)

Clearly, these functions satisfy (i)-(iii) and have been in active use for a very long
time [52, 92, 114]. A plot of these functions are shown in figure 6.2a below. Another
temporal basis function that has been used is the Approximate Prolate Spheroidal
Wave Function (APSWF) given by

sin a
T (t) =

sin(sω0 t)
sω0 t

sinh(a)

2
t
N ∆t − 1
2
t
N ∆t − 1

(6.17)

where ω0 is the highest end of the frequency band of interest, s is an oversampling
factor, N ∈ Z is the APSWF width parameter and a = πN (s−1)/s. These functions,
shown in figure 6.2b can be shown to be band limited (to ωmax = (2s − 1)ω0 ), but
are non causal, (the function value at a time instant depends on future values of
time) [117]. This problem of the APSWF is solved by an extrapolation procedure that
maps the forward part of the function using a least-squares approach, albeit at the
expense of band-limitedness [116].

Either of these two approaches has an interesting

consequence. Upon examination of equation (6.13c), it is clear that the solution of
either TDIE requires the computation of at least one derivative of the temporal basis
function. However, it can be shown that neither the Lagrange polynomials or the
extrapolated APSWF functions used support more than one derivative. This has
been reported to result in inaccuracies in the solution and instabilities in the marching
125

scheme. Some recent efforts to solve this effort attempt to use separable polynomial
approximations to the Green’s function to construct solutions [118, 119].

6.5

GMM for transient scattering

Traditionally the spatial basis functions used for the discretization of the TDIE are
the RWG functions. These functions have been carried over to time domain integral
equations (TDIE) from their frequency domain counterparts and have been shown to
be highly effective in a variety of problems. However, the drawbacks that are suffered
by the RWG functions in frequency domain are prevalent in the analysis of TDIE

(a) Lagrange Polynomial Basis Functions
Figure 6.2: Temporal basis functions for solution of the TDIE. (a) shows the Lagrange
polynomial basis functions which are not band limited as opposed to (b) approximate
prolate spheroidal functions that are band limited but require extrapolation to maintain causality.
126

(b) APS Functions
Figure 6.2 Continued: Temporal basis functions for solution of the TDIE. (a) shows
the Lagrange polynomial basis functions which are not band limited as opposed to
(b) approximate prolate spheroidal functions that are band limited but require extrapolation to maintain causality.

127

as well. While the close marriage between the basis function construction and the
tessellation helps to satisfies continuity constraints across the interior edges of the tessellation, this nature of construction is inherently restrictive. As has been discussed
before, even the use of mixed orders of polynomial basis functions is significantly
involved and usually requires the careful design of patches and “transition” areas. In
the frequency domain, the GMM was designed as an umbrella framework to include
arbitrary functions in the basis space. We also demonstrated that the basis functions
proposed result in a matrix system with stable condition numbers to very low frequencies. Here we extend the GMM scheme to the discretization of the TD-MFIE. We will
describe a scheme for the extension of the GMM to time domain using products of
GMM functions (for spatial basis functions) and interpolatory polynomials (for temporal basis functions). In the remainder of this chapter, we will first show that these
basis functions provide accurate scattering results over a wide variety of geometries.
Further, we will demonstrate several of the advantages of the GMM scheme including
(i) its ability to arbitrarily mix basis functions across the geometry: using different
basis functions in different areas of the tessellation can be trivially achieved under
the GMM framework by simply changing an input flag. We will present results that
show that the low frequency nature of the GMM discretized TD-MFIE is much more
stable than the standard RWG-TD-MFIE.

The actual implementation of the GMM to the MOT scheme described above
is straightforward. We retain the Lagrange polynomial interpolation in time and
the point matching temporal testing scheme. The spatial basis functions can be
replaced with the GMM functions exactly as done in the frequency domain case.The
spatial integrals are now evaluated with a “Green’s kernel” of 1/R as opposed to the
oscillatory Green’s kernel used for the frequency domain implementation.
128

6.6

Numerical Experiments

This section will illustrate a variety of numerical experiments that illustrate the
validity and utility of the GMM discretization scheme. We will assume that the
incident electric field is a Gaussian modulated plane wave of the form Ei (r, t) =
ˆ
ˆ
ˆ
pcos(2πf0 τ ) exp −(τ − tp )2 /2σ 2 , where σ = 6/2πfbw and τ = t − k · r/c. p and
ˆ
k denote the polarization and the direction of propagation of the plane wave respectively. The center frequency f0 and bandwidth fbw describe the frequency characteristics of the signal. Examination of the incident field shows that the energy content
is almost −160dB down at f0 + fω in comparison to the energy at f0 . The scattered
far field is computed from this signal and its Fourier transform is computed. The
backscatter RCS is then extracted as the ratio of the frequency spectrum of the far
field to that of the incident field. This allows us to evaluate the GMM discretization
of the TD-MFIE against both other time domain methods (such as RWG-MOT) and
frequency domain methods (by comparing RCS). The GMM basis functions are constructed starting from a triangular spatial discretization with an average edgelength
of 0.1λm , where λm is the lowest wavelength in the frequency bandwidth.

6.6.1

Validation

First we validate the GMM discretization scheme against both experiments and
against the traditional RWG discretization. In this section, we will provide a plethora
of results comparing the GMM-discretized TDIE compared with the standard RWG
based TDIE. Figure 6.3a shows the backscatter RCS obtained from two PEC spheres,
of radii 0.2m and 1.0m. The GMM discretizations result in NGM M = 388 and
NGM M = 200 unknowns respectively. The RCS is compared to the that obtained
from an MOT solver using RWG basis functions, resulting in NRW G = 576 and
NRW G = 192 unknowns respectively. The center frequency f0 = 30M Hz and band129

width is fbw = 15M Hz. As is clear from the figures, there is excellent match between
the backscatter RCS results.
Figure 6.3b compares the backscatter RCS for a NASA almond that fits in a box
ˆ
of dimensions 0.3m × 1.15m × 3m. The direction of propagation k = −ˆ and the
z
ˆ
ˆ
polarization p = x of the incident wave is shown inset in the figure, along with a
time-snapshot of the surface current. Again f0 = 1e8 and the bandwidth extends
to a low frequency end of 100Hz. The large bandwidth allows for the examination
of the low frequency behavior of the GMM scheme. The number of unknowns are
NGM M = 1528 and NRW G = 1656 respectively. The agreement between the RCS
shown provides further confidence in the GMM discretization scheme.
Next, figure 6.4a compares the backscatter RCS obtained from a thin PEC box of
size 20m × 1m × 20m using the GMM discretization NGM M = 1780 to that using
ˆ
standard RWG NRW G = 2640. The incident electric field is polarized along y and
ˆ
incident along x, on the thin edge of the box. The center frequency of the incident
field f0 = 15M Hz and extends at the low frequency end to 100Hz. The agreement of
the RCS in figure 6.4a further validates the GMM scheme. Finally, in figure 6.4b the
backscatter RCS obtained from a PEC cube of side 0.7m, as a function of the incident
wavelength, is compared against that obtained using both the RWG basis functions
ˆ
and experimental data reported in [109]. The direction of propagation k = −ˆ and
y
ˆ ˆ
polarization p = x of the incident field is again reported inset in the figure along with
a snapshot of the magnitude of surface current on the cube. Again, the agreement
of the RCS to experiment and the RWG scheme validates the GMM scheme for the
discretization of time domain integral equations.

6.6.2

Low frequency stability

The well conditioned nature of the GMM scheme at low frequencies is expected to
result in low frequency stability in time. As a first step to investigating this property,
130

(a) PEC Sphere
Figure 6.3: Backscatter RCS comparison for (a) two spheres of radii 0.2m and 0.5m
and (b) a thin NASA almond (0.3m×1.15m×3m) compared against RWG discretization. The inset figure shows snapshots of the magnitude of the current density on the
almond.

we observe the backscatter RCS of the GMM scheme for very wideband incident fields.
The same scattering problem is solved using both the GMM and RWG discretizations
and the RCS observed in both cases is compared to a frequency domain scheme (using
RWG basis functions) at a number of frequency points all the way down to the end of
the frequency band of the input signal. Figure 6.5a compares the backscatter RCS for
the thin box considered above (20m × 1m × 20m, NGM M = 1780, NRW G = 2640).
The incident field direction, polarization, center frequency and bandwidth are the
same as the example in figure 6.4a above. As is clear from the figure, the RCS
computed using the GMM discretization agrees with the frequency domain RCS down
to very low frequencies. However, the RCS computed using the RWG scheme departs
131

(b) NASA Almond
Figure 6.3 Continued: Backscatter RCS comparison for (a) two spheres of radii 0.2m
and 0.5m and (b) a thin NASA almond (0.3m × 1.15m × 3m) compared against RWG
discretization. The inset figure shows snapshots of the magnitude of the current
density on the almond.

132

(a) Comparison against RWG
Figure 6.4: Backscatter RCS comparison for (a) a thin box of size 20m × 1m × 20m
compared against RWG discretization and (b) a PEC cube of side 0.7m compared
against both RWG and experimental result. The inset figures show snapshots of the
magnitude of the current density on the scatterer.

from the frequency domain RCS at lower frequencies. While this needs to be further
investigated for the EFIE, this nevertheless does point to the low frequency stability
of the GMM discretization.
In figure 6.5b, we illustrate the utility of using a local surface Helmholtz decomposition. The use of a global Helmholtz decomposition, such as in the loop-star
scheme [93], would preclude the application to multiply connected domains, such
as the toroid shown inset in figure 6.5b. However, being based on a local surface
Helmholtz decomposition, the GMM scheme is directly applicable to multiply connected domains. Figure 6.5b compares the backscatter RCS due to the scattering
ˆ
of a −ˆ directed plane wave polarized along x from a PEC toroid of outer radius
z
133

(b) Comparison against experiment
Figure 6.4 Continued: Backscatter RCS comparison for (a) a thin box of size
20m × 1m × 20m compared against RWG discretization and (b) a PEC cube of
side 0.7m compared against both RWG and experimental result. The inset figures
show snapshots of the magnitude of the current density on the scatterer.

134

0.75m and inner radius 0.5m, discretized using 1503 GMM and 3900 RWG unknowns
respectively. The center frequency of the incident field f0 = 300M Hz and the bandwidth is fbw = 150M Hz. The agreement of the RCS demonstrates the validity of
using the GMM scheme for multiply connected domains.

(a) Thin box at very low frequencies
Figure 6.5: Backscatter RCS comparison for (a) thin box (20m × 1m × 20m) and
(b) a toroid (outer radius 0.75m and inner radius 0.5m). The inset figures show a
snapshot of the magnitude of the current density and direction of propagation and
incidence.

6.6.3

Mixtures of basis functions

Here, we present examples to illustrate the capability of the GMM scheme to mix
basis functions in time domain. To this end, we first use the GMM discretization
to generate the backscatter RCS due to scattering from a PEC box of dimensions
ˆ
1.0m × 4.0m × 1.0m aligned with the long axis along x, discretized uniformly at
135

(b) Multiply connected domains
Figure 6.5 Continued: Backscatter RCS comparison for (a) thin box (20m×1m×20m)
and (b) a toroid (outer radius 0.75m and inner radius 0.5m). The inset figures show
a snapshot of the magnitude of the current density and direction of propagation and
incidence.

136

0.2λm at the highest frequency. Following this, the box is re-meshed with an average
edge length of 0.2λm , 0.2m from the two narrow ends, and 0.4λm elsewhere. Then
the GMM discretization is applied to the new mesh using basis functions of order
p = 1 near the edge and p = 2 elsewhere. Figure 6.6a shows the excellent agreement
ˆ
between the RCS obtained from these two cases. The x polarized incident field is
incident along −ˆ and f0 = 60M Hz and fbw = 30M Hz.
z
ˆ
In figure 6.6b the backscatter RCS due to scattering of a x polarized Gaussian
modulated plane wave incident along −ˆ with center frequency 300M Hz and bandz
width 299.9M Hz on the above thin NASA almond (0.3m×1.15m×3m) aligned along
ˆ
the x axis. The almond is discretized with an average edge length of 0.1m. First the
RCS is computed using GMM basis functions of order p = 1 everywhere. Following
this, GMM basis functions complete to order p = 2 are used in all patches within
0.2m of the tip. Again, the RCS match each other very well.
Finally in figure 6.6c we mix three orders of basis functions. The inset in figure
ˆ
6.6c shows a PEC box of dimension 3m × 3m × 6m aligned along the x axis. An
ˆ
x-polarized incident field of center frequency 200M Hz and bandwidth 100M Hz is
incident on this box along −ˆ. The backscatter RCS scattered by the box is computed
z
using, first a uniform discretization of (average edge length) h = 0.1m everywhere and
basis functions of order p = 1 everywhere. This is followed by discretizing the box
using h = 0.1m near the corners, h = 0.2m near the edges and h = 0.25m elsewhere.
The RCS is then computed using p = 1 in patches of size h = 0.1m, p = 2 where
h = 0.2m and p = 3 where h = 0.25m. The two RCS are computed in figure 6.6c.
The RCS, again, show excellent agreement.
The obvious advantage of this scheme is the ability to reduce the number of
unknowns (the spatial unknown count dropped from 776 to 226 in the case of the
box in figure 6.6a and from 11000 to 2800 in figure 6.6c ). However, notwithstanding
this gain in computational time, we would like to stress the ease of implementing the
137

mixture of orders of basis functions in the GMM scheme. Indeed, in the code, this
change is effected simply by changing an input flag. This ability to trivially mix basis
functions is unique to the GMM scheme.

6.6.4

Non conformal meshes

As before, one of the important advantages of disassociating the definition of the
basis function from the underlying tessellation of the geometry is the ability to use
non conformal meshes. This was investigated in some detail in earlier chapters. Here
we consider scattering from a PEC box of dimensions 2.0m × 0.2m × 0.5m. Two

(a) PEC box (1m × 4m × 1m)
Figure 6.6: Backscatter RCS comparison of standard GMM vs. Mixed-order GMM
(a) Box using p = 1 near edges and p = 2 in larger patches near center (b) NASA
almond discretized with GMM basis functions complete to order p = 2 near edges
and p = 1 elsewhere and (c) dielectric box using p = 1 near edges, p = 3 near the
center and p = 2 in the transition region.
138

(b) NASA thin almond
Figure 6.6 Continued: Backscatter RCS comparison of standard GMM vs. Mixedorder GMM (a) Box using p = 1 near edges and p = 2 in larger patches near center
(b) NASA almond discretized with GMM basis functions complete to order p = 2
near edges and p = 1 elsewhere and (c) dielectric box using p = 1 near edges, p = 3
near the center and p = 2 in the transition region.

139

(c) Three orders of basis functions
Figure 6.6 Continued: Backscatter RCS comparison of standard GMM vs. Mixedorder GMM (a) Box using p = 1 near edges and p = 2 in larger patches near center
(b) NASA almond discretized with GMM basis functions complete to order p = 2
near edges and p = 1 elsewhere and (c) dielectric box using p = 1 near edges, p = 3
near the center and p = 2 in the transition region.

140

halves of the box are meshed at slightly different discretization rates to obtain a
ˆ
non conformal mesh. The backscattered RCS is computed due to scattering of a x
ˆ
polarized plane wave incident along z. The RCS obtained due to scattering from a
discretization constructed on this non conformal mesh is compared against that is
obtained by constructing a fully conformal mesh on the geometry. The agreement
of the RCS shown in figure 6.7 demonstrates the ability of the GMM scheme to
seamlessly handle non conformal meshes.

141

(a) Non conformal mesh for a box

(b) Backscatter RCS comparison
Figure 6.7: Backscatter RCS comparison for scattering from a non conformally
meshed PEC box. RCS obtained by using the GMM discretization on a non conformal
mesh (shown inset) compared against that obtained by discretizing on a (standard)
conformal mesh for a thin box

142

Chapter 7
Conclusions and Future Work
This thesis introduces a novel umbrella framework for the solution of integral equations through the method of moments. The technique expands the basis function
spaces that can be used to construct MoM solutions. The GMM framework encompasses a wide variety of schemes, such as meshless methods, Nystrom schemes, point
cloud techniques as well as standard RWG based models. We have introduced a
mechanism for the near complete decoupling of the basis function definition from the
underlying tessellation, in stark contrast to the existing methods, where the basis
function is usually designed to take advantage of an existing tessellation. The standard tessellation is replaced by an overlapping tessellation and a partition of unity
function. The implicit continuity imposed by the partition of unity function allows for
the above decoupling that then facilitates the use of a wide variety of approximation
functions. The use of such functions in a traditional RWG scheme requires extremely
careful manipulations.
In addition to the GMM framework, this thesis also proposes a local Helmhotlz
decomposition based basis function scheme has resulted in a well conditioned system
over a wide range of frequencies. This has important implications in the design
of solution schemes for complex geometries involving a wide range of discretization
143

scales. Results have been presented that demonstrate the accuracy and advantages
of the technique presented here.
The final two chapters of the thesis also describe the application of the GMM
discretization scheme to the analysis of scattering from dielectric bodies and to time
domain scattering. The advantages of the GMM technique are shown in terms of
the well conditioned nature of the matrix systems in the discretization of the PMCHWT system and the low frequency stability of the discretized time domain integral
equations.
The chief contributions of this work are detailed below.
• An umbrella framework for the discretization of electromagnetic integral equations. A theoretical foundation has been laid for the development of Partition
of Unity based basis functions and their error bounds have been rigorously derived for smooth geometries. A mechanism has been described for the possible
decoupling of the geometry description and the basis function definition.
• A novel basis function design has been proposed which relies on a local Helmholtz
decomposition and a mechanism has been described for its implementation using
a variety of approximation functions.
• A detailed implementation structure has been provided for the GMM starting
from commonly available meshes where a prescription has been provided for
– The construction of patches
– The definition of functions on these patches
– Design of a partition of unity function
– Design of approximation function
– Removal of extra degrees of freedom and
– Evaluation of the matrix elements.
144

• The technique has been applied to a wide variety of examples both validating
the technique and demonstrating its prowess.
• An important consequence of the basis function definition in terms of well conditioned nature of the matrix system has been demonstrated for a very wide
variety of geometries and formulations. Specifically, the system has been shown
to result in (i) well conditioned matrix systems for PEC and dielectric objects
over more than 6 decades of frequencies and (ii) low frequency stable discretizations in time domain.
• The technique has been extended to the discretization of the PMCHWT and
M¨ller integral equations for the analysis of scattering from homogeneous diu
electric bodies.
• The GMM has been applied to the analysis of transient scattering from PEC
objects and the ability of the technique to result in a stable system at very low
frequencies where the standard discretization breaks down has been demonstrated through several examples.

7.1

Future Work

The work that has been reported in this thesis opens the door for a variety of future
research directions. We identify some of them in the remainder of this chapter.

7.1.1

Mathematical proof for the well conditioned nature of
the GMM

The GMM discretized integral equations have been shown to be well conditioned
over an extensive range of geometries and frequency scales. In particular, it was
145

demonstrated that for both PEC and dielectric objects, the numerical system has
a near-constant condition number over more than 6 decades of frequency (or correspondingly length scales). This has very important consequences in the analysis of
multi-scale objects. In addition to well conditioned matrix systems, the ability to
trivially integrate multiple basis functions in the basis space will also enable the analysis of realistic geometries. This well conditioned nature of the discretized equations
have also been translated to time domain where the GMM discretization has been
shown to result in low frequency stable systems. Again, this result has important
and far reaching consequences in the analysis of transient scattering. The immediate
future work in this direction is to develop a rigorous set of proofs for the low frequency stability of the GMM discretization and potential bounds on the condition
number of the matrix system. This can start by examining the projections of the
basis functions on the eigen values of canonical structures and estimating condition
numbers. Further proofs of the well conditioned nature can be obtained by examining
the function spaces spanned by the basis functions involved in the GMM. A study of
this nature will provide solid footing on which to design basis functions for realistic
and challenging problems involving very large scales of discretization.

7.1.2

Generalized higher order patch construction

The GMM technique developed here has been implemented for flat triangulations.
However extension to smoother geometry descriptions is extremely straightforward
and will form the first of future efforts. Higher order patches are typically defined on
a tangent triangle using shifted Sylvester polynomials [74]. To extend the GMM to
such patch definitions using the structure described in this thesis is relatively easy.
The patches can be constructed starting from nodes, with the collection of “higher
order” triangles sharing that node forming a patch. The rest of the implementation
structure will remain essentially unchanged with the sole exception of the construction
146

of matrix elements. Numerical integration on the curved triangles will have to now
account for the curvature. Similarly the analytical integration routines will need to
be redesigned. However, interesting questions arise if higher order patches are not
constructed using such a combination of triangles. Another possible design is to
construct a surface such that it coincides with the original surface at a finite set of
points but maintains a continuous surface normal. This has the advantage of both
having a smooth surface description and, at the same time, allowing for continuity of
the normal component and thereby avoiding the need for the evaluation of singular
integrals on the line as described in 4. The use of higher order patches have been
shown to lead to significant gains in accuracy.

7.1.3

Further development of the transient GMM

The application of the GMM basis functions described in this work to the discretization of the time domain integral equations leaves a lot of interesting issues to be
resolved. The main among these are the convergence properties of the inverse of the
time domain operator. Further, the well conditioned nature of the discretized frequency domain operator leads to the expectation that the GMM discretization of the
time domain integral equation will lead to a stable discrete system. This has been
shown in several numerical examples in this thesis and a rigorous mathematical proof
is a natural follow-up. Finally, the advantages of using a GMM-type discretization
based on a partition of unity framework from both a stability and accuracy point of
view can be explored.

7.1.4

Application to multi-scale geometries

The well conditioned nature of the GMM at low frequencies motivates its application
to multi-scale geometries. In addition to the use of the fast multipole method that
147

has already been applied to the GMM, the accelerated Cartesian expansions can also
be implemented trivially within a GMM framework making it an excellent choice
for the solution of scattering problems from highly multi-scale geometries. Further,
straightforward parallelization of the GMM framework will enable the analysis of
extremely large and realistic structures.

7.1.5

Domain Decomposition

The ability to include arbitrary functions in the basis function space of the GMM
scheme can be exploited in a domain decomposition framework. Domain decomposition techniques have existed in literature for a while now and involve separating
the domain of interest into smaller regions each of which are analyzed separately.
The idea of domain decomposition is popular in constructing differential equation
based solvers, and has been present in integral equation solvers for a while although
it has only recently been used for analyzing mixed scale problems. This technique
can be extended to integral equations and can be very useful for solving highly multi
scale problems. The GMM framework is excellently suited to domain decomposition
problems. The freedom of choice of basis functions allows a flexibility to use basis functions best suited for each sub domain and the partition of unity framework
provides a natural mechanism for the interaction of these domains.

7.1.6

Application to other integral equations

Finally, the GMM mechanism has been developed in this work for electromagnetic
integral equations. Their extension to the solution of integral equation formulations
of acoustic scattering problem is expected to be straightforward and will become an
area of active future research.

148

BIBLIOGRAPHY

149

BIBLIOGRAPHY

[1] R. Mittra, Ed., Computer Techniques for Electromagnetics. Oxford: Pergamon
Press, 1973.
[2] L. F. Richardson, “The approximate arithmetical solution by finite differences
of physical problems involving differential equations, with an application to the
stress in a masonary dam,” Philosophical Transactions of the Royal Society of
London: Series A, vol. 210, pp. 307–357, 1901.
[3] R. Courant, K. O. Freidrichs, and L. H., “Uber die partiellen differenzengleichungen der mathematischen physik,” Annals of Mathematics, vol. 100, pp.
32–74, 1928.
[4] P. D. Lax, “Weak solutions of nonlinear hyperbolic equations and their numerical computation,” Communications of Pure and Applied Mathematics, vol. 7,
pp. 159–193, 1954.
[5] P. D. Lax and B. Wendroff, “Systems of conservation laws,” Commun. Pure
Appl. Math., vol. 13, pp. 217–237, 1960.
[6] A. Taflove and S. Hagness, Computational Electrodynamics: The Finite Difference Time Domain Method. Boston: Artech House, 2000.
[7] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element
Methods. Springer, 2002.
[8] J. M. Jin, The Finite Element Method in Electromagnetics. New York: Wiley,
2002.
[9] I. Babuska, “Error Bounds for the Finite Element Methods,” Numerische Mathematik, vol. 333, pp. 322–333, 1971.
[10] ——, “The finite element method with lagrange multipliers,” Numerische Mathematik, pp. 179–192, 1973.
150

[11] A. Bossavit and J.-C. Verite, “A mixed fem-biem method to solve 3-d eddycurrent problems,” Magnetics, IEEE Transactions on, vol. 18, no. 2, pp. 431–
435, Mar 1982.
[12] A. Bendali, “Numerical analysis of the exterior boundary value problem for the
time- harmonic maxwell equations by a boundary finite element method part 1:
The continuous problem,” Mathematics of Computation, vol. 43, no. 167, pp.
29–46, jul 1984.
[13] ——, “Numerical analysis of the exterior boundary value problem for the timeharmonic maxwell equations by a boundary finite element method part 2: The
discrete problem,” Mathematics of Computation, vol. 43, no. 167, pp. 47–68,
jul 1984.
[14] M. L. Barton and Z. J. Cendes, “New vector finite elements for threedimensional magnetic field computation,” Journal of Applied Physics, vol. 61,
pp. 3919–3921, Apr. 1987.
[15] R. Lee and A. C. Cangellaris, “A study of discretization error in the finite
element approximation of wave solution,” IEEE Transactions on Antennas and
Propagation, vol. 40, pp. 542–549, 1992.
[16] D. Botteldooren, “Acoustical finite-difference time-domain simulation in a
quasi-cartesian grid,” Journal of the Acoustical Society of America, vol. 95,
no. 5, pp. 2313–2319, 1994.
[17] ——, “Finite-difference time-domain simulation of low-frequency room acoustic
problems,” Journal of the Acoustical Society of America, vol. 98, no. 6, pp.
3302–3308, 1995.
[18] J.-M. Jin, “Finte element-boundary element methods for electromagnetic scattering,” Ph.D. dissertation, The University of Michigan, 1989.
[19] A. F. Peterson, Mapped Vector Basis Functions for Electromagnetic Integral
Equations. Morgan & Claypool, 2006.
[20] I. Perugia, D. Sch¨tzau, and P. Monk, “Stabilized interior penalty methods for
o
the time-harmonic maxwell equations,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 41, pp. 4675–4697, September 2002.
[21] C. Lu, J. Villa, B. Shanker, and L. Kempel, “Generalized finite element method
for analyzing scattering from non-lipschitzian domains,” in Antennas and Propagation Society International Symposium 2006, IEEE, 9-14 July 2006, pp. 1783–
1786.
[22] Y.-J. Li and J.-M. Jin, “A new dual-primal domain decomposition approach for
finite element simulation of 3-d large-scale electromagnetic problems,” Antennas
and Propagation, IEEE Transactions on, vol. 55, no. 10, pp. 2803–2810, Oct.
2007.
151

[23] R. F. Harrington, “Boundary integral formulations for homogeneous material
bodies,” Journal of Electromagnetic Waves and Applications, vol. 3, pp. 1–15,
1989.
[24] ——, Field Computation by Moment Method. MacMillan Publishing Company,
1968.
[25] J. R. Mautz and R. F. Harrington, “H-field, e-field, and combined field solutions
¨
for conducting bodies of revolution,” AEU, vol. 32, pp. 157–164, 1978.
[26] ——, “Electromagnetic scattering from a homogeneous body of revolution,”
¨
AEU, vol. 33, pp. 71–80, 1979.
[27] W. C. Chew, M. Tong, and B. Hu, Integral Equation Methods for Electromagnetic and Elastic Waves, C. A. Balanis, Ed. Morgan & Claypool, 2009.
[28] R. Kress, Linear Integral Equations. Springer-Verlag, 1989.
[29] A. J. Poggio and E. K. Miller, “Integral equation solutions of three-dimensional
scattering problems,” in Computer Techniques for Electromagnetics, R. Mittra,
Ed. Oxford: Pergamon Press, 1973.
[30] J.-C. Nedelec, Acoustic and Electromagnetic Equations: Integral Represenations for Harmonic Problems.
Springer, March 2001. [Online]. Available: http://www.amazon.ca/exec/obidos/redirect?tag=citeulike0920&path=ASIN/0387951555
[31] R. Adams, “Physical and analytical properties of a stabilized electric field integral equation,” IEEE T ANTENN PROPAG, vol. 52, no. 2, pp. 362–372,
2004.
[32] A. Buffa and S. H. Christiansen, “The electric Field integral equation on Lipschitz screens: definitions and numerical approximation,” Numerische Mathematik, vol. 94, pp. 229–267, 2003.
[33] A. Buffa and R. Hiptmair, “A coercive combined field integral equation for
electromagnetic scattering,” Society, vol. 42, no. 2, pp. 621–640, 2004.
[34] S. H. Christiansen, “Discrete fredholm properties and convergence estimates
for the electric field integral equation,” Math. Comput., vol. 73, no. 245, pp.
143–167, 2004.
[35] T. Ha-Duong, On retarded potential boundary integral equations and their discretisation in Topics in Computational Wave Propagation: Direct and Inverse
Problems. Springer-Verlag,, 2003, p. 301336.
[36] R. Hiptmair, “Higher order whitney forms,” Progress in Electromagnetics Research, vol. 32, pp. 271–299, 1988.
152

[37] G. Hsiao and R. Kleinman, “Mathematical foundations for error estimation in
numerical solutions of integral equations in electromagnetics,” Antennas and
Propagation, IEEE Transactions on, vol. 45, no. 3, pp. 316–328, March 1997.
[38] R. Coifman, V. Rokhlin, and S. Wandzura, “The fast multipole method for the
wave equation: A pedestrian prescription,” IEEE Antennas Propagat. Mag.,
vol. 35, no. 3, pp. 7–12, 1993.
[39] V. Rokhlin, “Diagonal forms of translation operators for the helmholtz equation
in three dimensions,” Applied and Computational Harmonic Analysis, vol. 1, pp.
82–93, 1993.
[40] H. Huang and B. Shanker, “A (n) procedure for evaluating the yukawa potential
and the ewald summation for peroidic systems,” Submited to J. Comp. Phys.
(Technical Report Available at https://www.egr.msu.edu/ece/Technicalpapers/),
2007.
[41] L. Greengard and W. D. Gropp, “A parallel version of the fast multipole
method,” Computers and Mathematics with Applications, vol. 20, pp. 63–71,
1990.
[42] L. Greengard and V. Rokhlin, “A new version of the fast multipole method for
the laplace equation in three dimensions,” Acta Numerica, vol. 6, pp. 229–269,
1997.
[43] ——, “A fast algorithm for particle simulations,” Journal of Computational
Physics, vol. 20, pp. 63–71, 1987.
[44] W. C. Chew, J. M. Jin, C. C. Lu, E. Michielssen, and J. M. Song, “Fast solution
methods in electromagnetics,” IEEE Transactions on Antennas and Propagation, vol. 45, pp. 533–543, 1997.
[45] J. M. Song and W. C. Chew, “Multilevel fast-multipole algorithm for solving
combined field integral equations of electromagnetic scattering,” Microwave and
Optical Technology Letters, vol. 10, no. 1, pp. 14–19, 1995.
[46] J. M. Song, C. C. Lu, and W. C. Chew, “Mlfma for electromagnetic scattering
by large complex objects,” IEEE Transactions on Antennas and Propagation,
vol. 45, pp. 1488–1493, 1997.
[47] W. Chew, E. Michielssen, J. M. Song, and J. M. Jin, Eds., Fast and Efficient
Algorithms in Computational Electromagnetics. Norwood, MA, USA: Artech
House, Inc., 2001.
[48] M. Vikram, H. Huang, B. Shanker, and T. Van, “A novel wideband fmm for
fast integral equation solution of multiscale problems in electromagnetics,” Antennas and Propagation, IEEE Transactions on, vol. 57, no. 7, pp. 2094–2104,
July 2009.
153

[49] A. Ergin, B. Shanker, and E. Michielssen, “The plane-wave time-domain algorithm for the fast analysis of transient wave phenomena,” Antennas and Propagation Magazine, IEEE, vol. 41, no. 4, pp. 39–52, Aug. 1999.
[50] A. A. Ergin, B. Shanker, and E. Michielssen, “Fast evaluation of threedimensional transient wave fields using diagonal translation operators,” Journal
of Computational Physics, vol. 146, no. 1, pp. 157–180, 1998.
[51] B. Shanker, A. A. Ergin, K. Ayg¨, and E. Michielssen, “Analysis of transient
n
electromagnetic scattering phenomena using a two-level lane wave time domain
algorithm,” IEEE Trans. Antennas and Propagat., vol. 48, pp. 510–523, 2000.
[52] B. Shanker, A. Ergin, M. Lu, and E. Michielssen, “Fast analysis of transient
electromagnetic scattering phenomena using the multilevel plane wave time
domain algorithm,” Antennas and Propagation, IEEE Transactions on, vol. 51,
no. 3, pp. 628–641, March 2003.
[53] B. Shanker, A. A. Ergin, and E. Michielssen, “Analysis of transient scattering from penetrable bodies using using the multilevel plane wave time domain
algorithm,” Journal of Optical Society of America, vol. 19, pp. 716–726, 2002.
[54] G. Kobidze, J. Gao, B. Shanker, and E. Michielssen, “A fast time domain integral equation based scheme for analyzing scattering from dispersive objects,”
Antennas and Propagation, IEEE Transactions on, vol. 53, no. 3, pp. 1215–
1226, March 2005.
[55] P. Jiang, K. Yegin, S. Li, B. Shanker, and E. Michielssen, “An improved plane
wave time domain algorithm for dissipative media,” in Antennas and Propagation Society International Symposium, 2003. IEEE, vol. 3, 22-27 June 2003, pp.
563–566vol.3.
[56] G. Kobidze and B. Shanker, “Integral equation based analysis of scattering from
3d inhomogeneous anisotropic bodies,” in Antennas and Propagation Society International Symposium, 2003. IEEE, vol. 1, 22-27 June 2003, pp. 633–636vol.1.
[57] A. A. Ergin, B. Shanker, and E. Michielssen, “Analysis of transient wave scattering from rigid bodies using a burton-miller approach,” Journal of the Acoustical
Society of America, vol. 106, no. 5, pp. 2396–2404, 1999.
[58] N.-W. Chen, B. Shanker, and E. Michielssen, “Integral-equation-based analysis of transient scattering from periodic perfectly conducting structures,” in
Microwaves, Antennas and Propagation, IEE Proceedings -, vol. 150, 10 April
2003, pp. 120–124.
[59] J. H. Richmond, “A wire-grid model for scattering by conducting bodies,” IEEE
Transactions on Antennas & Propagation, vol. AP-14, pp. 782–786, 1966.
154

[60] E. K. Miller and F. J. Deadrick, “Somc computational aspects of thin-wire modeling,” in Numerical and Asymprotic Techniques in Electromagnerics, R. Mittra,
Ed. Springer-Verlag, 1975.
[61] K. S. H. Lee, L. Marin, and J. P. Castillo, “Limitations of wire- grid modeling
of a closed surface,” IEEE Transactions on Electromagnetic Compatability, vol.
EMC-18, pp. 123–129, 1976.
[62] N. C. Albertsen, J. E. Hansen, and N. E. Jensen, “Computation of radiation
from wire antennas on conducting bodies,” IEEE Transactions on An, vol. AP22, pp. 200–206, 1974.
[63] D. L. Knepp and J. Goldhirsh, “Numerical analysis of electromagnetic radiation
properties of smooth conducting bodies,” IEEE Transactions on Antennas &
Propagation, vol. AP-20, pp. 383–388, 1972.
[64] N. N. Wang, J. H. Richmond, and M. C. Gilreath, “Sinusoidal reaction formulation for radiation and scattering from conducting surfaces,” IEEE Transactions
on Antennas & Propagation, vol. AP-23, pp. 376–382, 1975.
[65] E. H. Newman and D. M. Pozar, “Electromagnetic modeling of composite wire
and surface geometries,” IEEE Transactions on Antennas & Propagation, vol.
AP-26, pp. 784–789, 1978.
[66] A. Sankar and T. Tong, “Current computation on complex structures by finiteelement method,” Electronics Letters, vol. 11, pp. 481–482, 1975.
[67] J. J. H. Wang, “Numerical analysis of three-dimensional arbitrarily-shaped conducting scatters by trilateral surface cell modelling,” Radio Science, vol. 13, pp.
947–952, 1978.
[68] J. J. H. Wang and C. Papanicolopulos, “Surface patch modeling of scatterers of
arbitrary shapes,” in The digest of the International Symposium of the Antennas
and Propagation Society, 1979.
[69] G. Jeng and A. Wexler, “Finite element, boundary integral equation analysis
of scattering problems,” in URSI Symposium on Electromagnetic Theory, 1977.
[70] A. W. Glisson and D. R. Wilton, “Simple and efficient numerical methods
for problems of electromagnetic radiation and scattering from surface,” IEEE
Transactions on Antennas and Propagation, vol. 28, pp. 593–603, 1980.
[71] S. M. Rao, D. R. Wilton, and A. W. Glisson, “Electromagnetic scattering by
surfaces of arbitrary shape,” IEEE Transactions on Antennas and Propagation,
vol. 30, pp. 408–418, 1982.
[72] J. Nedelec, “Mixed finite elements in R3 ,” Numer. Math, pp. 93–315, 1980.
155

[73] P.-A. Raviart and J. M. Thomas, Mathematical aspects of finite element methods. Springer, 1975, vol. 606, ch. A mixed finite element method for 2nd order
elliptic problems, pp. 292–315.
[74] R. D. Graglia, D. R. Wilton, and A. F. Peterson, “Higher order interpolatory
vector bases for computational electromagnetics,” IEEE Transactions on Antennas and Propagation, vol. 45, pp. 329–342, 1997.
[75] W. Cai, T. Yu, H. Wang, and Y. Yu, “High-order mixed rwg basis functions for
electromagnetic applications,” IEEE Transactions on Microwave Theory and
Techniques, vol. 49, no. 7, pp. 1295–1303, Jul 2001.
[76] J. M. Melenk and I. Babuska, “The partition of unity finite element method
: Basic theory and applications,” Computer Methods in Applied Mechanical
Engineering, vol. 139, pp. 289–314, 1996.
[77] C. Lu and B. Shanker, “Generalized finite element method for vector electromagnetic problems,” IEEE Transactions on Antennas and Propagation, vol. 55,
no. 5, pp. 1369–1381, May 2007.
[78] N. Nair and B. Shanker, “Generalized method of moments: A novel
discretization scheme for integral equations,” to appear in IEEE Transactions on Antennas and Propagation; technical report available at
https://www.egr.msu.edu/ece/technical papers, 2010.
[79] D. Wilton and S. Govind, “Incorporation of edge conditions in moment method
solutions,” IEEE Transactions on Antennas and Propagation, vol. 25, no. 6, pp.
845–850, Nov 1977.
[80] W. Brown and D. Wilton, “Singular basis functions and curvilinear triangles in
the solution of the electric field integral equation,” Antennas and Propagation,
IEEE Transactions on, vol. 47, no. 2, pp. 347–353, Feb 1999.
[81] R. Graglia and G. Lombardi, “Singular higher order complete vector bases for
finite methods,” IEEE Transactions on Antennas and Propagation, vol. 52,
no. 7, pp. 1672–1685, July 2004.
[82] A. F. Peterson and M. M. Bibby, “High-order numerical solutions of the mfie for
the linear dipole,” IEEE Transactions on Antennas and Propagation, vol. 52,
pp. 2684–2691, 2004.
[83] M. M. Bibby, A. F. Peterson, and C. M. Coldwell, “High order representations for singular currents at corners,” IEEE Transactions on Antennas and
Propagation, vol. 56, pp. 2277–2287, 2008.
[84] N. Nair, C. Lu, and B. Shanker, “Differential and integral equation solvers
based on generalized moments and partitions of unity,” in Electromagnetics in
Advanced Applications, 2007. ICEAA 2007. International Conference on, Sept.
2007, pp. 974–977.
156

[85] W. C. Chew, Waves and Fields in Inhomogeneous Media.
Nostrand Reinhold, 1990.

New York: Van

[86] A. F. Peterson, S. L. Ray, and R. Mittra, Computational Methods for Electromagnetics, 1st ed., ser. IEEE/OUP series on electromagnetic wave theory. New
York: IEEE Press, 1998.
[87] J. A. Stratton, Electromagnetic theory, ser. International series in physics. New
York: McGraw-Hill, 1941.
[88] A. W. Maue, “The formulation of a general diffraction problem by an integral
equation,” Zeitschrift fur Physik A, vol. 126, pp. 601–618, 1949.
[89] L. N. Medgyesi-Mitschang, J. M. Putnam, and M. B. Gedera, “Generalized
method of moments for three-dimensional penetrable scatteres,” J. Opt. Soc.
Am. A, vol. 11, pp. 1383–1398, 1994.
[90] K. A. Michalski and D. Zheng, “Electromagnetic scattering and radiation by
surfaces of arbitrary shape in layered media,” IEEE Transactions on Antennas
and Propagation, vol. 38, pp. 335–344, 1990.
[91] ——, “Electromagnetic scattering and radiation by surfaces of arbitrary shape
in layered media. ii. implementation and results for contiguous half-spaces,”
IEEE Transactions on Antennas and Propagation, vol. 38, pp. 345–352, 1990.
[92] S. Rao and D. Wilton, “Transient scattering by conducting surfaces of arbitrary shape,” IEEE Transactions on Antennas and
Propagation, vol. 39, no. 1, pp. 56–61, 1991. [Online]. Available:
http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=64435
[93] G. Vecchi, “Loop-Star Decomposition of Basis Functions in the Discretization
of the EFIE,” IEEE Transactions on Antennas and Propagation, vol. 47, no. 2,
pp. 339–346, 1999.
[94] T. J. Cui and W. C. Chew, “Accurate Analysis of Wire Structures From VeryLow Frequency to Microwave Frequency,” IEEE Transactions on Antennas and
Propagation, vol. 50, no. 3, pp. 301–307, 2002.
[95] F. Andriulli, K. Cools, H. Bagci, F. Olyslager, A. Buffa, S. Christiansen, and
E. Michielssen, “A multiplicative calderon preconditioner for the electric field
integral equation,” IEEE Transactions on Antennas and Propagation, vol. 56,
no. 8, pp. 2398–2412, Aug. 2008.
[96] A. Bossavit, Computational Electromagnetism. Academic Press, 1998.
[97] R. J. Adams, N. J. Champagne, and S. Member, “A Numerical Implementation of a Modified Form of the Electric Field Integral Equation,” IEEE Transactions on Antennas and Propagation,
157

vol. 52, no. 9, pp. 2262–2266, September 2004. [Online]. Available:
http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1331612
[98] O. Tuncer, C. Lu, N. Nair, B. Shanker, and L. C. Kempel, “Further
development of the vector generalized finite element method and its hybridization with boundary integrals,” accepted for publication in IEEE
Transactions on Antennas and Propagation; technical report available at
https://www.egr.msu.edu/ece/technical papers, 2010.
[99] I. Babuska and J. M. Melenek, “The partition of unity method,” International
Journal for Numerical Methods in Engineering, vol. 40, pp. 727–858, 1997.
[100] R. Scharstein, “Helmholtz decomposition of surface electric current in electromagnetic scattering problems,” System Theory, 1991. Proceedings., TwentyThird Southeastern Symposium on, pp. 424–426, Mar 1991.
[101] N. Nair and B. Shanker, “An extension of the generalized method of moments
to higher order geometries,” in Proceedings of the 25th International Review of
Progress in Applied Computational Electromagnetics, 2009.
[102] A. Buffa and S. H. Christiansen, “A dual finite element complex on the barycentric refinement,” Mathematics of Computation, vol. 76, no. 260, pp. 1743–1769,
2007.
[103] C. Delgado, E. Garcia, F. Catedra, and R. Mittra, “Generation of characteristic
basis functions defined over large surfaces by using a multilevel approach,” IEEE
Transactions on Antennas and Propagation, vol. 57, no. 4, pp. 1299–1301, Apr.
2009.
[104] M. Khayat and D. Wilton, “Numerical evaluation of singular and near-singular
potential Integrals,” IEEE Transactions on Antennas and Propagation,
vol. 53, no. 10, pp. 3180–3190, October 2005. [Online]. Available:
http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=1514571
[105] L. Knockaert, “A general gauss theorem for evaluating singular integrals over
polyhedral domains,” Electromagnetics, vol. vol, pp. 11pp269–280, 1991.
[106] S. Jarvenpaa, M. Taskinen, and P. Yl¨-oijala, “Singularity extraction technique
a
for integral equation methods with higher order basis functions on plane triangles and tetrahedra,” International Journal for Numerical Methods in Engineering, vol. 58, pp. 1149–1165, 2003.
[107] 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
nystrom discretization,” Journal of Computational Physics, vol. 146, no. 2, pp.
627–663, 1998.
158

[108] M. Vikram and B. Shanker, “An incomplete review of fast multipole methods
from static to wideband as applied to problems in computational electromagnetics,” Applied Computational Electromagnetics Society Journal, vol. 27, p. 79,
2009.
[109] A. Chatterjee, J. Jin, and J. L. Volakis, “Edge-based finite elements and vector
abc’s applied to 3-d scattering,” IEEE Transactions on Antennas and Propagation, vol. 41, pp. 221–226, 1993.
[110] S.-C. Lee, M. N. Vouvakis, and J.-F. Lee, “A non-overlapping domain decomposition method with non-matching grids for modeling large finite antenna arrays,” J. Comput. Phys., vol. 203, no. 1, pp. 1–21, 2005.
[111] C. Muller, Foundations of the Mathematical Theory of Electromagnetic Waves.
Berlin: Springer Verlag, 1969.
[112] P. Yla-Oijala and M. Taskinen, “Well-Conditioned M¨ller Formulation for Elecu
tromagnetic Scattering by Dielectric Objects,” October, vol. 53, no. 10, pp.
3316–3323, 2005.
[113] F. P. Andriulli, K. Cools, F. Olyslager, and E. Michielssen, “Time Domain
´
CalderOn Identities and Their Application to the Integral Equation Analysis of
Scattering by PEC Objects Part II: Stability,” IEEE Transactions on Antennas
and Propagation, vol. 57, no. 8, pp. 2365–2375, August 2009. [Online]. Available:
http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5067343
[114] B. Shanker, A. A. Ergin, K. Agy¨n, and E. Michielssen, “Analysis of transient
u
electromagnetic scattering from closed surfaces using a combined field integral
equation,” IEEE Transactions on Antennas and Propagation, vol. 48, pp. 1064–
1074, 2000.
[115] A. A. Ergin, B. Shanker, and E. Michielssen, “Transient analysis of acoustic
scattering using marching-on-in-time with plane wave time domain algorithm,”
J. Acoustical Soc. Am., vol. 106, pp. 2405–2416, 1999.
[116] D. Weile, G. Pisharody, N.-W. Chen, B. Shanker, and E. Michielssen, “A novel
scheme for the solution of the time-domain integral equations of electromagnetics,” Antennas and Propagation, IEEE Transactions on, vol. 52, no. 1, pp.
283–295, Jan. 2004.
[117] J. J. Knab, “Interpolation of bandlimited functions using the approximate prolate series,” IEEE Trans. Information Theory, vol. 25, pp. 717–720, 1979.
[118] A. Pray, N. V. Nair, and B. Shanker, “A generalized discretization scheme for
transient integral equations,” in to appear in Proceedings of the International
Symposium on Antennas and Propagation, 2010.
159

[119] A. Pray, N. Nair, and B. Shanker, “Evaluation of retarded potential integrals
using an approximation to the green’s function,” To be submitted to IEEE
Transactions on Antennas and Propagation, 2010.

160