# VIE-FG-FFT for Analyzing EM Scattering from Inhomogeneous Nonmagnetic Dielectric Objects

Authors: Shu-Wen Chen; Hou-Xing Zhou; Wei Hong; Jia-Ye Xie
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

A new realization of the volume integral equation (VIE) in combination with the fast Fourier transform (FFT) is established by fitting Green’s function (FG) onto the nodes of a uniform Cartesian grid for analyzing EM scattering from inhomogeneous nonmagnetic dielectric objects. The accuracy of the proposed method is the same as that of the P-FFT and higher than that of the AIM and the IE-FFT especially when increasing the grid spacing size. Besides, the preprocessing time of the proposed method is obviously less than that of the P-FFT for inhomogeneous nonmagnetic dielectric objects. Numerical examples are provided to demonstrate the accuracy and efficiency of the proposed method.


---

## Body

## 1. Introduction

The volume integral equation (VIE) method [1] based on the method of moments (MoM) [2] is one of the efficient methods to analyze electromagnetic (EM) scattering from inhomogeneous dielectric objects. As is well known, for the traditional VIE-MoM, both the storage requirement and the computational complexity of a matrix-vector multiplication when an iterative method is applied are proportional to O
(
N
2
), where N denotes the number of unknowns. Therefore, the VIE-MoM is not suitable for the direct analysis of EM scattering from electrically large and inhomogeneous dielectric objects.One of approaches for improving the efficiency of the VIE-MoM is the VIE in combination with the fast Fourier transform (FFT), and it already has several implementations, such as the VIE-AIM [3, 4], the VIE-P-FFT [5, 6], and the VIE-IE-FFT [7, 8], which are simply called the FFT-based methods. These implementations are all transplanted from the corresponding versions [9–12] for the surface integral equation (SIE) [13]. Not long ago, a new realization, the FG-FFT, of the SIE in combination with the FFT for the electric field integral equation (EFIE) was proposed [14] and soon extended to the combined field integral equation (CFIE) [15].In this paper, the FG-FFT for the SIE will be extended to the VIE for analyzing EM scattering from inhomogeneous nonmagnetic dielectric objects, and resultant method is simply called the VIE-FG-FFT. The remainder of this paper is organized as follows. In Section2, the VIE-FG-FFT is presented in detail. In Section 3, some numerical examples are provided to demonstrate the accuracy and efficiency of the VIE-FG-FFT. Finally, the conclusion is given in Section 4. In this paper, the time convention e
j
ω
t is assumed and suppressed.


## 2. Formulation

### 2.1. The Volume Integral Equation

The permittivity and permeability of the free space are denoted byϵ
0 and μ
0, respectively. Let V denote the volumetric domain occupied by an inhomogeneous nonmagnetic dielectric object with relative permittivity ϵ
r and relative permeability μ
r
=
1 (meaning nonmagnetic).LetE
→
i be the incident electric field and E
→
s the scattered electric field; then the total electric field E
→
tot can be expressed as the sum of E
→
i and E
→
s:
(1)
E
→
tot
=
E
→
i
+
E
→
s
,
and the volume integral equation (VIE) on V for the total electric field can be rigorously expressed as [16]
(2)
E
→
tot
r
→
=
E
→
i
r
→
+
k
0
2
∫
V
‍
ε
r
r
→
′
-
1
E
→
tot
r
→
′
·
G
¯
¯
r
→
,
r
→
′
d
v
′
,
where r
→
∈
V and k
0
=
ω
μ
0
ε
0, which implies that
(3)
E
→
s
r
→
=
∫
V
k
0
2
+
∇
∇
·
G
r
→
,
r
→
′
κ
r
→
′
D
→
tot
r
→
′
ε
0
d
v
′
,
where D
→
tot is the electric flux density of E
→
tot and
(4)
κ
r
→
′
=
ε
r
r
→
′
-
1
ε
r
r
→
′
.
Therefore, we have
(5)
E
→
s
r
→
=
k
0
2
∫
V
G
r
→
,
r
→
′
κ
r
→
′
D
→
tot
r
→
′
ε
0
d
v
′
W
W
i
i
W
-
∇
∫
∂
V
G
r
→
,
r
→
′
κ
r
→
′
n
^
∂
V
·
D
→
tot
r
→
′
ε
0
d
s
′
WWiiW
+
∇
∫
V
G
r
→
,
r
→
′
κ
r
→
′
∇
′
·
D
→
tot
(
r
→
′
)
ε
0
WWWWWWWWWWW
∇
′
·
D
→
tot
r
→
′
ε
0
+
D
→
tot
r
→
′
ε
0
·
∇
′
κ
r
→
′
d
v
′
,
where ∂
V is the outer boundary surface of V.It should be pointed out that the second term in the right-hand side of (5) cannot be ignored in the strict sense (see the second paragraph of Section 2 in [17]). However, this term will force one to introduce “half” basis function, which will be seen in the following.


### 2.2. Buiding the MoM Model

