# Fitting Green’s Function FFT Acceleration Applied to Anisotropic Dielectric Scattering Problems

Authors: Shu-Wen Chen; Feng Lu; Yao Ma
Journal: International Journal of Antennas and Propagation
 (2015)

---

## Abstract

A volume integral equation based fast algorithm using the Fast Fourier Transform of fitting Green’s function (FG-FFT) is proposed in this paper for analysis of electromagnetic scattering from 3D anisotropic dielectric objects. For the anisotropic VIE model, geometric discretization is still implemented by tetrahedron cells and the Schaubert-Wilton-Glisson (SWG) basis functions are also used to represent the electric flux density vectors. Compared with other Fast Fourier Transform based fast methods, using fitting Green’s function technique has higher accuracy and can be applied to a relatively coarse grid, so the Fast Fourier Transform of fitting Green’s function is selected to accelerate anisotropic dielectric model of volume integral equation for solving electromagnetic scattering problems. Besides, the near-field matrix elements in this method are used to construct preconditioner, which has been proved to be effective. At last, several representative numerical experiments proved the validity and efficiency of the proposed method.


---

## Body

## 1. Introduction

The volume integral equation (VIE) model using the SWG basis function [1], based on the method of moments (MoM) [2], is one of the efficient methods to analyze the problem of electromagnetic scattering from 3D inhomogeneous anisotropic dielectric objects. Generally, for the conventional VIE model applied to anisotropic objects, both the storage requirement and the computational complexity of matrix-vector multiplication are proportional to O
(
N
2
), where N denotes the number of unknowns. Therefore, the conventional VIE-MoM has a huge amount of unknowns and a dense impedance matrix, which is not suitable for the direct analysis of electrically large anisotropic dielectric objects.One of the approaches for accelerating the VIE model is the MLFMA, such as [3–5]. Another approach for improving efficiency is the VIE-based equivalent dipole method [6]. The third approach for improving efficiency is the VIE in combination with the Fast Fourier Transform (FFT). At present, alternative FFT-based accelerating methods include the VIE-AIM [7], the VIE-P-FFT [8, 9], VIE-IE-FFT [10], and VIE-FG-FFT [11], which are simply called the FFT-based methods. The thought behind the IE-FFT is interpolating Green’s function onto the nodes of a uniform Cartesian grid by Lagrange polynomials. The main idea of Adaptive Integral Method (AIM) is to project the local basis functions onto node basis functions defined on the nodes of a uniform Cartesian grid. Its projection is obtained by matching the high order multipole moments, so it has no relation with Green’s function. In the P-FFT method, the projection is obtained by matching the electromagnetic potentials rather than the high order multipole moments, having higher accuracy compared with the AIM and IE-FFT, and can be applied to a relatively coarse grid, reducing the burden on the FFT.Until now, we have not seen these FFT methods applied to anisotropic VIE model. Through the previous studies of these fast algorithms, it has been known that the VIE-P-FFT and VIE-FG-FFT have the same higher accuracy than the VIE-AIM and VIE-IE-FFT, especially when increasing the Cartesian grid spacing size. Besides, from the study in [11], the VIE-FG-FFT preprocessing time is generally less than that of the VIE-P-FFT, which depends on inhomogeneous degree of the dielectric object. So, for general inhomogeneous anisotropic objects, the scheme of the FG-FFT is selected to do algorithm acceleration in this work.Briefly speaking, advantages of the FG-FFT can be generalized as requiring only Green’s function to participate in fitting procedure for calculating equivalent projection coefficients on Cartesian grids. The scheme is simple, totally independent of selected basis functions and dielectric permittivity tensor. As a result, the FG-FFT is more flexible to complicated dielectric structure (including anisotropic media). Afterwards, constructing preconditioner using near-field elements is investigated. The preconditioner can effectively speed up convergence of the CG solution. Several numerical experiments demonstrate the good accelerating effectiveness.In this paper, a new FG-FFT acceleration applied to the anisotropic VIE model is developed. Actually, it is a later extended work of the FG-FFT technology. This new realization can not only achieve high accuracy and reduce memory requirement but also have better adaptability to complicated dielectric materials (including anisotropic media). This paper is organized as follows: In Section2, the anisotropic VIE-FG-FFT frame and preconditioning technique are presented. In Section 3, some numerical examples are provided to prove the validity and efficiency of the proposed method. Finally, the conclusion is given in Section 4.


## 2. Formulation

### 2.1. The Volume Integral Equation for Anisotropic Media

The free space has the permittivityε
0 and permeability μ
0, respectively. Let V denote an inhomogeneous anisotropic dielectric volume with relative permittivity tensor ε
̿
r and relative permeability tensor μ
̿
r
=
I
̿. Here, I
̿ denotes unit tensor. Let E
→
i be the incident electric field and E
→
s the scattered electric field; then the total electric field E
→
t
o
t can be expressed as(1)
E
→
tot
=
E
→
i
+
E
→
s
.The scattered field due to equivalent volume polarization currentJ
→
V
(
r
→
) can be written as(2)
E
→
s
r
→
=
-
j
ω
μ
0
I
̿
+
1
k
0
2
∇
∇
·
∫
V
g
r
→
,
r
→
′
J
→
V
r
→
′
d
v
′
,where (3)
J
→
V
r
→
=
j
ω
κ
̿
r
→
D
→
tot
r
→
,
k
0
=
ω
μ
0
ε
0
,
η
0
=
μ
0
ε
0
.
κ
̿ is permittivity contrast ratio tensor, defined as (4)
κ
̿
r
→
=
I
̿
-
ε
̿
r
-
1
.


