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

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.

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

Let 𝑉 denote the volume domain of an inhomogeneous dielectric object with a relative permittivity 𝜀_𝑟(^⇀𝑟) and a relative permeability 𝜇_𝑟(^⇀𝑟)=1. The permittivity and the permeability of the free space are 𝜀₀(^⇀𝑟) and 𝜇₀(^⇀𝑟), respectively. Let ^⇀𝐸^𝑖 be the incident field and ^⇀𝐸^𝑠 the scattered field; then the total field ^⇀𝐸^tot can be expressed as the sum of ^⇀𝐸^𝑖 and ^⇀𝐸^𝑠: ^⇀𝐸^tot⇀𝑟=^⇀𝐸^𝑖⇀𝑟+^⇀𝐸^𝑠⇀𝑟. (1)

The scattered field due to the volume polarization current ^⇀𝐽_𝑉(^⇀𝑟) in 𝑉 can be written as ^⇀𝐸^𝑠⇀𝑟=−𝑗𝑘₀𝜂₀1𝐼+𝑘²₀∇∇⋅_𝑉𝑔^⇀𝑟,^⇀𝑟^′⇀𝐽_𝑉^⇀𝑟^′𝑑𝑣^′, (2) where ^⇀𝐽_𝑉^⇀𝑟𝑘=𝑗⁰𝜂₀𝜀_𝑟^⇀𝑟−1^⇀𝐸^tot=𝑗𝜔𝜅^⇀𝑟𝐷^⇀tot⇀𝑟,𝜅^⇀𝑟=𝜀_𝑟^⇀𝑟−1𝜀_𝑟^⇀𝑟,𝑘₀=𝜔√𝜇₀𝜀₀,𝜂₀=𝜇₀𝜀₀. (3)

2.2. Building the MoM Model

The electric flux density 𝐷^⇀tot(^⇀𝑟) can be selected as the unknown function because it is continuous along the normal of the medium interface of 𝑉.

After 𝑉 being discretized by tetrahedra, 𝐷^⇀tot(^⇀𝑟)/𝜀₀ can be expanded with the SWG functions {^⇀𝑓_𝑛(^⇀𝑟)} as follows: 𝐷^⇀tot⇀𝑟𝜀₀=^𝑁_𝑛=1𝑥_𝑛^⇀𝑓_𝑛^⇀𝑟, (4) where ^⇀𝑓_𝑛(^⇀𝑟) is defined on a pair of tetrahedra as ^⇀𝑓_𝑛^⇀𝑟=^⇀𝑓⁺_𝑛^⇀𝑟=𝑆_𝑛3𝑉⁺_𝑛^⇀𝑟−^⇀𝑟⁺_𝑛,free,^⇀𝑟∈𝑉⁺_𝑛,_⇀𝑓⁻_𝑛^⇀𝑟=𝑆_𝑛3𝑉⁻_𝑛^⇀𝑟−^⇀𝑟⁻_𝑛,free,^⇀𝑟∈𝑉⁻_𝑛, (5) 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 𝐷^⇀tot(^⇀𝑟)/𝜀₀ at the outer boundary of  𝑉, 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) where  𝐴=(𝑎_𝑖𝑗)_𝑁×𝑁 is a 𝑁×𝑁 matrix,  𝑥=[𝑥₁,𝑥₂,…,𝑥_𝑁]^𝑇 is the unknown vector,  𝑏=[𝑏₁,𝑏₂,…,𝑏_𝑁]^𝑇 represents the incident vector.

More precisely, 𝑎_𝑖𝑗=_𝑉𝑖^⇀𝑓_𝑖^⇀𝑟⋅𝜀_𝑟^−1⇀𝑟^⇀𝑓_𝑗^⇀𝑟−𝑘²₀×+∇∇⋅_𝑉𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟^′×𝑓_𝑗^⇀𝑟^′𝑑𝑣^′𝑏𝑑𝑣,_𝑛=_𝑉𝑛^⇀𝑓_𝑛^⇀𝑟⋅^⇀𝐸^𝑖⇀𝑟𝑑𝑣. (7)(8)