The electric flux densityD
→
tot can be chosen as the unknown function because it is continuous along the normal direction of the medium interface. After V is discretized by using tetrahedrons, D
→
tot
/
ε
0 can be expanded with the SWG functions [1]:
(6)
D
→
tot
r
→
ε
0
=
∑
m
=
1
N
F
‍
x
m
F
f
→
m
F
r
→
+
∑
n
=
1
N
H
‍
x
n
H
f
→
n
H
r
→
,
where “F” and “H” mean “full SWG function” and “half SWG function,” respectively, and the total number of the basis functions is N
=
N
F
+
N
H. A full SWG function f
→
n
F is defined on the union V
n
F
=
V
n
F
,
+
∪
V
n
F
,
- of a pair of tetrahedrons V
n
F
,
+ and V
n
F
,
- that share a common face as follows:
(7)
f
→
n
F
r
→
=
f
→
n
F
,
+
r
→
=
S
n
3
V
n
+
r
→
-
r
→
n
,
free
+
,
r
→
∈
V
n
F
,
+
f
→
n
F
,
-
r
→
=
S
n
3
V
n
-
r
→
-
r
→
n
,
free
-
,
r
→
∈
V
n
F
,
-
,
whose geometry is shown in Figure 1, while a half SWG function f
→
n
H is defined only on a single tetrahedron V
n
H. For convenience, we use [
V
n
F
] to denote the common face of V
n
F
,
+ and V
n
F
,
-, as shown in Figure 1.Figure 1
The geometry of the SWG function.In fact, half SWG functions are applied only at the outer boundary ofV, and hence the number N
H is just that of the triangular elements over ∂
V. In this paper, we have a convention that a half SWG function f
→
n
H is so defined on V
n
H that f
→
n
H
(
r
→
)
·
n
^
(
r
→
)
=
1 when r
→
∈
[
V
n
H
]
∶
=
V
n
H
∩
∂
V.When (6) is substituted into (1) and after the Galerkin procedure is applied, the VIE-MoM matrix equation can be built as follows:
(8)
A
X
=
B
,
where
(9)
A
=
a
m
n
F
F
N
F
×
N
F
a
m
n
F
H
N
F
×
N
H
a
m
n
H
F
N
H
×
N
F
a
m
n
H
H
N
H
×
N
H
;
B
=
b
m
F
N
F
×
1
b
m
H
N
H
×
1
is
the
incident
vector
;
X
=
x
n
F
N
F
×
1
x
n
H
N
H
×
1
is
the
unknown
vector
.Note that, in our VIE-MoM model, the restriction ofκ
(
r
→
) on a tetrahedral region will be considered as a constant, and the following symbol abbreviations will be applied:
(10)
κ
n
F
,
±
∶
=
κ
r
→
|
V
n
F
,
±
,
κ
n
H
∶
=
κ
r
→
|
V
n
H
.
Then the elements of matrix A have the following expressions:
(11)
a
m
n
F
F
=
∫
V
m
F
f
→
m
F
(
r
→
)
·
f
→
n
F
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
F
F
=
-
k
0
2
∫
V
m
F
d
v
f
→
m
F
(
r
→
)
·
∫
V
n
F
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
F
r
→
′
d
v
′
a
m
n
F
F
=
+
∫
V
m
F
d
v
∇
·
f
→
m
F
(
r
→
)
I
m
n
F
,
(
1
)
(
r
→
)
+
I
m
n
F
,
(
2
)
(
r
→
)
a
m
n
F
F
∶
=
a
m
n
F
F
,
1
-
k
0
2
a
m
n
F
F
,
2
+
a
m
n
F
F
,
3
,
(12)
a
m
n
H
F
=
∫
V
m
H
f
→
m
H
(
r
→
)
·
f
→
n
F
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
H
F
=
-
k
0
2
∫
V
m
H
d
v
f
→
m
H
(
r
→
)
·
∫
V
n
F
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
F
r
→
′
d
v
′
a
m
n
H
F
=
+
∫
V
m
H
d
v
∇
·
f
→
m
H
(
r
→
)
I
m
n
F
,
(
1
)
(
r
→
)
+
I
m
n
F
,
(
2
)
(
r
→
)
a
m
n
H
F
=
+
5
-
∫
[
V
n
H
]
d
s
I
m
n
F
,
(
1
)
(
r
→
)
+
I
m
n
F
,
(
2
)
(
r
→
)
a
m
n
H
F
∶
=
a
m
n
H
F
,
1
-
k
0
2
a
m
n
H
F
,
2
+
a
m
n
H
F
,
3
,
where
(13)
I
m
n
F
,
1
r
→
=
∫
V
n
F
G
r
→
,
r
→
′
κ
r
→
′
∇
′
·
f
→
n
F
r
→
′
d
v
′
,
I
m
n
F
,
2
r
→
=
∫
V
n
F
G
r
→
,
r
→
′
κ
n
F
,
-
-
κ
n
F
,
+
d
s
′
.
Note that a
m
n
F
F
,
1
=
0 in (11) when V
m
F and V
n
F are properly separated and that a
m
n
H
F
,
1
=
0 in (12) when V
m
H and V
n
F are properly separated:
(14)
a
m
n
F
H
=
∫
V
m
F
f
→
m
F
(
r
→
)
·
f
→
n
H
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
F
H
=
-
k
0
2
∫
V
m
F
d
v
f
→
m
F
(
r
→
)
·
∫
V
n
H
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
H
r
→
′
d
v
′
a
m
n
F
H
=
+
∫
V
m
F
d
v
∇
·
f
→
m
F
(
r
→
)
I
m
n
H
,
(
1
)
(
r
→
)
-
I
m
n
H
,
(
2
)
(
r
→
)
a
m
n
F
H
∶
=
a
m
n
F
H
,
1
-
k
0
2
a
m
n
F
H
,
2
+
a
m
n
F
H
,
3
,
(15)
a
m
n
H
H
=
∫
V
m
H
f
→
m
H
(
r
→
)
·
f
→
n
H
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
H
H
=
-
k
0
2
∫
V
m
H
d
v
f
→
m
H
(
r
→
)
·
∫
V
n
H
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
H
r
→
′
d
v
′
a
m
n
H
H
=
+
∫
V
m
H
d
v
∇
·
f
→
m
H
r
→
I
m
n
H
,
1
r
→
-
I
m
n
H
,
2
r
→
a
m
n
H
H
=
+
5
-
∫
V
m
H
d
s
∇
·
f
→
m
H
r
→
I
m
n
H
,
1
r
→
-
I
m
n
H
,
2
r
→
a
m
n
H
H
∶
=
a
m
n
H
H
,
1
-
k
0
2
a
m
n
H
H
,
2
+
a
m
n
H
H
,
3
,
where
(16)
I
m
n
H
,
1
r
→
=
∫
V
n
H
G
r
→
,
r
→
′
κ
r
→
′
∇
′
·
f
→
n
H
r
→
′
d
v
′
,
I
m
n
H
,
2
r
→
=
∫
V
n
H
d
s
′
G
r
→
,
r
→
′
κ
r
→
′
.
Note that a
m
n
F
H
,
1
=
0 in (14) when V
m
F and V
n
H are properly separated and that a
m
n
H
H
,
1
=
0 in (15) when V
m
H and V
n
H are properly separated.


### 2.3. The Frame for the VIE-FG-FFT

In this section, the FG-FFT technology is introduced into the VIE-MoM for both reducing memory requirement and improving computational efficiency.The entire MoM matrixA can be split into two parts: the near-field matrix A
near and the far-field matrix A
far, and
(17)
A
=
A
-
A
far
+
A
far
≃
A
near
+
A
far
,
where the A
near is a sparse matrix (the identification of a near-element will be presented in the second paragraph of Section 3), which is obtained by forcing all “far elements” of A
-
A
far to be equal to zero, and A
far can ultimately expressed in such a form as follows:
(18)
A
far
=
-
k
0
2
Π
→
f
·
G
Π
→
c
T
+
Π
f
G
Π
c
T
,
where Π
→
f, Π
→
c, Π
f, and Π
c are all sparse matrices which will be constructed in Section 2.4 and where the head mark “→” implies matrix elements being 3D vectors; G is a triple Toeplitz matrix related to Green’s function; the superscript T indicates matrix transpose.When an iterative solver is applied, the matrix-vector product will be performed by means of(19)
A
x
≃
A
near
x
+
A
far
x
,
where A
near
x is directly calculated, while A
far
x can be speeded up by means of the FFT through (18). In this way, the VIE-FG-FFT can reduce the memory requirement and the computational complexity to O
(
N
) and O
(
N
log
⁡
(
N
)
) theoretically (also see [9] for more detailed analysis).


### 2.4. Fitting Green’s Function