### 2.2. VIE Model for Anisotropic Media

The electric flux densityD
→
t
o
t can be selected as the unknown function. After V is discretized by tetrahedrons, D
→
t
o
t
/
ε
0 can be expanded with the SWG basis functions (5)
D
→
tot
r
→
ε
0
=
∑
n
=
1
N
x
n
f
→
n
r
→
,where “N” means the number of the SWG functions. The details about the SWG function are referred to in [1]. Use the Galerkin test procedure; then the corresponding MoM matrix equation can be built below: (6)
A
x
=
B
,where A
=
(
a
m
n
)
N
×
N is a N
×
N matrix; x
=
[
x
1
,
x
2
,
…
,
x
N
]
T is the unknown vector; B
=
[
b
1
,
b
2
,
…
,
b
N
]
T is the incident vector.More specifically, the matrix elementa
m
n gotten from f
→
m testing the basis function f
→
n is (7)
a
m
n
=
∫
V
m
f
→
m
r
→
·
ε
̿
-
1
r
→
·
f
→
n
r
→
-
k
0
2
∫
V
m
f
→
m
r
→
d
v
·
∫
V
n
g
r
→
,
r
→
′
κ
̿
r
→
′
·
f
→
n
r
→
′
d
v
′
-
∫
V
m
f
→
m
r
→
d
v
·
∇
∇
·
∫
V
n
g
r
→
,
r
→
′
κ
̿
r
→
′
·
f
→
n
r
→
′
d
v
′
,
(8)
b
n
=
∫
V
n
f
→
n
r
→
·
E
→
i
r
→
d
v
.After expanding (7), a
m
n can be separated into the algebraic sum of 6 terms(9)
a
m
n
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
+
I
m
n
5
-
I
m
n
6
,where(10)
I
m
n
1
=
∫
V
m
f
→
m
r
→
·
ε
̿
r
-
1
r
→
f
→
n
r
→
d
v
,
(11)
I
m
n
2
=
k
0
2
∫
V
m
f
→
m
r
→
d
v
·
∫
V
n
g
r
→
,
r
→
′
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
,
(12)
I
m
n
3
=
∫
V
m
∇
·
f
→
m
r
→
d
v
∫
∂
V
n
d
s
′
g
r
→
,
r
→
′
m
^
V
n
·
κ
̿
r
→
′
f
→
n
r
→
′
,
(13)
I
m
n
4
=
∫
V
m
∇
·
f
→
m
r
→
d
v
∫
V
n
g
r
→
,
r
→
′
∇
′
·
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
,
(14)
I
m
n
5
=
∫
∂
V
m
d
s
m
^
V
m
·
f
→
m
r
→
∫
∂
V
n
d
s
′
g
r
→
,
r
→
′
m
^
V
n
·
κ
̿
r
→
′
f
→
n
r
→
′
,
(15)
I
m
n
6
=
∫
∂
V
m
d
s
m
^
V
m
·
f
→
m
r
→
∫
V
n
g
r
→
,
r
→
′
∇
′
·
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
,where ∂
V
m
∂
V
n represents the outer boundary surface of V
m
V
n and m
^
V
m
V
n is outer normal vector of the volume element V
m
V
n. The element V
m
V
n can represent any full or half SWG basis function, shown in Figure 1. It is noted that the above 6 terms do not always exist in a matrix element. On one hand, when f
→
n is a full SWG, the inner integral of (12) and (14) is zero for ordinary VIE model; now it is no longer zero for anisotropic VIE model. On the other hand, when f
→
n is a half SWG, the inner integral of (12) and (14) includes one triangular surface integral in ordinary VIE model; however the inner integral becomes the sum of four triangular surface integrals in anisotropic VIE model. If the superscripts “F” and “H” mean “full SWG function” and “half SWG function,” respectively, according to different combinations of full and half SWG basis functions, there exist four types of matrix elements for anisotropic VIE model:(16)
a
m
n
F
F
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
,
a
m
n
F
H
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
,
a
m
n
H
F
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
+
I
m
n
5
-
I
m
n
6
,
a
m
n
H
H
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
+
I
m
n
5
-
I
m
n
6
.For anisotropic media, matrix element expression a
m
n
F
F adds one term -
I
m
n
3, while a
m
n
H
F adds two terms -
I
m
n
3 and I
m
n
5, compared with that in [11]. No matter which type of matrix elements, the value of I
m
n
(
1
) has the following basic law: (17)
I
m
n
1
=
0
,
when the basis function and test function do not overlap
≠
0
,
when the basis function and test function overlap or partly overlap.Figure 1
Full SWG and half SWG.


### 2.3. The FG-FFT Frame for Anisotropic Media and Fitting Green’s Function