Filling the matrix 𝐴 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) where𝐼⁽¹⁾_𝑖𝑗=_𝑉𝑖𝜀_𝑟^−1⇀𝑟^⇀𝑓_𝑖^⇀𝑟⋅^⇀𝑓_𝑗^⇀𝑟𝐼𝑑𝑣,_(2)𝑖𝑗=𝑘²_0𝑉𝑖^⇀𝑓_𝑖^⇀𝑟𝑑𝑣⋅_𝑉𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟^′⇀𝑓_𝑗^⇀𝑟^′𝑑𝑣^′,𝐼⁽³⁾_𝑖𝑗=_𝑉𝑖∇⋅^⇀𝑓_𝑖^⇀𝑟×𝑑𝑣_𝜕𝑉𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟^′𝑚_𝑉𝑗⋅^⇀𝑓_𝑗^⇀𝑟^′𝑑𝑠^′,𝐼⁽⁴⁾_𝑖𝑗=_𝑉𝑖∇⋅^⇀𝑓_𝑖^⇀𝑟×𝑑𝑣_𝑉𝑗𝑔^⇀𝑟,^⇀𝑟^′∇^′⋅𝜅^⇀𝑟^′⇀𝑓_𝑗^⇀𝑟^′𝑑𝑣^′,𝐼⁽⁵⁾_𝑖𝑗=_𝜕𝑉𝑖𝑑𝑠𝑚_𝑉𝑖⋅^⇀𝑓_𝑖^⇀𝑟×_𝜕𝑉𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟^′𝑚_𝑉𝑗⋅^⇀𝑓_𝑗^⇀𝑟^′𝑑𝑠^′,𝐼⁽⁶⁾_𝑖𝑗=_𝜕𝑉𝑖𝑑𝑠𝑚_𝑉𝑖⋅^⇀𝑓_𝑖^⇀𝑟×_𝑉𝑗𝑔^⇀𝑟,^⇀𝑟^′∇^′⋅𝜅^⇀𝑟^′⇀𝑓_𝑗^⇀𝑟^′𝑑𝑣^′, (10a)(10b)(10c)(10d)(10e)(10f)where 𝜕𝑉_𝑖 is the outer surface of  𝑉_𝑖 and  𝑚_𝑉𝑖 represents the unit outward normal vector of 𝜕𝑉_𝑖.

Usually, 𝑉_𝑖 consists of two tetrahedron subdomains; that is, 𝑉_𝑖=𝑉⁺_𝑖∪𝑉⁻_𝑖. (11) In such a case, 𝑉_𝑖 is called a complete support set, and 𝜕𝑉_𝑖 consists of 6 triangle subsurfaces. It is readily known that if 𝑉_𝑖 is complete, all the integrals in (10a)–(10f) having 𝜕𝑉_𝑖 as their integration regions are equal to zero. When 𝑉_𝑖=𝑉⁺_𝑖 or 𝑉_𝑖=𝑉⁻_𝑖, 𝑉_𝑖 is called an incomplete support set or a half support set. Clearly, when 𝑉_𝑖 is incomplete, 𝜕𝑉_𝑖 consists of 4 triangle subsurfaces. In such a case, every integral in (10a)–(10f) having 𝜕𝑉_𝑖 as its integration region can be reduced as one on the triangle sub-surface 𝑆_𝑖 of 𝜕𝑉_𝑖 (see Figure 1).