In this section, matricesΠ
→
f, Π
→
c, Π
f, and Π
c in (18) will be constructed. First, let a uniform Cartesian grid enclose the given volumetric region V. Use h
x, h
y, and h
z to denote the three grid spacing sizes in the directions x
^, y
^, and z
^, respectively. In this paper, however, the convention h
∶
=
h
x
=
h
y
=
h
z is always selected.In a uniform Cartesian grid, anexpansion box (or simplybox) C is defined as a cube-like collection composed of (
M
x
+
1
)
×
(
M
y
+
1
)
×
(
M
z
+
1
) nodes. When M
x
=
M
y
=
M
z
=
M, C includes M
+
1
3 nodes and M is called itsexpansion order (or simplyorder), written as |
C
|. Figure 2 illustrates a box of order 2.Figure 2
An expansion box of order 2.Assume thatC
n is a box centered at c
n of radius r
n and q
→ is a fixed point within C
n. Let Green’s function be represented into
(20)
G
r
→
,
q
→
=
∑
v
∈
C
n
‍
π
v
,
C
n
q
→
G
r
→
,
v
,
where the coefficients π
v
,
C
n
q
→ are to be determined and r
→ is an arbitrary point outside the box C
n, as shown in Figure 3.Figure 3
A 3D representation of the matching Green’s function.Now we select a spherical surface with the center atc
n and the radius R
n slightly larger than r
n, and select a group of sample points {
r
→
t
}
t
=
1
T on the surface [18], where the number T ought to be significantly greater than that of the nodes within C
n. In this paper, R
n
-
r
n
=
0.15 wavelength, which was adopted in [15]. Then we match Green’s function values at these sample points. The resultant systems of matrix equations can be solved by the least-square method [19]. In this paper, the order of a box is selected to be 2, and at this time the number of sample points on the testing spherical surface is selected to be 120.Without loss of generality, assume thatV
m
F
,
± and V
m
H (V
n
F
,
± and V
n
H) are in box C
m (C
n). Calculations of MoM matrix elements ultimately come down to calculations on the Gaussian points. We specify using N
G Gaussian points in a tetrahedral region and M
G Gaussian points in a triangular region. For convenience, the correspondences between integral regions and Gaussian points and weights that will be used are all listed in Table 1. When the testing element and the source element are separated by a proper distance, all the first terms in the right-hand side of (11), (12), (14), and (15) will vanish. At this case, we can easily obtain the following approximate expressions:
(21)
a
m
n
F
F
,
2
=
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
f
→
m
F
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
F
,
c
‍
‍
,
a
m
n
F
F
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
F
,
f
G
u
,
v
π
v
,
C
n
f
→
n
F
,
c
,
a
m
n
H
F
,
2
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
→
u
,
C
m
f
→
m
H
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
F
,
c
,
a
m
n
H
F
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
H
,
f
G
u
,
v
π
v
,
C
n
f
→
n
F
,
c
,
a
m
n
F
H
,
2
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
→
u
,
C
m
f
→
m
F
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
H
,
c
,
a
m
n
F
H
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
F
,
f
G
u
,
v
π
v
,
C
n
f
→
n
H
,
c
,
a
m
n
H
H
,
2
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
→
u
,
C
m
f
→
n
H
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
H
,
c
,
a
m
n
H
H
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
H
,
f
G
(
u
,
v
)
π
v
,
C
n
f
→
n
H
,
c
,
where the coefficients are calculated by using the following formulae:
(22)
π
→
u
,
C
m
f
→
m
F
,
f
=
∑
i
=
1
N
G
‍
w
i
f
→
m
F
,
+
p
→
i
+
π
u
,
C
m
p
→
i
+
+
f
→
m
F
,
-
p
→
i
-
π
u
,
C
m
p
→
i
-
,
π
u
,
C
m
f
→
m
F
,
f
=
∑
i
=
1
N
G
‍
w
i
∇
·
f
→
m
F
,
+
(
p
→
i
+
)
π
u
,
C
m
p
→
i
+
WWWWii
.
W
+
∇
·
f
→
m
F
,
-
p
→
i
-
π
u
,
C
m
p
→
i
-
,
π
→
v
,
C
n
f
→
n
F
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
F
,
+
f
→
n
F
,
+
q
→
j
+
π
v
,
C
n
q
→
j
+
WWWWii
.
W
+
κ
n
F
,
-
f
→
n
F
,
-
q
→
j
-
π
v
,
C
n
q
→
j
-
,
π
v
,
C
n
f
→
n
F
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
F
,
+
∇
·
f
→
n
F
,
+
q
→
j
+
π
v
,
C
n
q
→
j
+
WWWWWiW
+
κ
n
F
,
-
∇
·
f
→
n
F
,
-
q
→
j
-
π
v
,
C
n
q
→
j
-
π
v
,
C
n
f
→
n
F
,
c
=
+
∑
j
=
1
M
G
w
^
j
κ
n
F
,
-
-
κ
n
F
,
+
π
v
,
C
n
t
→
j
‍
,
π
→
u
,
C
m
f
→
m
H
,
f
=
∑
i
=
1
N
G
‍
w
i
f
→
m
H
p
→
i
π
u
,
C
m
p
→
i
,
π
u
,
C
m
f
→
m
H
,
f
=
∑
i
=
1
N
G
‍
w
i
∇
·
f
→
m
H
p
→
i
π
u
,
C
n
p
→
i
-
∑
i
=
1
M
G
‍
w
^
i
π
u
,
C
m
s
→
i
,
π
→
v
,
C
n
f
→
n
H
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
H
f
→
n
H
q
→
j
π
v
,
C
n
q
→
j
,
π
v
,
C
n
f
→
n
H
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
H
∇
·
f
→
n
H
q
→
j
π
v
,
C
n
q
→
j
-
∑
j
=
1
M
G
‍
w
^
j
κ
n
H
π
v
,
C
n
t
→
j
.Table 1
Correspondences between integral regions and Gaussian points and weights.
Integral region
V
m
F
,
±
V
m
H
[
V
m
F
]
V
n
F
,
±
V
n
H
[
V
n
H
]
Gaussian points
{
p
→
i
±
}
i
=
1
N
G
{
p
→
i
}
i
=
1
N
G
{
s
→
i
}
i
=
1
M
G
{
q
→
j
±
}
j
=
1
N
G
{
q
→
j
}
j
=
1
N
G
{
t
→
j
}
j
=
1
M
G
Gaussian weights
{
w
i
}
i
=
1
N
G
{
w
^
i
}
i
=
1
M
G
{
w
j
}
j
=
1
N
G
{
w
^
j
}
j
=
1
M
GSubstituting (21) into (11)–(16), we can easily obtain formula (18). It should be pointed out that performing matrix-vector product A
far
x once based on (18) requires at least 8 FFTs.


### 2.5. On the Preprocessing Time