The FG-FFT frame for the anisotropic media MoM model is similar to that in ordinary media [11]. The entire impedance matrix A can also be separated into near-field matrix A
n
e
a
r and far-field matrix A
f
a
r; thus the matrix-vector product of an iterative solver will be performed by(18)
A
x
≈
A
near
x
+
A
far
x
.Similarly, near-field matrix A
n
e
a
r is obtained by -
A
f
a
r; A
n
e
a
r
x could be calculated directly. However, A
f
a
r
x could be speeded up by means of the FFT, which can reduce the memory requirement and computational complexity to O
(
N
) and O
(
N
log
⁡
(
N
)
), theoretically.Let a uniform Cartesian grid enclose the given domainV. Three grid spacing sizes in the directions x
^, y
^, and z
^ are h
x, h
y, and h
z, respectively. For simplicity, define h
≔
h
x
=
h
y
=
h
z. In the Cartesian grid, an expansion box C
n is defined as a cube-like collection composed of M
x
+
1
×
M
y
+
1
×
M
z
+
1 grid nodes. When M
=
M
x
=
M
y
=
M
z, C
n includes M
+
1
3 nodes, and M is called its expansion order. Figure 1 describes a box of order 2. Assume that a box C
n has center at c
n and radius r
n; Green’s function from source q
→ can be represented into(19)
G
r
→
,
q
→
=
∑
u
∈
C
n
π
u
,
C
n
q
→
G
r
→
,
u
→
,where the coefficients π
u
,
C
n
q
→ are to be determined and r
→ is an arbitrary point outside the box C
n, as shown in Figure 2.Figure 2
A 3D representation of the matching Green’s function values.Through selecting a group of sample pointsr
→
t
t
=
1
T on a spherical surface [12] with a bit larger radius R
n, the coefficients in (19) can be determined. When the number of sample points T is significantly greater than that of nodes in C
n, through matching Green’s function values at these sample points, the resultant systems of matrix equations can be solved by the least-square method. Without loss of generality, assume that integral domain V
m
F
,
± and V
m
H (V
n
F
,
± and V
n
H) are in the box C
m (C
n). Calculations of MoM matrix elements ultimately come down to calculating projections on these Gaussian points. Due to different inner and outer integrals, for anisotropic model, the 6 types of projection coefficients on the grid nodes are deduced as follows: (20)
T
→
V
n
F
,
±
V
n
H
1
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
κ
̿
p
→
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
κ
̿
p
→
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
→
u
,
C
n
T
V
n
F
,
±
V
n
H
1
,
T
→
V
n
F
,
±
V
n
H
2
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
→
u
,
C
n
T
V
n
F
,
±
V
n
H
2
,
T
V
n
F
,
±
V
n
H
3
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
∇
′
·
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
∇
·
κ
̿
p
→
i
f
→
n
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
∇
·
κ
-
-
p
→
i
f
→
n
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
V
n
F
,
±
V
n
H
3
,
T
V
n
F
,
±
V
n
H
4
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
∇
′
·
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
∇
·
f
→
n
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
∇
·
f
→
n
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
V
n
F
,
±
V
n
H
4
,
T
∂
V
n
F
,
±
∂
V
n
H
5
=
∫
∂
V
n
F
,
±
∂
V
n
H
d
s
′
g
r
→
,
r
→
′
m
^
V
i
·
κ
̿
r
→
′
f
→
i
r
→
′
=
∑
s
=
1
4
∑
i
=
1
M
G
V
n
∑
u
=
1
M
C
n
w
i
′
m
^
V
i
·
κ
̿
s
→
i
f
→
n
s
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
s
=
1
4
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
M
G
V
n
s
w
i
′
m
^
V
i
·
κ
̿
s
→
i
f
→
n
s
→
i
π
u
,
C
n
p
→
i
=
∑
s
=
1
4
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
S
,
T
∂
V
n
F
,
±
∂
V
n
H
5
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
∂
V
n
F
,
±
∂
V
n
H
5
,
T
∂
V
n
H
6
=
∫
∂
V
n
H
d
s
′
g
r
→
,
r
→
′
m
^
V
i
·
f
→
i
r
→
′
=
∑
i
=
1
M
G
V
n
∑
u
=
1
M
C
n
w
i
′
m
^
V
i
·
f
→
n
s
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
M
G
V
n
w
i
′
m
^
V
i
·
f
→
n
s
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
∂
V
n
H
6
.In the above 6 formulas, there are three different integral regions: tetrahedral integral region V
m
(
n
)
F
,
±, V
m
(
n
)
H, and triangular integral region V
m
(
n
)
F
H
s, corresponding to their Gaussian points p
→
i
i
=
1
N
G
F
,
±, p
→
i
i
=
1
N
G
H, s
→
i
i
=
1
M
G and weights w
i
i
=
1
N
G
F
,
±, w
i
i
=
1
N
G
H, w
i
′
i
=
1
M
G, respectively. Besides difference between κ
-
-
(
r
→
) and κ
(
r
→
), note that T
(
5
)contains four triangular surface integrals for anisotropic model, while in [11] for isotropic model, T
(
5
)contains only one triangular surface integral. Obviously, (21)
T
→
V
n
F
V
n
H
1
=
T
→
V
n
F
,
+
1
+
T
→
V
n
F
,
-
1
T
→
V
n
H
1
,
T
→
V
n
F
V
n
H
2
=
T
→
V
n
F
,
+
2
+
T
→
V
n
F
,
-
2
T
→
V
n
H
2
,
T
V
n
F
V
n
H
3
=
T
V
n
F
,
+
3
+
T
V
n
F
,
-
3
T
V
n
H
3
,
T
V
n
F
V
n
H
4
=
T
V
n
F
,
+
4
+
T
V
n
F
,
-
4
T
V
n
H
4
,
T
V
n
F
V
n
H
5
=
T
V
n
F
,
+
5
+
T
V
n
F
,
-
5
T
V
n
H
5
,
T
∂
V
n
H
6
=
T
∂
V
n
H
6
Only about half SWG
.So, the ultimate π
→ or π matrix of the FFT uplink and downlink can be gotten:(22)
π
→
u
,
C
n
T
V
n
F
V
n
H
1
=
π
→
u
,
C
n
T
V
n
F
,
+
V
n
H
1
+
π
→
u
,
C
n
T
V
n
F
,
-
V
n
H
1
π
→
u
,
C
n
T
V
n
H
1
,
π
→
u
,
C
n
T
V
n
F
V
n
H
2
=
π
→
u
,
C
n
T
V
n
F
,
+
V
n
H
2
+
π
→
u
,
C
n
T
V
n
F
,
-
V
n
H
2
π
→
u
,
C
n
T
V
n
H
2
,
π
u
,
C
n
T
V
n
F
V
n
H
3
=
π
u
,
C
n
T
V
n
F
,
+
V
n
H
3
+
π
u
,
C
n
T
V
n
F
,
-
V
n
H
3
π
u
,
C
n
T
V
n
H
3
,
π
u
,
C
n
T
V
n
F
V
n
H
4
=
π
u
,
C
n
T
V
n
F
,
+
V
n
H
4
+
π
u
,
C
n
T
V
n
F
,
-
V
n
H
4
π
u
,
C
n
T
V
n
H
4
,
π
u
,
C
n
T
V
n
F
V
n
H
5
=
π
u
,
C
n
T
V
n
F
,
+
V
n
H
5
+
π
u
,
C
n
T
V
n
F
,
-
V
n
H
5
π
u
,
C
n
T
V
n
H
5
,
π
u
,
C
n
T
∂
V
n
H
6
=
π
u
,
C
n
T
∂
V
n
H
6
Only about half SWG
.Subscript u
=
1,2
,
…
,
M
C
n. At last, the far-field matrix elements can be expressed in such following form:(23)
a
m
n
F
F
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
F
2
G
u
,
v
π
→
v
,
C
n
T
V
n
F
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
F
3
,
a
m
n
F
H
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
F
2
G
u
,
v
π
→
v
,
C
n
T
V
n
H
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
H
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
H
3
,
a
m
n
H
F
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
H
2
G
u
,
v
π
→
v
,
C
n
T
V
n
F
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
F
3
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
F
5
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
F
3
,
a
m
n
H
H
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
H
2
G
u
,
v
π
→
v
,
C
n
T
V
n
H
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
H
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
H
3
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
H
5
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
H
3
.Compared with ordinary media, a
m
n
F
F adds one term -
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5, while a
m
n
H
F adds two terms -
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5 and ∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
F
5. So, it is more complicated relatively.