It can be found that every integral in (10a)–(10d) and (10f) for a fixed element 𝑎_𝑖𝑗 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 both 𝑉_𝑖 and 𝑉_𝑗 are complete; then 𝐼⁽²⁾_𝑖𝑗=𝐼⁽²⁾_𝑖+,𝑗++𝐼⁽²⁾_{𝑖+,𝑗−}+𝐼⁽²⁾_{𝑖−,𝑗+}+𝐼⁽²⁾_{𝑖−,𝑗−}, (12) where 𝐼⁽²⁾_{𝑖±,𝑗±} has 𝑉^±_𝑖 and 𝑉^±_𝑗 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 𝐼⁽²⁾_𝑖+,𝑗+=𝑘²₀𝑆_𝑖3𝑉⁺_𝑖𝑆_𝑗3𝑉⁺_𝑗×_{𝑉+𝑖𝑉+𝑗}^⇀𝑟−^⇀𝑟⁺_𝑖,free⋅^⇀𝑟^′−^⇀𝑟⁺_𝑗,free×𝜅^⇀𝑟^′𝑔^⇀𝑟,^⇀𝑟^′𝑑𝑣𝑑𝑣^′. (13) After expanding (^⇀𝑟−^⇀𝑟⁺_𝑖,free)⋅(^⇀𝑟^′−^⇀𝑟⁺_𝑗,free), we obtain 𝐼⁽²⁾_𝑖+,𝑗+=𝑘²₀𝑆_𝑖𝑆_𝑗9×𝑄₁⁽²⁾−^⇀𝑟⁺_𝑗,free⋅𝑄^⇀₂⁽²⁾−^⇀𝑟⁺_𝑖,free⋅𝑄^⇀₃⁽²⁾+^⇀𝑟⁺_𝑖,free⋅^⇀𝑟⁺_𝑗,free𝑄₄⁽²⁾, (14) where𝑄₁⁽²⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖^𝑉+_𝑗^⇀𝑟⋅^⇀𝑟^′𝜅^⇀𝑟^′𝑔^⇀𝑟,^⇀𝑟^′𝑑𝑣𝑑𝑣^′,𝑄₂⁽²⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖^𝑉+_𝑗^⇀𝑟𝜅^⇀𝑟^′𝑔^⇀𝑟,^⇀𝑟^′𝑑𝑣𝑑𝑣^′,𝑄₃⁽²⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖^𝑉+_𝑗^⇀𝑟^′𝜅^⇀𝑟^′𝑔^⇀𝑟,^⇀𝑟^′𝑑𝑣𝑑𝑣^′,𝑄₄⁽²⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖^𝑉+_𝑗𝜅^⇀𝑟^′𝑔^⇀𝑟,^⇀𝑟^′𝑑𝑣𝑑𝑣^′. (15a)(15b)(15c)(15d)

Note that the formulae (15a)–(15d) are only related with the global numbers of the tetrahedra 𝑉⁺_𝑖 and 𝑉⁺_𝑗. 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 𝐼⁽¹⁾_𝑖𝑗, 𝐼⁽²⁾_𝑖𝑗, and 𝐼⁽⁴⁾_𝑖𝑗 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 𝑧-axis and the polarization direction in the 𝑥-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 is 0.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.

Figure 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 Figure 3(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).

(a) The cuboid

(b) The cube

Figure 3 

The partitioned homogeneous cuboid and the inhomogeneous cube.

Example C (inhomogeneous cube). An  inhomogeneous cube is considered in Figure 3(b). The relative permittivity of the cube varies according to 𝜀_𝑟^⇀𝑟=𝑥+𝑦+𝑧9+1.5. (16) 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 𝑍RE=^Benchmark−𝑍^AM_𝐹𝑍^Benchmark_𝐹, (17) where ‖⋅‖_𝐹 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 Table 2. 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.

It 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 as 𝑂(16𝑛²) bytes when using double precision, where 𝑛 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.

Appendix

A. Calculation Formulae

In the main text, the calculation formulae of 𝐼⁽²⁾_𝑖+,𝑗+ have been presented; see (14) and (15a), (15b), (15c), and (15d). Here, the calculation formulae of 𝐼⁽¹⁾_𝑖+,𝑗+ and 𝐼⁽³⁾_𝑖+,𝑗+~𝐼⁽⁶⁾_𝑖+,𝑗+ will be provided.