Here, the preprocessing time means the time taken by generating the coefficients in (21). Clearly, the preprocessing time is proportional to the corresponding computational complexity. For general cases, from the coefficient formulae in Section 2.4, it can be evaluated that the number N
FG-FFT of required float multiplications is about
(23)
N
FG-FFT
W
=
N
F
22
N
G
+
5
M
G
+
24
+
N
H
11
N
G
+
2
M
G
+
15
.In the VIE-P-FFT, the corresponding coefficients are generated by utilizing “projecting coefficients,” which requires solving the least-square problem with multiple right-hand terms. If the order of a box is selected to be 2, containing 27 grid points, and the number of the sample points on the testing spherical surface is selected to beT, then the matrix of the least-square problem in whether the VIE-P-FFT or the VIE-FG-FFT is a T
×
27 complex matrix. Further, if the SVD process of the matrix is ignored, then the least-square solution for a right-hand term requires about T
2
+
2
7
2 complex multiplications or about 3
(
T
2
+
2
7
2
) real multiplications because a complex multiplication requires at least 3 real multiplications. It is easily known from the projection scheme introduced in [20] that, for generating the corresponding coefficients, the number N
P-FFT of required float multiplications is about
(24)
N
P-FFT
=
24
N
T
2
+
2
7
2
,
which is usually much larger than N
FG-FFT.It can be concluded from the above analysis that the preprocessing time of the VIE-FG-FFT is obviously less than that of the VIE-P-FFT for inhomogeneous nonmagnetic dielectric objects.


## 2.1. The Volume Integral Equation

The permittivity and permeability of the free space are denoted byϵ
0 and μ
0, respectively. Let V denote the volumetric domain occupied by an inhomogeneous nonmagnetic dielectric object with relative permittivity ϵ
r and relative permeability μ
r
=
1 (meaning nonmagnetic).LetE
→
i be the incident electric field and E
→
s the scattered electric field; then the total electric field E
→
tot can be expressed as the sum of E
→
i and E
→
s:
(1)
E
→
tot
=
E
→
i
+
E
→
s
,
and the volume integral equation (VIE) on V for the total electric field can be rigorously expressed as [16]
(2)
E
→
tot
r
→
=
E
→
i
r
→
+
k
0
2
∫
V
‍
ε
r
r
→
′
-
1
E
→
tot
r
→
′
·
G
¯
¯
r
→
,
r
→
′
d
v
′
,
where r
→
∈
V and k
0
=
ω
μ
0
ε
0, which implies that
(3)
E
→
s
r
→
=
∫
V
k
0
2
+
∇
∇
·
G
r
→
,
r
→
′
κ
r
→
′
D
→
tot
r
→
′
ε
0
d
v
′
,
where D
→
tot is the electric flux density of E
→
tot and
(4)
κ
r
→
′
=
ε
r
r
→
′
-
1
ε
r
r
→
′
.
Therefore, we have
(5)
E
→
s
r
→
=
k
0
2
∫
V
G
r
→
,
r
→
′
κ
r
→
′
D
→
tot
r
→
′
ε
0
d
v
′
W
W
i
i
W
-
∇
∫
∂
V
G
r
→
,
r
→
′
κ
r
→
′
n
^
∂
V
·
D
→
tot
r
→
′
ε
0
d
s
′
WWiiW
+
∇
∫
V
G
r
→
,
r
→
′
κ
r
→
′
∇
′
·
D
→
tot
(
r
→
′
)
ε
0
WWWWWWWWWWW
∇
′
·
D
→
tot
r
→
′
ε
0
+
D
→
tot
r
→
′
ε
0
·
∇
′
κ
r
→
′
d
v
′
,
where ∂
V is the outer boundary surface of V.It should be pointed out that the second term in the right-hand side of (5) cannot be ignored in the strict sense (see the second paragraph of Section 2 in [17]). However, this term will force one to introduce “half” basis function, which will be seen in the following.


## 2.2. Buiding the MoM Model