### 2.4. Near-Field Preconditioning Technique

There are a few preconditioning techniques successfully developed for some fast algorithms, solving MoM linear systems arising from surface integral equation (SIE). However, the preconditioning techniques applied to FFT-based VIE are less studied. This section investigates convergence effect of near-field preconditioning technique for the anisotropic VIE-FG-FFT. In the FFT-based methods, only near-field matrix elements are stored in the memory, and because of the characteristics of 3D Green’s function, near-field elements are generally larger than far-field matrix elements, so near-field matrix can be seen as sparsification of the global matrix by abandoning small matrix elements. The near-field matrix could be used to construct DILU preconditioner [13]. However, different from the MLFMA, the near-field matrix elements in the FG-FFT are not real; namely, the near-field matrix is built with FFT revision. Experiments demonstrate it is not right to directly use the revised near-field matrix to construct preconditioner. Now we use real near-field matrix A
n
e
a
r, which can be separated into three parts:(24)
A
near
=
L
near
+
D
near
+
U
near
,where L
n
e
a
r is lower triangular part of A
n
e
a
r; U
n
e
a
r is upper triangular part of A
n
e
a
r; D
n
e
a
r is main diagonal part of A
n
e
a
r. Hence, the new preconditioning matrix is(25)
M
=
D
near
+
L
near
D
near
-
1
D
near
+
U
near
.Namely, (26)
M
=
D
near
+
L
near
I
+
D
near
-
1
U
near
.This preconditioner is easy to implement and additional time cost of each CG iteration is extremely small. Numerical examples have proved the good effectiveness in Section 3.


## 2.1. The Volume Integral Equation for Anisotropic Media

The free space has the permittivityε
0 and permeability μ
0, respectively. Let V denote an inhomogeneous anisotropic dielectric volume with relative permittivity tensor ε
̿
r and relative permeability tensor μ
̿
r
=
I
̿. Here, I
̿ denotes unit tensor. Let E
→
i be the incident electric field and E
→
s the scattered electric field; then the total electric field E
→
t
o
t can be expressed as(1)
E
→
tot
=
E
→
i
+
E
→
s
.The scattered field due to equivalent volume polarization currentJ
→
V
(
r
→
) can be written as(2)
E
→
s
r
→
=
-
j
ω
μ
0
I
̿
+
1
k
0
2
∇
∇
·
∫
V
g
r
→
,
r
→
′
J
→
V
r
→
′
d
v
′
,where (3)
J
→
V
r
→
=
j
ω
κ
̿
r
→
D
→
tot
r
→
,
k
0
=
ω
μ
0
ε
0
,
η
0
=
μ
0
ε
0
.
κ
̿ is permittivity contrast ratio tensor, defined as (4)
κ
̿
r
→
=
I
̿
-
ε
̿
r
-
1
.


## 2.2. VIE Model for Anisotropic Media

