# A Method for Fast Generating the VIE-MoM Matrix with the SWG Functions

Authors: Shu-Wen Chen; Hou-Xing Zhou; Wei-Dong Li; Guang Hua; Wei Hong
Journal: International Journal of Antennas and Propagation
 (2013)

---

## Abstract

A method for fast generating the VIE-MoM matrix with the SWG functions is proposed. This method makes the integrals representing the shared data independent of any quadratures
and hence can ensure both the efficiency and accuracy. In the VIE-MoM, the time cost of filling the VIE-MoM matrix is prominent especially
for an electrically large object. With the proposed method, the time cost of filling the VIE-MoM matrix can
be reduced by about 85% compared with the traditional basis-basis matching method. Numerical examples
are provided to demonstrate the efficiency and accuracy of the proposed method.


---

## Body

## 1. Introduction

Study of electromagnetic scattering characteristic of 3D inhomogeneous dielectric objects has significance in many important fields [1]. The method of moments (MoM) [2] is frequently used to solve electromagnetic scattering problems. No matter based on the surface integral equation (SIE) [3] or on the volume integral equation (VIE) [4], the MoM matrixes are always dense. Therefore, filling the MoM matrix will consume a large amount of CPU time especially for an electrically large object.For the VIE-MoM with the Schauert-Wilton-Glisson (SWG) functions, Schaubert et al. did prophecy in [4]. Considerable savings in computation time can be achieved by expressing the VIE-MoM matrix elements into four independent scalar integrals over each tetrahedron. But so far, there is no publication for the concrete realization of this vision. For the SIE-MoM with the Rao-Wilton-Glisson (RWG) functions [3], Zhang proposed a triangle-triangle scheme in [5] to accelerate filling the SIE-MoM matrix. The interactions between a pair of triangle patches associate with up to 9 matrix elements, which is the foundation of the triangle-triangle scheme. Zhang's method has many applications [6–9]. Actually, the ideas of both Schauert and Zhang are similar and rely on using the midpoint formula for the testing integrals. However, the midpoint formula cannot reach the required accuracy in many practical problems.In this paper, a method for fast generating the VIE-MoM matrix with the SWG functions is proposed. This method makes the integrals representing the shared data independent of any quadratures and hence can ensure both the efficiency and accuracy. With the proposed method, the time cost of filling the VIE-MoM matrix can be reduced by about 85% compared with the traditional basis-basis matching method. Numerical examples are provided to demonstrate the efficiency and accuracy of the proposed method.


## 2. Integral Formulae of the VIE-MoM Matrix Elements

### 2.1. The Volume Integral Equation

LetV denote the volume domain of an inhomogeneous dielectric object with a relative permittivity εr(r⃑) and a relative permeability μr(r⃑)=1. The permittivity and the permeability of the free space are ε0(r⃑) and μ0(r⃑), respectively. Let E⃑i be the incident field and E⃑s the scattered field; then the total field E⃑tot can be expressed as the sum of E⃑i and E⃑s:
(1)E⃑tot(r⃑)=E⃑i(r⃑)+E⃑s(r⃑).The scattered field due to the volume polarization currentJ⃑V(r⃑) in V can be written as
(2)E⃑s(r⃑)=-jk0η0(I̿+1k02∇∇·)∫Vg(r⃑,r⃑′)J⃑V(r⃑′)dv′,
where
(3)J⃑V(r⃑)=jk0η0(εr(r⃑)-1)E⃑tot=jωκ(r⃑)D⃑tot(r⃑),κ(r⃑)=(εr(r⃑)-1)εr(r⃑),k0=ωμ0ε0,η0=μ0ε0.


### 2.2. Building the MoM Model