The electric flux densityD
→
tot can be chosen as the unknown function because it is continuous along the normal direction of the medium interface. After V is discretized by using tetrahedrons, D
→
tot
/
ε
0 can be expanded with the SWG functions [1]:
(6)
D
→
tot
r
→
ε
0
=
∑
m
=
1
N
F
‍
x
m
F
f
→
m
F
r
→
+
∑
n
=
1
N
H
‍
x
n
H
f
→
n
H
r
→
,
where “F” and “H” mean “full SWG function” and “half SWG function,” respectively, and the total number of the basis functions is N
=
N
F
+
N
H. A full SWG function f
→
n
F is defined on the union V
n
F
=
V
n
F
,
+
∪
V
n
F
,
- of a pair of tetrahedrons V
n
F
,
+ and V
n
F
,
- that share a common face as follows:
(7)
f
→
n
F
r
→
=
f
→
n
F
,
+
r
→
=
S
n
3
V
n
+
r
→
-
r
→
n
,
free
+
,
r
→
∈
V
n
F
,
+
f
→
n
F
,
-
r
→
=
S
n
3
V
n
-
r
→
-
r
→
n
,
free
-
,
r
→
∈
V
n
F
,
-
,
whose geometry is shown in Figure 1, while a half SWG function f
→
n
H is defined only on a single tetrahedron V
n
H. For convenience, we use [
V
n
F
] to denote the common face of V
n
F
,
+ and V
n
F
,
-, as shown in Figure 1.Figure 1
The geometry of the SWG function.In fact, half SWG functions are applied only at the outer boundary ofV, and hence the number N
H is just that of the triangular elements over ∂
V. In this paper, we have a convention that a half SWG function f
→
n
H is so defined on V
n
H that f
→
n
H
(
r
→
)
·
n
^
(
r
→
)
=
1 when r
→
∈
[
V
n
H
]
∶
=
V
n
H
∩
∂
V.When (6) is substituted into (1) and after the Galerkin procedure is applied, the VIE-MoM matrix equation can be built as follows:
(8)
A
X
=
B
,
where
(9)
A
=
a
m
n
F
F
N
F
×
N
F
a
m
n
F
H
N
F
×
N
H
a
m
n
H
F
N
H
×
N
F
a
m
n
H
H
N
H
×
N
H
;
B
=
b
m
F
N
F
×
1
b
m
H
N
H
×
1
is
the
incident
vector
;
X
=
x
n
F
N
F
×
1
x
n
H
N
H
×
1
is
the
unknown
vector
.Note that, in our VIE-MoM model, the restriction ofκ
(
r
→
) on a tetrahedral region will be considered as a constant, and the following symbol abbreviations will be applied:
(10)
κ
n
F
,
±
∶
=
κ
r
→
|
V
n
F
,
±
,
κ
n
H
∶
=
κ
r
→
|
V
n
H
.
Then the elements of matrix A have the following expressions:
(11)
a
m
n
F
F
=
∫
V
m
F
f
→
m
F
(
r
→
)
·
f
→
n
F
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
F
F
=
-
k
0
2
∫
V
m
F
d
v
f
→
m
F
(
r
→
)
·
∫
V
n
F
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
F
r
→
′
d
v
′
a
m
n
F
F
=
+
∫
V
m
F
d
v
∇
·
f
→
m
F
(
r
→
)
I
m
n
F
,
(
1
)
(
r
→
)
+
I
m
n
F
,
(
2
)
(
r
→
)
a
m
n
F
F
∶
=
a
m
n
F
F
,
1
-
k
0
2
a
m
n
F
F
,
2
+
a
m
n
F
F
,
3
,
(12)
a
m
n
H
F
=
∫
V
m
H
f
→
m
H
(
r
→
)
·
f
→
n
F
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
H
F
=
-
k
0
2
∫
V
m
H
d
v
f
→
m
H
(
r
→
)
·
∫
V
n
F
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
F
r
→
′
d
v
′
a
m
n
H
F
=
+
∫
V
m
H
d
v
∇
·
f
→
m
H
(
r
→
)
I
m
n
F
,
(
1
)
(
r
→
)
+
I
m
n
F
,
(
2
)
(
r
→
)
a
m
n
H
F
=
+
5
-
∫
[
V
n
H
]
d
s
I
m
n
F
,
(
1
)
(
r
→
)
+
I
m
n
F
,
(
2
)
(
r
→
)
a
m
n
H
F
∶
=
a
m
n
H
F
,
1
-
k
0
2
a
m
n
H
F
,
2
+
a
m
n
H
F
,
3
,
where
(13)
I
m
n
F
,
1
r
→
=
∫
V
n
F
G
r
→
,
r
→
′
κ
r
→
′
∇
′
·
f
→
n
F
r
→
′
d
v
′
,
I
m
n
F
,
2
r
→
=
∫
V
n
F
G
r
→
,
r
→
′
κ
n
F
,
-
-
κ
n
F
,
+
d
s
′
.
Note that a
m
n
F
F
,
1
=
0 in (11) when V
m
F and V
n
F are properly separated and that a
m
n
H
F
,
1
=
0 in (12) when V
m
H and V
n
F are properly separated:
(14)
a
m
n
F
H
=
∫
V
m
F
f
→
m
F
(
r
→
)
·
f
→
n
H
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
F
H
=
-
k
0
2
∫
V
m
F
d
v
f
→
m
F
(
r
→
)
·
∫
V
n
H
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
H
r
→
′
d
v
′
a
m
n
F
H
=
+
∫
V
m
F
d
v
∇
·
f
→
m
F
(
r
→
)
I
m
n
H
,
(
1
)
(
r
→
)
-
I
m
n
H
,
(
2
)
(
r
→
)
a
m
n
F
H
∶
=
a
m
n
F
H
,
1
-
k
0
2
a
m
n
F
H
,
2
+
a
m
n
F
H
,
3
,
(15)
a
m
n
H
H
=
∫
V
m
H
f
→
m
H
(
r
→
)
·
f
→
n
H
(
r
→
)
ε
r
(
r
→
)
d
v
a
m
n
H
H
=
-
k
0
2
∫
V
m
H
d
v
f
→
m
H
(
r
→
)
·
∫
V
n
H
G
r
→
,
r
→
′
κ
r
→
′
f
→
n
H
r
→
′
d
v
′
a
m
n
H
H
=
+
∫
V
m
H
d
v
∇
·
f
→
m
H
r
→
I
m
n
H
,
1
r
→
-
I
m
n
H
,
2
r
→
a
m
n
H
H
=
+
5
-
∫
V
m
H
d
s
∇
·
f
→
m
H
r
→
I
m
n
H
,
1
r
→
-
I
m
n
H
,
2
r
→
a
m
n
H
H
∶
=
a
m
n
H
H
,
1
-
k
0
2
a
m
n
H
H
,
2
+
a
m
n
H
H
,
3
,
where
(16)
I
m
n
H
,
1
r
→
=
∫
V
n
H
G
r
→
,
r
→
′
κ
r
→
′
∇
′
·
f
→
n
H
r
→
′
d
v
′
,
I
m
n
H
,
2
r
→
=
∫
V
n
H
d
s
′
G
r
→
,
r
→
′
κ
r
→
′
.
Note that a
m
n
F
H
,
1
=
0 in (14) when V
m
F and V
n
H are properly separated and that a
m
n
H
H
,
1
=
0 in (15) when V
m
H and V
n
H are properly separated.


## 2.3. The Frame for the VIE-FG-FFT

In this section, the FG-FFT technology is introduced into the VIE-MoM for both reducing memory requirement and improving computational efficiency.The entire MoM matrixA can be split into two parts: the near-field matrix A
near and the far-field matrix A
far, and
(17)
A
=
A
-
A
far
+
A
far
≃
A
near
+
A
far
,
where the A
near is a sparse matrix (the identification of a near-element will be presented in the second paragraph of Section 3), which is obtained by forcing all “far elements” of A
-
A
far to be equal to zero, and A
far can ultimately expressed in such a form as follows:
(18)
A
far
=
-
k
0
2
Π
→
f
·
G
Π
→
c
T
+
Π
f
G
Π
c
T
,
where Π
→
f, Π
→
c, Π
f, and Π
c are all sparse matrices which will be constructed in Section 2.4 and where the head mark “→” implies matrix elements being 3D vectors; G is a triple Toeplitz matrix related to Green’s function; the superscript T indicates matrix transpose.When an iterative solver is applied, the matrix-vector product will be performed by means of(19)
A
x
≃
A
near
x
+
A
far
x
,
where A
near
x is directly calculated, while A
far
x can be speeded up by means of the FFT through (18). In this way, the VIE-FG-FFT can reduce the memory requirement and the computational complexity to O
(
N
) and O
(
N
log
⁡
(
N
)
) theoretically (also see [9] for more detailed analysis).


## 2.4. Fitting Green’s Function