The electric flux densityD
→
t
o
t can be selected as the unknown function. After V is discretized by tetrahedrons, D
→
t
o
t
/
ε
0 can be expanded with the SWG basis functions (5)
D
→
tot
r
→
ε
0
=
∑
n
=
1
N
x
n
f
→
n
r
→
,where “N” means the number of the SWG functions. The details about the SWG function are referred to in [1]. Use the Galerkin test procedure; then the corresponding MoM matrix equation can be built below: (6)
A
x
=
B
,where A
=
(
a
m
n
)
N
×
N is a N
×
N matrix; x
=
[
x
1
,
x
2
,
…
,
x
N
]
T is the unknown vector; B
=
[
b
1
,
b
2
,
…
,
b
N
]
T is the incident vector.More specifically, the matrix elementa
m
n gotten from f
→
m testing the basis function f
→
n is (7)
a
m
n
=
∫
V
m
f
→
m
r
→
·
ε
̿
-
1
r
→
·
f
→
n
r
→
-
k
0
2
∫
V
m
f
→
m
r
→
d
v
·
∫
V
n
g
r
→
,
r
→
′
κ
̿
r
→
′
·
f
→
n
r
→
′
d
v
′
-
∫
V
m
f
→
m
r
→
d
v
·
∇
∇
·
∫
V
n
g
r
→
,
r
→
′
κ
̿
r
→
′
·
f
→
n
r
→
′
d
v
′
,
(8)
b
n
=
∫
V
n
f
→
n
r
→
·
E
→
i
r
→
d
v
.After expanding (7), a
m
n can be separated into the algebraic sum of 6 terms(9)
a
m
n
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
+
I
m
n
5
-
I
m
n
6
,where(10)
I
m
n
1
=
∫
V
m
f
→
m
r
→
·
ε
̿
r
-
1
r
→
f
→
n
r
→
d
v
,
(11)
I
m
n
2
=
k
0
2
∫
V
m
f
→
m
r
→
d
v
·
∫
V
n
g
r
→
,
r
→
′
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
,
(12)
I
m
n
3
=
∫
V
m
∇
·
f
→
m
r
→
d
v
∫
∂
V
n
d
s
′
g
r
→
,
r
→
′
m
^
V
n
·
κ
̿
r
→
′
f
→
n
r
→
′
,
(13)
I
m
n
4
=
∫
V
m
∇
·
f
→
m
r
→
d
v
∫
V
n
g
r
→
,
r
→
′
∇
′
·
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
,
(14)
I
m
n
5
=
∫
∂
V
m
d
s
m
^
V
m
·
f
→
m
r
→
∫
∂
V
n
d
s
′
g
r
→
,
r
→
′
m
^
V
n
·
κ
̿
r
→
′
f
→
n
r
→
′
,
(15)
I
m
n
6
=
∫
∂
V
m
d
s
m
^
V
m
·
f
→
m
r
→
∫
V
n
g
r
→
,
r
→
′
∇
′
·
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
,where ∂
V
m
∂
V
n represents the outer boundary surface of V
m
V
n and m
^
V
m
V
n is outer normal vector of the volume element V
m
V
n. The element V
m
V
n can represent any full or half SWG basis function, shown in Figure 1. It is noted that the above 6 terms do not always exist in a matrix element. On one hand, when f
→
n is a full SWG, the inner integral of (12) and (14) is zero for ordinary VIE model; now it is no longer zero for anisotropic VIE model. On the other hand, when f
→
n is a half SWG, the inner integral of (12) and (14) includes one triangular surface integral in ordinary VIE model; however the inner integral becomes the sum of four triangular surface integrals in anisotropic VIE model. If the superscripts “F” and “H” mean “full SWG function” and “half SWG function,” respectively, according to different combinations of full and half SWG basis functions, there exist four types of matrix elements for anisotropic VIE model:(16)
a
m
n
F
F
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
,
a
m
n
F
H
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
,
a
m
n
H
F
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
+
I
m
n
5
-
I
m
n
6
,
a
m
n
H
H
=
I
m
n
1
-
I
m
n
2
-
I
m
n
3
+
I
m
n
4
+
I
m
n
5
-
I
m
n
6
.For anisotropic media, matrix element expression a
m
n
F
F adds one term -
I
m
n
3, while a
m
n
H
F adds two terms -
I
m
n
3 and I
m
n
5, compared with that in [11]. No matter which type of matrix elements, the value of I
m
n
(
1
) has the following basic law: (17)
I
m
n
1
=
0
,
when the basis function and test function do not overlap
≠
0
,
when the basis function and test function overlap or partly overlap.Figure 1
Full SWG and half SWG.


## 2.3. The FG-FFT Frame for Anisotropic Media and Fitting Green’s Function