A.1. The Formula of 𝐼⁽¹⁾_𝑖+,𝑗+ When 𝑉⁺_𝑖=𝑉⁺_𝑗

Consider 𝐼⁽¹⁾_𝑖+,𝑗+=𝑆_𝑖𝑆_𝑗9𝑄₁⁽¹⁾−^⇀𝑟⁺_𝑖,free+^⇀𝑟⁺_𝑗,free⋅𝑄^⇀₂⁽¹⁾+^⇀𝑟⁺_𝑖,free⋅^⇀𝑟⁺_𝑗,free𝑄₃⁽¹⁾, (A.1)

where𝑄₁⁽¹⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖𝜀_𝑟^−1⇀𝑟^⇀𝑟⋅^⇀𝑟𝑄𝑑𝑣,₂₍₁₎=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖𝜀_𝑟^−1⇀𝑟^⇀𝑟𝑄𝑑𝑣,₃₍₁₎=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖𝜀_𝑟^−1⇀𝑟𝑑𝑣. (A.2a)(A.2b)(A.2c)

A.2. The Formula of 𝐼⁽³⁾_𝑖+,𝑗+

Consider 𝐼⁽³⁾_𝑖+,𝑗+=𝑆_𝑖𝑄⁽³⁾, (A.3)

where 𝑄⁽³⁾=1𝑉⁺_𝑖𝑉⁺_𝑖𝑑𝑣_𝑆𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟𝑑𝑠^′. (A.4)

A.3. The Formula of 𝐼⁽⁴⁾_𝑖+,𝑗+

Consider 𝐼⁽⁴⁾_𝑖+,𝑗+=𝑆_𝑖𝑆_𝑗𝑄₁⁽⁴⁾−^⇀𝑟⁺_𝑗,free⋅𝑄^⇀₂⁽⁴⁾, (A.5)

where𝑄₁⁽⁴⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖×𝑑𝑣_𝑉+𝑗𝑔^⇀𝑟,^⇀𝑟^′×𝜅^⇀𝑟^′+13∇^′𝜅^⇀𝑟^′⋅^⇀𝑟^′𝑑𝑣^′,𝑄₂⁽⁴⁾=1𝑉⁺_𝑖𝑉⁺_𝑗^𝑉+_𝑖𝑑𝑣_𝑉+𝑗𝑔^⇀𝑟,^⇀𝑟^′13∇^′𝜅^⇀𝑟^′𝑑𝑣^′. (A.6a)(A.6b)

A.4. The Formula of 𝐼⁽⁵⁾_𝑖+,𝑗+

Consider 𝐼⁽⁵⁾_𝑖+,𝑗+=_𝑆𝑖𝑑𝑠_𝑆𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟^′𝑑𝑠^′. (A.7)

A.5. The Formula of 𝐼⁽⁶⁾_𝑖+,𝑗+

Consider 𝐼⁽⁶⁾_𝑖+,𝑗+=𝑆_𝑗𝑄₁⁽⁶⁾−^⇀𝑟_𝑗,free⋅𝑄^⇀₂⁽⁶⁾, (A.8)

where𝑄₁⁽⁶⁾=1𝑉⁺_𝑗^𝑆_𝑖×𝑑𝑠_𝑉+𝑗𝑔^⇀𝑟,^⇀𝑟^′𝜅^⇀𝑟^′+13∇^′𝜅^⇀𝑟^′⋅^⇀𝑟^′𝑑𝑣^′,𝑄₂⁽⁶⁾=1𝑉⁺_𝑗^𝑆_𝑖𝑑𝑠_𝑉+𝑗𝑔^⇀𝑟,^⇀𝑟^′13∇^′𝜅^⇀𝑟^′𝑑𝑣^′. (A.9a)(A.9b)

Acknowledgment

This work was supported by the National Basic Research Program of China (nos. 2009CB320203 and 2010CB327400).