In this section, matricesΠ
→
f, Π
→
c, Π
f, and Π
c in (18) will be constructed. First, let a uniform Cartesian grid enclose the given volumetric region V. Use h
x, h
y, and h
z to denote the three grid spacing sizes in the directions x
^, y
^, and z
^, respectively. In this paper, however, the convention h
∶
=
h
x
=
h
y
=
h
z is always selected.In a uniform Cartesian grid, anexpansion box (or simplybox) C is defined as a cube-like collection composed of (
M
x
+
1
)
×
(
M
y
+
1
)
×
(
M
z
+
1
) nodes. When M
x
=
M
y
=
M
z
=
M, C includes M
+
1
3 nodes and M is called itsexpansion order (or simplyorder), written as |
C
|. Figure 2 illustrates a box of order 2.Figure 2
An expansion box of order 2.Assume thatC
n is a box centered at c
n of radius r
n and q
→ is a fixed point within C
n. Let Green’s function be represented into
(20)
G
r
→
,
q
→
=
∑
v
∈
C
n
‍
π
v
,
C
n
q
→
G
r
→
,
v
,
where the coefficients π
v
,
C
n
q
→ are to be determined and r
→ is an arbitrary point outside the box C
n, as shown in Figure 3.Figure 3
A 3D representation of the matching Green’s function.Now we select a spherical surface with the center atc
n and the radius R
n slightly larger than r
n, and select a group of sample points {
r
→
t
}
t
=
1
T on the surface [18], where the number T ought to be significantly greater than that of the nodes within C
n. In this paper, R
n
-
r
n
=
0.15 wavelength, which was adopted in [15]. Then we match Green’s function values at these sample points. The resultant systems of matrix equations can be solved by the least-square method [19]. In this paper, the order of a box is selected to be 2, and at this time the number of sample points on the testing spherical surface is selected to be 120.Without loss of generality, assume thatV
m
F
,
± and V
m
H (V
n
F
,
± and V
n
H) are in box C
m (C
n). Calculations of MoM matrix elements ultimately come down to calculations on the Gaussian points. We specify using N
G Gaussian points in a tetrahedral region and M
G Gaussian points in a triangular region. For convenience, the correspondences between integral regions and Gaussian points and weights that will be used are all listed in Table 1. When the testing element and the source element are separated by a proper distance, all the first terms in the right-hand side of (11), (12), (14), and (15) will vanish. At this case, we can easily obtain the following approximate expressions:
(21)
a
m
n
F
F
,
2
=
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
f
→
m
F
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
F
,
c
‍
‍
,
a
m
n
F
F
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
F
,
f
G
u
,
v
π
v
,
C
n
f
→
n
F
,
c
,
a
m
n
H
F
,
2
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
→
u
,
C
m
f
→
m
H
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
F
,
c
,
a
m
n
H
F
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
H
,
f
G
u
,
v
π
v
,
C
n
f
→
n
F
,
c
,
a
m
n
F
H
,
2
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
→
u
,
C
m
f
→
m
F
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
H
,
c
,
a
m
n
F
H
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
F
,
f
G
u
,
v
π
v
,
C
n
f
→
n
H
,
c
,
a
m
n
H
H
,
2
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
→
u
,
C
m
f
→
n
H
,
f
·
G
u
,
v
π
→
v
,
C
n
f
→
n
H
,
c
,
a
m
n
H
H
,
3
=
∑
u
∈
C
m
‍
∑
v
∈
C
n
‍
π
u
,
C
m
f
→
m
H
,
f
G
(
u
,
v
)
π
v
,
C
n
f
→
n
H
,
c
,
where the coefficients are calculated by using the following formulae:
(22)
π
→
u
,
C
m
f
→
m
F
,
f
=
∑
i
=
1
N
G
‍
w
i
f
→
m
F
,
+
p
→
i
+
π
u
,
C
m
p
→
i
+
+
f
→
m
F
,
-
p
→
i
-
π
u
,
C
m
p
→
i
-
,
π
u
,
C
m
f
→
m
F
,
f
=
∑
i
=
1
N
G
‍
w
i
∇
·
f
→
m
F
,
+
(
p
→
i
+
)
π
u
,
C
m
p
→
i
+
WWWWii
.
W
+
∇
·
f
→
m
F
,
-
p
→
i
-
π
u
,
C
m
p
→
i
-
,
π
→
v
,
C
n
f
→
n
F
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
F
,
+
f
→
n
F
,
+
q
→
j
+
π
v
,
C
n
q
→
j
+
WWWWii
.
W
+
κ
n
F
,
-
f
→
n
F
,
-
q
→
j
-
π
v
,
C
n
q
→
j
-
,
π
v
,
C
n
f
→
n
F
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
F
,
+
∇
·
f
→
n
F
,
+
q
→
j
+
π
v
,
C
n
q
→
j
+
WWWWWiW
+
κ
n
F
,
-
∇
·
f
→
n
F
,
-
q
→
j
-
π
v
,
C
n
q
→
j
-
π
v
,
C
n
f
→
n
F
,
c
=
+
∑
j
=
1
M
G
w
^
j
κ
n
F
,
-
-
κ
n
F
,
+
π
v
,
C
n
t
→
j
‍
,
π
→
u
,
C
m
f
→
m
H
,
f
=
∑
i
=
1
N
G
‍
w
i
f
→
m
H
p
→
i
π
u
,
C
m
p
→
i
,
π
u
,
C
m
f
→
m
H
,
f
=
∑
i
=
1
N
G
‍
w
i
∇
·
f
→
m
H
p
→
i
π
u
,
C
n
p
→
i
-
∑
i
=
1
M
G
‍
w
^
i
π
u
,
C
m
s
→
i
,
π
→
v
,
C
n
f
→
n
H
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
H
f
→
n
H
q
→
j
π
v
,
C
n
q
→
j
,
π
v
,
C
n
f
→
n
H
,
c
=
∑
j
=
1
N
G
‍
w
j
κ
n
H
∇
·
f
→
n
H
q
→
j
π
v
,
C
n
q
→
j
-
∑
j
=
1
M
G
‍
w
^
j
κ
n
H
π
v
,
C
n
t
→
j
.Table 1
Correspondences between integral regions and Gaussian points and weights.
Integral region
V
m
F
,
±
V
m
H
[
V
m
F
]
V
n
F
,
±
V
n
H
[
V
n
H
]
Gaussian points
{
p
→
i
±
}
i
=
1
N
G
{
p
→
i
}
i
=
1
N
G
{
s
→
i
}
i
=
1
M
G
{
q
→
j
±
}
j
=
1
N
G
{
q
→
j
}
j
=
1
N
G
{
t
→
j
}
j
=
1
M
G
Gaussian weights
{
w
i
}
i
=
1
N
G
{
w
^
i
}
i
=
1
M
G
{
w
j
}
j
=
1
N
G
{
w
^
j
}
j
=
1
M
GSubstituting (21) into (11)–(16), we can easily obtain formula (18). It should be pointed out that performing matrix-vector product A
far
x once based on (18) requires at least 8 FFTs.


## 2.5. On the Preprocessing Time

Here, the preprocessing time means the time taken by generating the coefficients in (21). Clearly, the preprocessing time is proportional to the corresponding computational complexity. For general cases, from the coefficient formulae in Section 2.4, it can be evaluated that the number N
FG-FFT of required float multiplications is about
(23)
N
FG-FFT
W
=
N
F
22
N
G
+
5
M
G
+
24
+
N
H
11
N
G
+
2
M
G
+
15
.In the VIE-P-FFT, the corresponding coefficients are generated by utilizing “projecting coefficients,” which requires solving the least-square problem with multiple right-hand terms. If the order of a box is selected to be 2, containing 27 grid points, and the number of the sample points on the testing spherical surface is selected to beT, then the matrix of the least-square problem in whether the VIE-P-FFT or the VIE-FG-FFT is a T
×
27 complex matrix. Further, if the SVD process of the matrix is ignored, then the least-square solution for a right-hand term requires about T
2
+
2
7
2 complex multiplications or about 3
(
T
2
+
2
7
2
) real multiplications because a complex multiplication requires at least 3 real multiplications. It is easily known from the projection scheme introduced in [20] that, for generating the corresponding coefficients, the number N
P-FFT of required float multiplications is about
(24)
N
P-FFT
=
24
N
T
2
+
2
7
2
,
which is usually much larger than N
FG-FFT.It can be concluded from the above analysis that the preprocessing time of the VIE-FG-FFT is obviously less than that of the VIE-P-FFT for inhomogeneous nonmagnetic dielectric objects.