The FG-FFT frame for the anisotropic media MoM model is similar to that in ordinary media [11]. The entire impedance matrix A can also be separated into near-field matrix A
n
e
a
r and far-field matrix A
f
a
r; thus the matrix-vector product of an iterative solver will be performed by(18)
A
x
≈
A
near
x
+
A
far
x
.Similarly, near-field matrix A
n
e
a
r is obtained by -
A
f
a
r; A
n
e
a
r
x could be calculated directly. However, A
f
a
r
x could be speeded up by means of the FFT, which can reduce the memory requirement and computational complexity to O
(
N
) and O
(
N
log
⁡
(
N
)
), theoretically.Let a uniform Cartesian grid enclose the given domainV. Three grid spacing sizes in the directions x
^, y
^, and z
^ are h
x, h
y, and h
z, respectively. For simplicity, define h
≔
h
x
=
h
y
=
h
z. In the Cartesian grid, an expansion box C
n is defined as a cube-like collection composed of M
x
+
1
×
M
y
+
1
×
M
z
+
1 grid nodes. When M
=
M
x
=
M
y
=
M
z, C
n includes M
+
1
3 nodes, and M is called its expansion order. Figure 1 describes a box of order 2. Assume that a box C
n has center at c
n and radius r
n; Green’s function from source q
→ can be represented into(19)
G
r
→
,
q
→
=
∑
u
∈
C
n
π
u
,
C
n
q
→
G
r
→
,
u
→
,where the coefficients π
u
,
C
n
q
→ are to be determined and r
→ is an arbitrary point outside the box C
n, as shown in Figure 2.Figure 2
A 3D representation of the matching Green’s function values.Through selecting a group of sample pointsr
→
t
t
=
1
T on a spherical surface [12] with a bit larger radius R
n, the coefficients in (19) can be determined. When the number of sample points T is significantly greater than that of nodes in C
n, through matching Green’s function values at these sample points, the resultant systems of matrix equations can be solved by the least-square method. Without loss of generality, assume that integral domain V
m
F
,
± and V
m
H (V
n
F
,
± and V
n
H) are in the box C
m (C
n). Calculations of MoM matrix elements ultimately come down to calculating projections on these Gaussian points. Due to different inner and outer integrals, for anisotropic model, the 6 types of projection coefficients on the grid nodes are deduced as follows: (20)
T
→
V
n
F
,
±
V
n
H
1
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
κ
̿
p
→
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
κ
̿
p
→
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
→
u
,
C
n
T
V
n
F
,
±
V
n
H
1
,
T
→
V
n
F
,
±
V
n
H
2
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
f
→
i
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
→
u
,
C
n
T
V
n
F
,
±
V
n
H
2
,
T
V
n
F
,
±
V
n
H
3
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
∇
′
·
κ
̿
r
→
′
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
∇
·
κ
̿
p
→
i
f
→
n
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
∇
·
κ
-
-
p
→
i
f
→
n
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
V
n
F
,
±
V
n
H
3
,
T
V
n
F
,
±
V
n
H
4
=
∫
V
n
F
,
±
V
n
H
g
r
→
,
r
→
′
∇
′
·
f
→
n
r
→
′
d
v
′
=
∑
i
=
1
N
G
V
n
∑
u
=
1
M
C
n
w
i
∇
·
f
→
n
p
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
N
G
V
n
w
i
∇
·
f
→
n
p
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
V
n
F
,
±
V
n
H
4
,
T
∂
V
n
F
,
±
∂
V
n
H
5
=
∫
∂
V
n
F
,
±
∂
V
n
H
d
s
′
g
r
→
,
r
→
′
m
^
V
i
·
κ
̿
r
→
′
f
→
i
r
→
′
=
∑
s
=
1
4
∑
i
=
1
M
G
V
n
∑
u
=
1
M
C
n
w
i
′
m
^
V
i
·
κ
̿
s
→
i
f
→
n
s
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
s
=
1
4
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
M
G
V
n
s
w
i
′
m
^
V
i
·
κ
̿
s
→
i
f
→
n
s
→
i
π
u
,
C
n
p
→
i
=
∑
s
=
1
4
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
S
,
T
∂
V
n
F
,
±
∂
V
n
H
5
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
∂
V
n
F
,
±
∂
V
n
H
5
,
T
∂
V
n
H
6
=
∫
∂
V
n
H
d
s
′
g
r
→
,
r
→
′
m
^
V
i
·
f
→
i
r
→
′
=
∑
i
=
1
M
G
V
n
∑
u
=
1
M
C
n
w
i
′
m
^
V
i
·
f
→
n
s
→
i
π
u
,
C
n
p
→
i
g
r
→
,
r
→
u
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
∑
i
=
1
M
G
V
n
w
i
′
m
^
V
i
·
f
→
n
s
→
i
π
u
,
C
n
p
→
i
=
∑
u
=
1
M
C
n
g
r
→
,
r
→
u
π
u
,
C
n
T
∂
V
n
H
6
.In the above 6 formulas, there are three different integral regions: tetrahedral integral region V
m
(
n
)
F
,
±, V
m
(
n
)
H, and triangular integral region V
m
(
n
)
F
H
s, corresponding to their Gaussian points p
→
i
i
=
1
N
G
F
,
±, p
→
i
i
=
1
N
G
H, s
→
i
i
=
1
M
G and weights w
i
i
=
1
N
G
F
,
±, w
i
i
=
1
N
G
H, w
i
′
i
=
1
M
G, respectively. Besides difference between κ
-
-
(
r
→
) and κ
(
r
→
), note that T
(
5
)contains four triangular surface integrals for anisotropic model, while in [11] for isotropic model, T
(
5
)contains only one triangular surface integral. Obviously, (21)
T
→
V
n
F
V
n
H
1
=
T
→
V
n
F
,
+
1
+
T
→
V
n
F
,
-
1
T
→
V
n
H
1
,
T
→
V
n
F
V
n
H
2
=
T
→
V
n
F
,
+
2
+
T
→
V
n
F
,
-
2
T
→
V
n
H
2
,
T
V
n
F
V
n
H
3
=
T
V
n
F
,
+
3
+
T
V
n
F
,
-
3
T
V
n
H
3
,
T
V
n
F
V
n
H
4
=
T
V
n
F
,
+
4
+
T
V
n
F
,
-
4
T
V
n
H
4
,
T
V
n
F
V
n
H
5
=
T
V
n
F
,
+
5
+
T
V
n
F
,
-
5
T
V
n
H
5
,
T
∂
V
n
H
6
=
T
∂
V
n
H
6
Only about half SWG
.So, the ultimate π
→ or π matrix of the FFT uplink and downlink can be gotten:(22)
π
→
u
,
C
n
T
V
n
F
V
n
H
1
=
π
→
u
,
C
n
T
V
n
F
,
+
V
n
H
1
+
π
→
u
,
C
n
T
V
n
F
,
-
V
n
H
1
π
→
u
,
C
n
T
V
n
H
1
,
π
→
u
,
C
n
T
V
n
F
V
n
H
2
=
π
→
u
,
C
n
T
V
n
F
,
+
V
n
H
2
+
π
→
u
,
C
n
T
V
n
F
,
-
V
n
H
2
π
→
u
,
C
n
T
V
n
H
2
,
π
u
,
C
n
T
V
n
F
V
n
H
3
=
π
u
,
C
n
T
V
n
F
,
+
V
n
H
3
+
π
u
,
C
n
T
V
n
F
,
-
V
n
H
3
π
u
,
C
n
T
V
n
H
3
,
π
u
,
C
n
T
V
n
F
V
n
H
4
=
π
u
,
C
n
T
V
n
F
,
+
V
n
H
4
+
π
u
,
C
n
T
V
n
F
,
-
V
n
H
4
π
u
,
C
n
T
V
n
H
4
,
π
u
,
C
n
T
V
n
F
V
n
H
5
=
π
u
,
C
n
T
V
n
F
,
+
V
n
H
5
+
π
u
,
C
n
T
V
n
F
,
-
V
n
H
5
π
u
,
C
n
T
V
n
H
5
,
π
u
,
C
n
T
∂
V
n
H
6
=
π
u
,
C
n
T
∂
V
n
H
6
Only about half SWG
.Subscript u
=
1,2
,
…
,
M
C
n. At last, the far-field matrix elements can be expressed in such following form:(23)
a
m
n
F
F
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
F
2
G
u
,
v
π
→
v
,
C
n
T
V
n
F
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
F
3
,
a
m
n
F
H
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
F
2
G
u
,
v
π
→
v
,
C
n
T
V
n
H
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
H
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
H
3
,
a
m
n
H
F
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
H
2
G
u
,
v
π
→
v
,
C
n
T
V
n
F
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
F
3
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
F
5
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
F
3
,
a
m
n
H
H
=
-
k
0
2
∑
u
∈
C
m
∑
v
∈
C
n
π
→
u
,
C
m
T
V
m
H
2
G
u
,
v
π
→
v
,
C
n
T
V
n
H
1
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
H
5
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
H
3
+
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
H
5
-
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
H
3
.Compared with ordinary media, a
m
n
F
F adds one term -
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
F
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5, while a
m
n
H
F adds two terms -
∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
4
G
u
,
v
π
v
,
C
n
T
V
n
F
5 and ∑
u
∈
C
m
∑
v
∈
C
n
π
u
,
C
m
T
V
m
H
6
G
u
,
v
π
v
,
C
n
T
V
n
F
5. So, it is more complicated relatively.