The electric flux densityD⃑tot(r⃑) can be selected as the unknown function because it is continuous along the normal of the medium interface of V.AfterV being discretized by tetrahedra, D⃑tot(r⃑)/ε0 can be expanded with the SWG functions {f⃑n(r⃑)} as follows:
(4)D⃑tot(r⃑)ε0=∑n=1N‍xnf⃑n(r⃑),
where f⃑n(r⃑) is defined on a pair of tetrahedra as
(5)f⃑n(r⃑)={f⃑n+(r⃑)=|Sn|3|Vn+|(r⃑-r⃑n,free+),r⃑∈Vn+,f⃑n-(r⃑)=|Sn|3|Vn-|(r⃑-r⃑n,free-),r⃑∈Vn-,
whose geometry is shown in Figure 1. The details about the SWG functions are referred to in [4]. It should be noted that, for modeling D⃑tot(r⃑)/ε0 at the outer boundary of  V, the so-called half-SWG functions are also needed.Figure 1
The geometry of the SWG function.When (4) is substituted into (1) and after the testing procedure, the corresponding MoM equation can be built below:
(6)Ax=b,
where  A=(aij)N×N is a N×N matrix,x=[x1,x2,…,xN]T is the unknown vector,b=[b1,b2,…,bN]T represents the incident vector.More precisely,(7)aij=∫Vif⃑i(r⃑)·[∫Vjg(r⃑,r⃑′)κ(r⃑′)fj(r⃑′)dv′εr-1(r⃑)f⃑j(r⃑)-(k02+∇∇·)===========×∫Vjg(r⃑,r⃑′)κ(r⃑′)===========×fj(r⃑′)dv′∫Vjg(r⃑,r⃑′)κ(r⃑′)fj(r⃑′)dv′]dv,(8)bn=∫Vnf⃑n(r⃑)·E⃑i(r⃑)dv.Filling the matrixA by using (7) (the basis-basis scheme) includes a lot of redundant calculations. This disadvantage can be overcome by the new method presented in the following.


### 2.3. Fast Filling the VIE-MoM Matrix

Rewrite (7) as
(9)aij=Iij(1)-Iij(2)-Iij(3)+Iij(4)+Iij(5)-Iij(6),
where(10a)Iij(1)=∫Viεr-1(r⃑)f⃑i(r⃑)·f⃑j(r⃑)dv,(10b)Iij(2)=k02∫Vif⃑i(r⃑)dv·∫Vjg(r⃑,r⃑′)κ(r⃑′)f⃑j(r⃑′)dv′,(10c)Iij(3)=∫Vi[∇·f⃑i(r⃑)]dv×∫∂Vjg(r⃑,r⃑′)κ(r⃑′)[m^Vj·f⃑j(r⃑′)]ds′,(10d)Iij(4)=∫Vi[∇·f⃑i(r⃑)]dv×∫Vjg(r⃑,r⃑′)[∇′·κ(r⃑′)f⃑j(r⃑′)]dv′,(10e)Iij(5)=∫∂Vids[m^Vi·f⃑i(r⃑)]×∫∂Vjg(r⃑,r⃑′)κ(r⃑′)[m^Vj·f⃑j(r⃑′)]ds′,(10f)Iij(6)=∫∂Vids[m^Vi·f⃑i(r⃑)]×∫Vjg(r⃑,r⃑′)[∇′·κ(r⃑′)f⃑j(r⃑′)]dv′,where ∂Vi is the outer surface of  Vi and  m^Vi represents the unit outward normal vector of ∂Vi.Usually,Vi consists of two tetrahedron subdomains; that is,
(11)Vi=Vi+∪Vi-.
In such a case, Vi is called a complete support set, and ∂Vi consists of 6 triangle subsurfaces. It is readily known that if Vi is complete, all the integrals in (10a)–(10f) having ∂Vi as their integration regions are equal to zero. When Vi=Vi+ or Vi=Vi-, Vi is called an incomplete support set or a half support set. Clearly, when Vi is incomplete, ∂Vi consists of 4 triangle subsurfaces. In such a case, every integral in (10a)–(10f) having ∂Vi as its integration region can be reduced as one on the triangle sub-surface Si of ∂Vi (see Figure 1).It can be found that every integral in (10a)–(10d) and (10f) for a fixed element aij includes a lot of calculations shared by other elements. Removing these redundant calculations can greatly improve the filling efficiency of the matrix. Here, we consider only (10b), and the others will be presented in the appendix.Without loss of generality, we assume that bothVi and Vj are complete; then
(12)Iij(2)=Ii+,j+(2)+Ii+,j-(2)+Ii-,j+(2)+Ii-,j-(2),
where Ii±,j±(2) has Vi± and Vj± as its outer and inner integration regions, respectively. The calculation formulae of the four terms in the right hand of (12) are very analogous, so only one needs to be provided. For doing this, we now consider
(13)Ii+,j+(2)=k02|Si|3|Vi+||Sj|3|Vj+|×∫Vi+∫Vj+(r⃑-r⃑i,free+)·(r⃑′-r⃑j,free+)×κ(r⃑′)g(r⃑,r⃑′)dvdv′.
After expanding (r⃑-r⃑i,free+)·(r⃑′-r⃑j,free+), we obtain
(14)Ii+,j+(2)=k02|Si||Sj|9×(Q⃑3(2)Q1(2)-r⃑j,free+·Q⃑2(2)-r⃑i,free+·Q⃑3(2)+r⃑i,free+·r⃑j,free+Q4(2)Q⃑2(2)r⃑i,free+),
where(15a)Q1(2)=1|Vi+||Vj+|∫Vi+∫Vj+r⃑·r⃑′κ(r⃑′)g(r⃑,r⃑′)dvdv′,(15b)Q2(2)=1|Vi+||Vj+|∫Vi+∫Vj+r⃑κ(r⃑′)g(r⃑,r⃑′)dvdv′,(15c)Q3(2)=1|Vi+||Vj+|∫Vi+∫Vj+r⃑′κ(r⃑′)g(r⃑,r⃑′)dvdv′,(15d)Q4(2)=1|Vi+||Vj+|∫Vi+∫Vj+κ(r⃑′)g(r⃑,r⃑′)dvdv′.Note that the formulae (15a)–(15d) are only related with the global numbers of the tetrahedra Vi+ and Vj+. Therefore, these numerical results can be shared by multiple matrix elements, which is the reason why the redundant calculations can be removed. Generally, the interactions between a pair of tetrahedra associate with up to 16 matrix elements.Similar analysis can also be done for the remaining integrals in (10a)–(10f), and the corresponding calculation formulae are listed in the appendix. In most cases, in the MoM model discussed, triangle subdomains are far less than tetrahedron subdomains, and hence calculation amount corresponding to the formulae (10c), (10e), and (10f) is much less than that of the formulae (10a), (10b), and (10d). Therefore, the time cost related with Iij(1), Iij(2), and Iij(4) is particularly prominent.It should be pointed out that the formulae (15a)–(15d), (A.2a)–(A.2c), (A.4), (A.6a), (A.6b), (A.9a), and (A.9b) provide the shared data for the corresponding matrix elements and are independent of any quadratures. Therefore, the proposed method can ensure both the efficiency and accuracy.


## 2.1. The Volume Integral Equation

LetV denote the volume domain of an inhomogeneous dielectric object with a relative permittivity εr(r⃑) and a relative permeability μr(r⃑)=1. The permittivity and the permeability of the free space are ε0(r⃑) and μ0(r⃑), respectively. Let E⃑i be the incident field and E⃑s the scattered field; then the total field E⃑tot can be expressed as the sum of E⃑i and E⃑s:
(1)E⃑tot(r⃑)=E⃑i(r⃑)+E⃑s(r⃑).The scattered field due to the volume polarization currentJ⃑V(r⃑) in V can be written as
(2)E⃑s(r⃑)=-jk0η0(I̿+1k02∇∇·)∫Vg(r⃑,r⃑′)J⃑V(r⃑′)dv′,
where
(3)J⃑V(r⃑)=jk0η0(εr(r⃑)-1)E⃑tot=jωκ(r⃑)D⃑tot(r⃑),κ(r⃑)=(εr(r⃑)-1)εr(r⃑),k0=ωμ0ε0,η0=μ0ε0.


## 2.2. Building the MoM Model

The electric flux densityD⃑tot(r⃑) can be selected as the unknown function because it is continuous along the normal of the medium interface of V.AfterV being discretized by tetrahedra, D⃑tot(r⃑)/ε0 can be expanded with the SWG functions {f⃑n(r⃑)} as follows:
(4)D⃑tot(r⃑)ε0=∑n=1N‍xnf⃑n(r⃑),
where f⃑n(r⃑) is defined on a pair of tetrahedra as
(5)f⃑n(r⃑)={f⃑n+(r⃑)=|Sn|3|Vn+|(r⃑-r⃑n,free+),r⃑∈Vn+,f⃑n-(r⃑)=|Sn|3|Vn-|(r⃑-r⃑n,free-),r⃑∈Vn-,
whose geometry is shown in Figure 1. The details about the SWG functions are referred to in [4]. It should be noted that, for modeling D⃑tot(r⃑)/ε0 at the outer boundary of  V, the so-called half-SWG functions are also needed.Figure 1
The geometry of the SWG function.When (4) is substituted into (1) and after the testing procedure, the corresponding MoM equation can be built below:
(6)Ax=b,
where  A=(aij)N×N is a N×N matrix,x=[x1,x2,…,xN]T is the unknown vector,b=[b1,b2,…,bN]T represents the incident vector.More precisely,(7)aij=∫Vif⃑i(r⃑)·[∫Vjg(r⃑,r⃑′)κ(r⃑′)fj(r⃑′)dv′εr-1(r⃑)f⃑j(r⃑)-(k02+∇∇·)===========×∫Vjg(r⃑,r⃑′)κ(r⃑′)===========×fj(r⃑′)dv′∫Vjg(r⃑,r⃑′)κ(r⃑′)fj(r⃑′)dv′]dv,(8)bn=∫Vnf⃑n(r⃑)·E⃑i(r⃑)dv.Filling the matrixA by using (7) (the basis-basis scheme) includes a lot of redundant calculations. This disadvantage can be overcome by the new method presented in the following.


## 2.3. Fast Filling the VIE-MoM Matrix

Rewrite (7) as
(9)aij=Iij(1)-Iij(2)-Iij(3)+Iij(4)+Iij(5)-Iij(6),
where(10a)Iij(1)=∫Viεr-1(r⃑)f⃑i(r⃑)·f⃑j(r⃑)dv,(10b)Iij(2)=k02∫Vif⃑i(r⃑)dv·∫Vjg(r⃑,r⃑′)κ(r⃑′)f⃑j(r⃑′)dv′,(10c)Iij(3)=∫Vi[∇·f⃑i(r⃑)]dv×∫∂Vjg(r⃑,r⃑′)κ(r⃑′)[m^Vj·f⃑j(r⃑′)]ds′,(10d)Iij(4)=∫Vi[∇·f⃑i(r⃑)]dv×∫Vjg(r⃑,r⃑′)[∇′·κ(r⃑′)f⃑j(r⃑′)]dv′,(10e)Iij(5)=∫∂Vids[m^Vi·f⃑i(r⃑)]×∫∂Vjg(r⃑,r⃑′)κ(r⃑′)[m^Vj·f⃑j(r⃑′)]ds′,(10f)Iij(6)=∫∂Vids[m^Vi·f⃑i(r⃑)]×∫Vjg(r⃑,r⃑′)[∇′·κ(r⃑′)f⃑j(r⃑′)]dv′,where ∂Vi is the outer surface of  Vi and  m^Vi represents the unit outward normal vector of ∂Vi.Usually,Vi consists of two tetrahedron subdomains; that is,
(11)Vi=Vi+∪Vi-.
In such a case, Vi is called a complete support set, and ∂Vi consists of 6 triangle subsurfaces. It is readily known that if Vi is complete, all the integrals in (10a)–(10f) having ∂Vi as their integration regions are equal to zero. When Vi=Vi+ or Vi=Vi-, Vi is called an incomplete support set or a half support set. Clearly, when Vi is incomplete, ∂Vi consists of 4 triangle subsurfaces. In such a case, every integral in (10a)–(10f) having ∂Vi as its integration region can be reduced as one on the triangle sub-surface Si of ∂Vi (see Figure 1).It can be found that every integral in (10a)–(10d) and (10f) for a fixed element aij includes a lot of calculations shared by other elements. Removing these redundant calculations can greatly improve the filling efficiency of the matrix. Here, we consider only (10b), and the others will be presented in the appendix.Without loss of generality, we assume that bothVi and Vj are complete; then
(12)Iij(2)=Ii+,j+(2)+Ii+,j-(2)+Ii-,j+(2)+Ii-,j-(2),
where Ii±,j±(2) has Vi± and Vj± as its outer and inner integration regions, respectively. The calculation formulae of the four terms in the right hand of (12) are very analogous, so only one needs to be provided. For doing this, we now consider
(13)Ii+,j+(2)=k02|Si|3|Vi+||Sj|3|Vj+|×∫Vi+∫Vj+(r⃑-r⃑i,free+)·(r⃑′-r⃑j,free+)×κ(r⃑′)g(r⃑,r⃑′)dvdv′.
After expanding (r⃑-r⃑i,free+)·(r⃑′-r⃑j,free+), we obtain
(14)Ii+,j+(2)=k02|Si||Sj|9×(Q⃑3(2)Q1(2)-r⃑j,free+·Q⃑2(2)-r⃑i,free+·Q⃑3(2)+r⃑i,free+·r⃑j,free+Q4(2)Q⃑2(2)r⃑i,free+),
where(15a)Q1(2)=1|Vi+||Vj+|∫Vi+∫Vj+r⃑·r⃑′κ(r⃑′)g(r⃑,r⃑′)dvdv′,(15b)Q2(2)=1|Vi+||Vj+|∫Vi+∫Vj+r⃑κ(r⃑′)g(r⃑,r⃑′)dvdv′,(15c)Q3(2)=1|Vi+||Vj+|∫Vi+∫Vj+r⃑′κ(r⃑′)g(r⃑,r⃑′)dvdv′,(15d)Q4(2)=1|Vi+||Vj+|∫Vi+∫Vj+κ(r⃑′)g(r⃑,r⃑′)dvdv′.Note that the formulae (15a)–(15d) are only related with the global numbers of the tetrahedra Vi+ and Vj+. Therefore, these numerical results can be shared by multiple matrix elements, which is the reason why the redundant calculations can be removed. Generally, the interactions between a pair of tetrahedra associate with up to 16 matrix elements.Similar analysis can also be done for the remaining integrals in (10a)–(10f), and the corresponding calculation formulae are listed in the appendix. In most cases, in the MoM model discussed, triangle subdomains are far less than tetrahedron subdomains, and hence calculation amount corresponding to the formulae (10c), (10e), and (10f) is much less than that of the formulae (10a), (10b), and (10d). Therefore, the time cost related with Iij(1), Iij(2), and Iij(4) is particularly prominent.It should be pointed out that the formulae (15a)–(15d), (A.2a)–(A.2c), (A.4), (A.6a), (A.6b), (A.9a), and (A.9b) provide the shared data for the corresponding matrix elements and are independent of any quadratures. Therefore, the proposed method can ensure both the efficiency and accuracy.


## 3. Testing the Efficiency

In this section, three examples are provided to test the efficiency of the presented method. The symmetric quadrature with 4 integral points is applied to each integral on tetrahedra, and the symmetric quadrature with 3 integral points is applied to each integral on triangles. The computations include the serial computation and the parallel computation based on OpenMP [10]. It is assumed that the excitation is plane wave that has the propagation direction along the positive z-axis and the polarization direction in the x-axis. Also, λ always denotes the wavelength in the free space.Example A (homogeneous ball).
A homogeneous ball with the relative permittivity 2.0 is considered. The radius of the ball is0.4λ. The volume domain is discretized by 4291 tetrahedra, and the number of unknowns is 8925.
The bistatic RCS curves for the fixed angleϕ=0 are shown in Figure 2. It can be seen from Figure 2 that the results from the MoM using the proposed method, the MoM using the basis-basis scheme, and the Mie series [11] are consistent, which implies that our computer program is correct. Time cost comparison is recorded in Table 1 (the part with A).Table 1
The data from the three examples.
Examples
Method
1 thread
8 threads
Time cost
Time ratio
Reduced by
Time cost
Time ratio
Reduced by
A
Basis-basis scheme
2038.4 s
8.748
88.6%
261.3 s
8.539
88.3%
Proposed method
233.0 s
30.6 s
B
Basis-basis scheme
718.3 s
7.537
86.7%
91.7 s
7.455
86.5%
Proposed method
95.3 s
12.3 s
C
Basis-basis scheme
324.8 s
6.867
85.4%
41.2 s
6.754
85.1%
Proposed method
47.3 s
6.1 sFigure 2
The bistatic RCS curves of the homogeneous ball.Example B (partitioned homogeneous cuboid).
A partitioned homogeneous cuboid with 3 different materials is considered in Figure3(a). The volume domain is discretized by 2511 tetrahedra, and the number of unknowns is 5312. Time cost comparison is recorded in Table 1 (the part with B).The partitioned homogeneous cuboid and the inhomogeneous cube.
(a)
The cuboid
(b)
The cubeExample C (inhomogeneous cube).
An  inhomogeneous cube is considered in Figure3(b). The relative permittivity of the cube varies according to
(16)εr(r⃑)=(x+y+z)9+1.5.
The volume domain is discretized by 1635 tetrahedra, and the number of unknowns is 3486. Time cost comparison is recorded in Table 1 (the part with C).
It can be seen from the above examples that the efficiency of the proposed method is very satisfactory whether in serial computation case or in parallel computation case. With the proposed method in this paper, the time cost of filling VIE-MoM matrix is reduced by about 85% compared with the traditional basis-basis scheme.


## 4. Testing the Accuracy

The proposed method makes the integrals representing the shared data independent of any quadratures and hence can be easily combined with the proper quadratures to achieve better numerical accuracy. We will verify this fact by using quadratures of different algebraic degrees (or orders).Symmetric quadrature rules on triangles and tetrahedra will be adopted [12]. For brevity, we will only mention the algebraic degrees on tetrahedra. The MoM matrix calculated by using the quadrature of algebraic degree 14 (contains 236 points) on the tetrahedron is regarded as the benchmark. The relative error (RE) of the MoM matrix is defined by
(17)RE=∥ZBenchmark-ZAM∥F∥ZBenchmark∥F,
where ∥·∥F denotes the Frobenius norm of the matrix [13] and the superscript “AM” means “approximation matrix” calculated by quadratures of lower algebraic degrees.Here, the homogeneous ball with the relative permittivity 2.0 in Example A is considered once again. The REs obtained by using quadratures of different algebraic degrees are recorded in Table2. It can be clearly seen from Table 2 that the accuracy of the MoM matrix can be improved with increasing the algebraic degree of the quadrature.Table 2
The relative error changing with the algebraic degree.
Algebraic degree
1
2
3
4
5
6
Number of points
1 (midpoint formula)
4
8
14
14
24
Relative error
1.61721
e
-
1
2.32651
e
-
2
2.11931
e
-
2
1.82897
e
-
2
9.42591
e
-
3
4.78024
e
-
3It should be pointed out that the amount of memory required by the MoM matrix (1215.45 MB in this example) itself is not changed due to the application of fast filling scheme. The interactions between one pair of tetrahedron subdomains can contribute to up to 16 matrix elements, and the calculations of shared data need to be temporarily stored, such as the values of Green's function, which requires proper auxiliary memory. In each thread, each time, the interactions between only one pair of tetrahedron subdomains need small storage in which the part required by the values of Green's function is dominant, and this small storage can be evaluated asO(16n2) bytes when using double precision, where n is the number of integral points on a tetrahedron subdomain. Clearly, this small storage can be ignored compared with the amount of storage required by the MoM matrix. Besides, when the algebraic degree of the quadrature is increased, the auxiliary storage for helping the fast calculation will be increased, but the acceleration effect will be more apparent due to shared data increase.


## 5. Conclusions

In this paper, a fast method for generating the VIE-MoM matrix using the SWG functions has been presented in details. By means of the inherent structure property of the SWG functions, those integrals having tetrahedra and/or triangles as their outer or inner integration regions can be rewritten as a sum of several terms that include integrals which can be shared by many matrix elements, significantly reducing time cost of filling the VIE-MoM matrix.The prominent feature of the proposed method is that the integrals representing the shared data are independent of any quadratures, which is the reason for that the proposed method can ensure both the efficiency and accuracy.It must be emphasized that the foundation of the proposed method is the structure property of the SWG functions. If other forms of basis functions such as the higher order basis functions are needed in some MoM models, additional approaches need to be explored.



---

*Source: 480790-2013-09-04.xml*