## 3. Numerical Results

In this section, several examples are provided to demonstrate the validity, accuracy, and efficiency of the VIE-FG-FFT. In all the examples, the expansion orderM is always chosen as 2, and the Cartesian grid spacing sizes in different directions are always selected to be the same as each other; namely, h
∶
=
h
x
=
h
y
=
h
z. Besides, λ denotes the wavelength in free space. Our computing platform is a DELL T5400 workstation with 8 cores of clock frequency 3 GHz, and the FFT codes are from the FFTW [21].Assume that the testing function and the source function are located within boxesC
m and C
n, respectively. Then the corresponding matrix element is identified as anear element if and only if the distance between the center of C
m and that of C
n is smaller than the sum of the radius of C
m and that of C
n.


### 3.1. Example A: Homogeneous Ball

A dielectric spherical shell with the relative permittivity 2.0 is considered, as shown in Figure4. This object is discretized by using 23174 tetrahedrons with the average edge length 0.1
λ, producing 51231 unknowns. The Cartesian grid spacing size is chosen as h
=
0.1
λ.Figure 4
Bistatic RCS curves of the dielectric spherical shell in Example A.The bistatic RCS curve obtained by the FG-FFT is compared with the Mie series solution in Figure4. It can be seen that the two curves coincide very well, which demonstrates the validity of the FG-FFT.


### 3.2. Example B: Homogeneous Slab

A homogeneous hollow dielectric slab with the relative permittivity 2.5 is considered, as shown in Figure5. The direction of the incident wave is (
θ
in
,
φ
in
)
=
(
9
0
∘
,
4
5
∘
). This object is discretized by using 19301 tetrahedrons with the average edge length 0.1
λ, producing 42862 unknowns. The Cartesian grid spacing size is chosen as h
=
0.1
λ, 0.15
λ, and 0.2
λ, respectively.Figure 5
The dimensions of the hollow dielectric slab in Example B.The bistation RCS curves obtained by the direct MoM, the FG-FFT, and the other FFT-based methods (the P-FFT, the AIM, and the IE-FFT) are all shown in Figure6, and the RMSEs are recorded in Table 2, which are calculated by the following formula:
(25)
RMSE
=
1
N
∑
m
=
1
N
‍
RC
S
FFT-based
-
RC
S
MoM
2
,
where N is the number of sampling theta azimuth angles. In this example, N
=
100.Table 2
RMSEs of the RCSs in Example B.
Method
h
(
λ
)
RMSE
h
(
λ
)
RMSE
h
(
λ
)
RMSE
FG-FFT
0.1
0.0943
0.15
0.2004
0.2
0.5062
AIM
0.1
0.3895
0.15
0.9826
0.2
2.7141
IE-FFT
0.1
0.4415
0.15
1.0912
0.2
2.9753
P-FFT
0.1
0.0986
0.15
0.2048
0.2
0.5124Bistatic RCS curves of the homogeneous dielectric slab in Example B.
(a)
h
=
0.1
λ
(b)
h
=
0.15
λ
(c)
h
=
0.2
λIt can be seen from Figure6 and Table 2 that when the grid spacing size is h
=
0.1
λ, the four RCS curves obtained by the FG-FFT, P-FFT, IE-FFT, and AIM agree well with that by the direct MoM. However, when increasing the grid spacing size gradually, both the FG-FFT and P-FFT can keep consistent with the MoM, while both the P-FFT and IE-FFT do not.It can be concluded from the above experiments that compared with both the AIM and IE-FFT, both the FG-FFT and P-FFT are more accurate and not sensitive to the Cartesian grid spacing size.


### 3.3. Example C: Partitioned Homogeneous Slab

A partitioned homogeneous hollow dielectric slab with 4 different relative permittivities is considered, as shown in Figure7. The direction of the incident wave is (
θ
in
,
φ
in
)
=
(
9
0
∘
,
4
5
∘
). This object is discretized by using 17714 tetrahedrons with the average edge length 0.1
λ, producing 39927 unknowns. The Cartesian grid spacing size is chosen as h
=
0.1
λ, 0.15
λ, and 0.2
λ, respectively.Figure 7
The dimensions of the hollow dielectric slab in Example C.The bistatic RCS curves obtained by the direct MoM, the FG-FFT, and other FFT-based methods are all shown in Figure8, and the RMSEs are recorded in Table 3. In this example, N
=
100 for the RMSE.Table 3
RMSEs of the RCSs in Example C.
Method
h
(
λ
)
RMSE
h
(
λ
)
RMSE
h
(
λ
)
RMSE
FG-FFT
0.1
0.0704
0.15
0.1579
0.2
0.4854
AIM
0.1
0.3771
0.15
0.5648
0.2
1.2602
IE-FFT
0.1
0.4529
0.15
0.6406
0.2
1.3502
P-FFT
0.1
0.0738
0.15
0.1614
0.2
0.4936Bistatic RCS curves of the partitioned homogeneous dielectric slab in Example C.
(a)
h
=
0.1
λ
(b)
h
=
0.15
λ
(c)
h
=
0.2
λIt is again seen from Figure8 and Table 3 that the accuracy of both the FG-FFT and P-FFT is higher than that of both the IE-FFT and AIM.


### 3.4. Example D: Inhomogeneous Slab

Here the dielectric slab in Figure5 is again considered. The direction of the incident wave is kept unchanged. But at this time, the relative permittivity
(26)
ε
r
(
x
,
y
,
z
)
=
ε
r
′
(
x
,
y
,
z
)
-
j
ε
r
′′
(
x
,
y
,
z
)
is defined by
(27)
ε
r
′
(
x
,
y
,
z
)
=
0.2
x
+
y
+
z
10.1
+
2.0
,
ε
r
′′
(
x
,
y
,
z
)
=
0.1
,
where x
,
y
,
z
≥
0. Now, the average edge length used in the discretization is selected within {
0.08
λ
,
0.1
λ
,
0.15
λ
}.Different discretization granularities correspond to different numbers of unknowns, as shown in Table4. The bistatic RCS curves obtained by the FG-FFT and P-FFT for different discretization granularities are all shown in Figure 9, and the preprocessing time for the FG-FFT and P-FFT is recorded in Table 4.Table 4
The Preprocessing time in Example D.
Method
Edge length (λ)
Unknowns
Preprocessing time (s)
FG-FFT
0.15
10875
18.42
P-FFT
84.36
FG-FFT
0.1
39927
261.27
P-FFT
512.69
FG-FFT
0.08
61376
934.07
P-FFT
1326.58Figure 9
Bistatic RCS curves of the inhomogeneous dielectric slab (h
=
0.2
λ) in Example D.It can be seen from Figure9 that the RCS curves for different discretization granularities are almost the same as each other. It can be known from Table 4 that the preprocessing time of the FG-FFT is obviously less than that of the P-FFT for different discretization granularities.