## 2.4. Near-Field Preconditioning Technique

There are a few preconditioning techniques successfully developed for some fast algorithms, solving MoM linear systems arising from surface integral equation (SIE). However, the preconditioning techniques applied to FFT-based VIE are less studied. This section investigates convergence effect of near-field preconditioning technique for the anisotropic VIE-FG-FFT. In the FFT-based methods, only near-field matrix elements are stored in the memory, and because of the characteristics of 3D Green’s function, near-field elements are generally larger than far-field matrix elements, so near-field matrix can be seen as sparsification of the global matrix by abandoning small matrix elements. The near-field matrix could be used to construct DILU preconditioner [13]. However, different from the MLFMA, the near-field matrix elements in the FG-FFT are not real; namely, the near-field matrix is built with FFT revision. Experiments demonstrate it is not right to directly use the revised near-field matrix to construct preconditioner. Now we use real near-field matrix A
n
e
a
r, which can be separated into three parts:(24)
A
near
=
L
near
+
D
near
+
U
near
,where L
n
e
a
r is lower triangular part of A
n
e
a
r; U
n
e
a
r is upper triangular part of A
n
e
a
r; D
n
e
a
r is main diagonal part of A
n
e
a
r. Hence, the new preconditioning matrix is(25)
M
=
D
near
+
L
near
D
near
-
1
D
near
+
U
near
.Namely, (26)
M
=
D
near
+
L
near
I
+
D
near
-
1
U
near
.This preconditioner is easy to implement and additional time cost of each CG iteration is extremely small. Numerical examples have proved the good effectiveness in Section 3.


## 3. Numerical Results

In this section, several anisotropic examples are provided to demonstrate the validity, accuracy, and efficiency of the VIE-FG-FFT for anisotropic media. In all the examples, the expansion orderM is always chosen as 2 (generally 2 is enough) and the Cartesian grid spacing sizes in three directions are always selected to be the same. Besides, it is assumed that the incident wave excitation is plane wave which has the propagation direction along the z+-axis and the polarization direction in the x-axis. The parameters of computation platform are Xeon E5450 CPU and 12 G memory.Example 1.
The scattering from an anisotropic spherical shell is calculated by the VIE-FG-FFT. The inner and outer radius of the spherical shell are 0.9λ and 1.0λ, respectively, shown in Figure 3.
The permittivity tensor of the spherical shell has nonzero elements:ε
r
,
x
x
=
ε
r
,
y
y
=
2.2
-
j
0.2, ε
r
,
z
z
=
1.5
-
j
0.1. The object is modeled by 17360 tetrahedrons with the average edge length 0.1λ, and the number of unknowns N is 38231. The Cartesian grid spacing size is set to be 0.1λ and 0.2λ. The bistatic RCS curves of the VIE-FG-FFT are compared with the analytical solution, shown in Figure 4. It can be seen that the three curves coincide very well, which demonstrates the validity of the new method and its stability (insensitivity to grid spacing sizes).Figure 3
An anisotropic spherical shell.Figure 4
The bistatic RCS curves of an anisotropic spherical shell.Example 2.
An anisotropic hollow lossy dielectric box is considered in this example. The three-dimensional sizes of the box are shown in Figure5. The relative permittivity tensor has nonzero elements: (27)
ε
r
,
x
x
=
ε
r
,
y
y
=
2.5
-
j
0.2
,
ε
r
,
z
z
=
1.8
-
j
0.1
,
ε
r
,
x
y
=
2.6
-
j
0.2
,
ε
r
,
z
x
=
2.4
-
j
0.3
.The box is modeled by 19301 tetrahedrons with the average edge length 0.1λ, and the number of unknowns N is 42826. The Cartesian grid spacing size is set to be 0.1λ. The two bistatic RCS curves obtained by the FG-FFT with or without DILU preconditioner are shown in Figure 6. They are all compared with that of conventional VIE with domain decomposition technique [14]. It can be seen that the three curves coincide very well, which demonstrates the validity of the proposed method. It can be seen from Table 1 that the DILU preconditioner can decrease the total solving time and speed up the convergence of CG solution. This example cannot be directly calculated by conventional VIE-MoM on this computational platform, because of memory limited, shown in Table 2.Table 1
Comparison of the CG solving time.
Method
One time of iteration (s)
Number of iterations
Total time of iteration (s)
FG-FFT
5.69
108
614.5
FG-FFT + DILU
5.95
85
505.8Table 2
Comparison of the memory requirement.
Method
VIE-FG-FFT (Pre)
VIE-MoM
Memory
913 MB
Memory limit exceededFigure 5
The three-dimensional sizes of a hollow dielectric box.Figure 6
The bi-station RCS curves whenh
=
0.1
λ.Example 3.
Another anisotropic example is provided. An inhomogeneous anisotropic solid dielectric box is considered here. The three-dimensional sizes of the box are shown in Figure7. The relative permittivity tensors have nonzero elements:(28)
ε
r
1
,
x
x
=
ε
r
1
,
y
y
=
2.2
-
j
0.2
,
ε
r
1
,
z
z
=
1.5
-
j
0.1
,
ε
r
2
,
x
x
=
ε
r
2
,
y
y
=
2.4
-
j
0.3
,
ε
r
2
,
z
z
=
1.6
-
j
0.2
,
ε
r
3
,
x
x
=
ε
r
3
,
y
y
=
2.6
-
j
0.2
,
ε
r
3
,
z
z
=
1.7
-
j
0.1
.The box is modeled by 17935 tetrahedrons with the average edge length 0.1λ, and the number of unknowns is 37470. The Cartesian grid spacing size is set to be 0.1λ. The two bistatic RCS curves obtained by the FG-FFT with or without DILU preconditioner and the RCS curve (blue) of FEM-based commercial software are shown in Figure 8. It can also be seen that the three curves coincide very well, which demonstrates the validity of the proposed method. It can also be seen from Table 3 that the DILU preconditioner has good effectiveness on speeding up convergence of CG solution. This example also exceeds memory limit of our computational platform if using conventional VIE-MoM, shown in Table 4.Table 3
Comparison of the CG solving time.
Method
One time of iteration (s)
Number of iterations
Total time of iteration (s)
FG-FFT
3.63
69
250.5
FG-FFT + DILU
4.07
48
195.4Table 4
Comparison of the memory requirement.
Method
VIE-FG-FFT (Pre)
VIE-MoM
Memory
668 MB
Memory limit exceededFigure 7
The three-dimensional sizes of a solid dielectric box.Figure 8
The bi-station RCS curves whenh
=
0.1
λ.


## 4. Conclusion

The Fast Fourier Transform of fitting Green’s function technique has high accuracy and requires less pretreatment time for dealing with the volume integral equations, so in this paper, it is selected to accelerate anisotropic VIE model. Actually, this work perfects the application system of the Fast Fourier Transform of fitting Green’s function technique from surface integral equation (SIE) to volume integral equation of isotropic media (isotropic VIE model), again to volume integral equation of anisotropic media (anisotropic VIE model). In this work, the proposed method has been compared with analytic result, anisotropic VIE with domain decomposition technique, and commercial software; the numerical examples demonstrate its validity and obvious memory advantages for solving inhomogeneous, electrically large, anisotropic, dielectric objects. Besides, preconditioning technique has been investigated as a complementary work. The DILU preconditioner can speed up the convergence of CG solution and decrease the total solution time. In the future, this method can further investigate other preconditioners and how to be applied to high order volume basis functions.



---

*Source: 123739-2015-12-31.xml*