## 3.1. Example A: Homogeneous Ball

A dielectric spherical shell with the relative permittivity 2.0 is considered, as shown in Figure4. This object is discretized by using 23174 tetrahedrons with the average edge length 0.1
λ, producing 51231 unknowns. The Cartesian grid spacing size is chosen as h
=
0.1
λ.Figure 4
Bistatic RCS curves of the dielectric spherical shell in Example A.The bistatic RCS curve obtained by the FG-FFT is compared with the Mie series solution in Figure4. It can be seen that the two curves coincide very well, which demonstrates the validity of the FG-FFT.


## 3.2. Example B: Homogeneous Slab

A homogeneous hollow dielectric slab with the relative permittivity 2.5 is considered, as shown in Figure5. The direction of the incident wave is (
θ
in
,
φ
in
)
=
(
9
0
∘
,
4
5
∘
). This object is discretized by using 19301 tetrahedrons with the average edge length 0.1
λ, producing 42862 unknowns. The Cartesian grid spacing size is chosen as h
=
0.1
λ, 0.15
λ, and 0.2
λ, respectively.Figure 5
The dimensions of the hollow dielectric slab in Example B.The bistation RCS curves obtained by the direct MoM, the FG-FFT, and the other FFT-based methods (the P-FFT, the AIM, and the IE-FFT) are all shown in Figure6, and the RMSEs are recorded in Table 2, which are calculated by the following formula:
(25)
RMSE
=
1
N
∑
m
=
1
N
‍
RC
S
FFT-based
-
RC
S
MoM
2
,
where N is the number of sampling theta azimuth angles. In this example, N
=
100.Table 2
RMSEs of the RCSs in Example B.
Method
h
(
λ
)
RMSE
h
(
λ
)
RMSE
h
(
λ
)
RMSE
FG-FFT
0.1
0.0943
0.15
0.2004
0.2
0.5062
AIM
0.1
0.3895
0.15
0.9826
0.2
2.7141
IE-FFT
0.1
0.4415
0.15
1.0912
0.2
2.9753
P-FFT
0.1
0.0986
0.15
0.2048
0.2
0.5124Bistatic RCS curves of the homogeneous dielectric slab in Example B.
(a)
h
=
0.1
λ
(b)
h
=
0.15
λ
(c)
h
=
0.2
λIt can be seen from Figure6 and Table 2 that when the grid spacing size is h
=
0.1
λ, the four RCS curves obtained by the FG-FFT, P-FFT, IE-FFT, and AIM agree well with that by the direct MoM. However, when increasing the grid spacing size gradually, both the FG-FFT and P-FFT can keep consistent with the MoM, while both the P-FFT and IE-FFT do not.It can be concluded from the above experiments that compared with both the AIM and IE-FFT, both the FG-FFT and P-FFT are more accurate and not sensitive to the Cartesian grid spacing size.


## 3.3. Example C: Partitioned Homogeneous Slab

A partitioned homogeneous hollow dielectric slab with 4 different relative permittivities is considered, as shown in Figure7. The direction of the incident wave is (
θ
in
,
φ
in
)
=
(
9
0
∘
,
4
5
∘
). This object is discretized by using 17714 tetrahedrons with the average edge length 0.1
λ, producing 39927 unknowns. The Cartesian grid spacing size is chosen as h
=
0.1
λ, 0.15
λ, and 0.2
λ, respectively.Figure 7
The dimensions of the hollow dielectric slab in Example C.The bistatic RCS curves obtained by the direct MoM, the FG-FFT, and other FFT-based methods are all shown in Figure8, and the RMSEs are recorded in Table 3. In this example, N
=
100 for the RMSE.Table 3
RMSEs of the RCSs in Example C.
Method
h
(
λ
)
RMSE
h
(
λ
)
RMSE
h
(
λ
)
RMSE
FG-FFT
0.1
0.0704
0.15
0.1579
0.2
0.4854
AIM
0.1
0.3771
0.15
0.5648
0.2
1.2602
IE-FFT
0.1
0.4529
0.15
0.6406
0.2
1.3502
P-FFT
0.1
0.0738
0.15
0.1614
0.2
0.4936Bistatic RCS curves of the partitioned homogeneous dielectric slab in Example C.
(a)
h
=
0.1
λ
(b)
h
=
0.15
λ
(c)
h
=
0.2
λIt is again seen from Figure8 and Table 3 that the accuracy of both the FG-FFT and P-FFT is higher than that of both the IE-FFT and AIM.


## 3.4. Example D: Inhomogeneous Slab

Here the dielectric slab in Figure5 is again considered. The direction of the incident wave is kept unchanged. But at this time, the relative permittivity
(26)
ε
r
(
x
,
y
,
z
)
=
ε
r
′
(
x
,
y
,
z
)
-
j
ε
r
′′
(
x
,
y
,
z
)
is defined by
(27)
ε
r
′
(
x
,
y
,
z
)
=
0.2
x
+
y
+
z
10.1
+
2.0
,
ε
r
′′
(
x
,
y
,
z
)
=
0.1
,
where x
,
y
,
z
≥
0. Now, the average edge length used in the discretization is selected within {
0.08
λ
,
0.1
λ
,
0.15
λ
}.Different discretization granularities correspond to different numbers of unknowns, as shown in Table4. The bistatic RCS curves obtained by the FG-FFT and P-FFT for different discretization granularities are all shown in Figure 9, and the preprocessing time for the FG-FFT and P-FFT is recorded in Table 4.Table 4
The Preprocessing time in Example D.
Method
Edge length (λ)
Unknowns
Preprocessing time (s)
FG-FFT
0.15
10875
18.42
P-FFT
84.36
FG-FFT
0.1
39927
261.27
P-FFT
512.69
FG-FFT
0.08
61376
934.07
P-FFT
1326.58Figure 9
Bistatic RCS curves of the inhomogeneous dielectric slab (h
=
0.2
λ) in Example D.It can be seen from Figure9 that the RCS curves for different discretization granularities are almost the same as each other. It can be known from Table 4 that the preprocessing time of the FG-FFT is obviously less than that of the P-FFT for different discretization granularities.


## 4. Conclusions

In this paper, a new realization of the VIE combined with the fast Fourier transform (VIE-FG-FFT) has been established to solve EM scattering from inhomogeneous nonmagnetic dielectric objects. The proposed method has been compared with several existing popular FFT-based methods, including the VIE-P-FFT, VIE-AIM, and VIE-IE-FFT. The accuracy of the VIE-FG-FFT is almost the same as that of the VIE-P-FFT and higher than that of the VIE-AIM and VIE-IE-FFT especially when increasing the grid spacing size. Besides, the preprocessing time of the VIE-FG-FFT is obviously less than that of the P-FFT for inhomogeneous nonmagnetic dielectric objects.



---

*Source: 720947-2014-12-24.xml*