Fast Integral Equation Solution of Scattering of Multiscale Objects by Domain Decomposition Method with Mixed Basis Functions

. Nonconformal nonoverlapping domain decomposition method (DDM) with mixed basis functions is presented to realize fast integral equation solution of electromagnetic scattering of multiscale objects. The original multiscale objects are decomposed into severalclosedsubdomains.Thehigherorderhierarchicalvectorbasisfunctionsareusedintheelectricallylargesmoothsubdomains tosignificantlyreducethenumberofunknowns,whiletraditionalRao-Wilton-Glissonbasisfunctionsareusedforsubdomains withtinystructures.Awell-posedmatrixissuccessfullyderivedbythepresentDDM.Besides,thenonconformalpropertyofDDM allowsflexiblemeshgenerationforcomplicatedobjects.Numericalresultsarepresentedtovalidatetheproposedmethodand illustrateitsadvantages.


Introduction
As an effective full wave method, integral equation method (IEM) is one of the most widespread methods for solving electromagnetic (EM) scattering and radiation problems.Different from partial differential equations, integral equations satisfy the radiation condition automatically, and their numerical solution does not require an absorption boundary condition.Besides, the unknowns of IEM are only distributed on surface for perfect electric conductor (PEC), and the number of unknowns is much less than the one by finite element method and finite difference method.However, for electrically large structures, traditional IEM with low order basis like Rao-Wilton-Glisson (RWG) basis [1] still leads into a large number of unknowns, which makes it challenging to store and solve the matrix equation.To circumvent this difficulty, many fast methods were developed, including the multilevel fast multipole algorithm (MLFMA) [2], the adaptive integral method [3], integral equation fast Fourier transformation (IE-FFT) [4], hierarchical matrix methods [5,6], and its parallel version [7].
These fast methods greatly reduced the computational complexity and memory requirements.Another kind of fast methods is to reduce the dimension of the matrix by using higher order basis functions.Higher order basis functions can be categorized into two kinds, the interpolatory type and the hierarchical type [8][9][10].The higher order interpolatory basis functions interpolate the value of field on a large patch with many interpolation points.The limitation of this basis is the expansion order which must be kept constant for all patches.Besides, the interpolatory type basis requires conformal mesh.The higher order hierarchical vector (HOHV) basis function does not need constant expansion order on different patches.It defines both low order and higher order basis on the same patch, and actually the basis of order  is a subset of basis of order  + 1.
With the HOHV basis function, a well-conditioned matrix can be obtained for conventional structures [9].However, when modeling multiscale structures, the high contrast (electrically large and electrically small) sized meshes exist together, so the proper equilibration technique of MoM matrix equation is needed to decrease the condition number 2 International Journal of Antennas and Propagation [11].At the same time, ill-shaped meshes are inevitable and the orthogonality of HOHV basis function is affected.This makes the matrix equation ill-conditioned and causes convergence difficulty for Krylov subspace methods [12].Therefore, HOHV basis function cannot be directly applied to model multiscale structures.
Recently, a novel nonconformal, nonoverlapping domain decomposition method (DDM) based on IEM has been developed by Peng et al. successfully [13].This IE-DDM is based on the philosophy of "divide and conquer." It divides original problems into several closed subdomains and enforces transmission conditions (TCs) on touching faces to maintain the continuity of currents across interfaces.Because of the nonconformal property of IE-DDM, each subdomain can be meshed independently, and high mesh quality can be easily realized in subdomains.In addition, IE-DDM also provides an effective preconditioner and makes the system matrix for multiscale problems well-posed.
In this paper, mixed basis functions are incorporated into the framework of IE-DDM (named as IE-DDM/MBFs) to further enhance the ability of IE-DDM for multiscale objects.Here, higher order hierarchical vector basis functions defined over quadrilaterals are used to expand surface current on electrically large smooth surfaces, and the traditional RWG basis functions defined over triangles are used to expand surface current on electrically small structures.Compared to the traditional higher order method, the present method realizes flexible mesh generation and fast convergence for multiscale EM problems.
This paper is organized as follows: the derivation of IE-DDM/MBFs is presented in Section 2. Numerical examples of EM scattering from PEC structures are presented in Section 3. Finally, conclusions are drawn in Section 4.

Derivation of IE-DDM/MBFs
Consider EM scattering of a PEC object Ω. Γ is the exterior surface of Ω.According to the surface equivalence principle, the surface equivalent current is represented as  = n ×  tol , where  tol =  inc +  sca .Based on the relationship between current source and scattering field, the scattered field can be represented as follows: In (1), (,   ) = exp(− 0 )/4 is the scalar Green function in free space, where  = | −   | is the distance between the source point and field point. 0 is the intrinsic impedance in free space. 0 is the wave number in free space.By using electric field boundary condition or magnetic field boundary condition, we can derive the electric field integral equation (EFIE) or the magnetic field integral equation (MFIE).Based on the linear combination of EFIE and MFIE, the combined field integral equation (CFIE) can be expressed as Here,  is the principle value of ,  is the combination factor, and  = 1 − .As shown in Figure 1, IE-DDM decomposes the original domain Ω into two closed nonoverlapping subdomains Ω 1 and Ω 2 , with Ω = Ω 1 ∪ Ω 2 .In Figure 1, Γ 1 and Γ 2 are the exterior surfaces of Ω 1 and Ω 2 , respectively.Furthermore, surfaces Γ 1 and Γ 2 can be decomposed as and Γ + 1 and Γ − 2 are the touching faces between Ω 1 and Ω 2 .As shown in Figure 2,  + 1 and  − 2 are the surface current on touching faces Γ + 1 and Γ − 2 , respectively.The decomposed problem is equivalent to the original problem by enforcing the transmission condition on touching faces: A remarkable advantage of IE-DDM is nonconformal mesh.
In other words, mesh density and shape can be chosen according to the geometrical characteristic of each subdomain.This also means mixed basis functions can be used to make the choice of basis functions compatible with each subdomain.Based on IE-DDM/MBFs, well-conditioned system matrices are attained in all subdomains, which leads to fast and stable convergence for multiscale objects.In this paper, mixed basis functions are firstly incorporated into the IE-DDM to form the IE-DDM/MBFs method.Here, the HOHV basis functions [9] defined over quadrilaterals are used to expand surface current on electrically large smooth surfaces, and traditional RWG basis functions [1] defined over triangles are used to expand surface current on electrically small structures.For example, the HOHV basis functions  HO based on quadrilateral surfaces are used to discretize the surface current on Ω 1 , and the RWG basis functions  RWG are used to discretize the current on Ω 2 .By applying Galerkin methods, IE-DDM/MBFs matrix equations can be derived as where In ( 5) and ( 6 1 are the cement matrices enforcing the transmission condition.Detailed expression of these matrix blocks is International Journal of Antennas and Propagation where ⟨, ⟩ denotes the inner product of two vector functions calculated as follows: For nonconformal IE-DDM, it is important to compute the cement matrix fast and accurately.For nonconformal meshes on the touching faces between adjacent subdomains, the union mesh technique [14,15] is adopted, as shown in Figure 2. The procedure to compute the cement matrix is described as follows. (a) Determine the common region between two nonconformal overlapping patches, triangle patch   and quadrilateral patch   , shown in Figure 2. A fast algorithm to determine the common region is critical to achieve high efficiency.Here, an octree algorithm is developed to determine the common region.The detailed process of octree can be descripted as follows.
(1) Enclose the original decomposed objects with a large cube and then partition the cube into eight smaller uniform sized cubes.Decompose each subcube into eight smaller uniform sized cubes and implement this process recursively until the length of the fines cube is about 0.5 .(2) Group the basis functions on the touching faces into the corresponding cubes and sign the corresponding boxes.The connected triangles on the touching faces can be found in each box itself and neighbor boxes.Based on this octree algorithm, the complexity of determining the common region is only (log), while the one of determining the common region by direct technique is ( 2 ), where  is the total number of patches on touching faces.
(b) Computation of cement matrices: the polygon shown in Figure 2 is first divided into many small triangles, and then the Gaussian quadrature formula is used over each triangle.The integration results are finally attained based on the summation of the integration result over each triangle.
The comparison of the CPU time required for determining the common region by the octree algorithm and direct technique is shown in Figure 3, in which the number of patches on the touching faces varies from 1,000 to 13,000.Obviously, the octree algorithm is very attractive, especially with the unknown of patches dramatically increasing.To realize a fast iterative solution of (4), a preconditioner is often required.Here, the inverse of matrix  is used as the preconditioner.Equation ( 4) is then written as follows: Because it is difficult to directly calculate the inverse of matrix , especially for electrically large problems, the innerouter iteration based on Krylov subspace method can be used to solve (10).In this paper, the generalized minimal residual (GMRES) [12] method is used in both inner and outer iterations.The convergence tolerances of outer and inner iterations are 0.01 and 0.001, respectively.

Numerical Results
In this section, numerical examples of EM scattering from PEC objects are presented to prove the validity of the present method.To illustrate its efficiency, the proposed method is also compared with traditional MoM with RWG and HOHV basis functions.

Example for Validation.
In the first numerical example, EM scattering from a PEC sphere is investigated to demonstrate the validity of IE-DDM/MBFs.The radius of sphere is 1.5 m.The incident wave illuminates the sphere from the negative direction of the -axis, and its frequency is 0.3 GHz.The sphere is divided into two subdomains, whose dimensions and meshes are shown in Figure 4.The RWG basis functions with mesh size of 0.1  are used in the first region, and the second-order HOHV basis functions with mesh size of 0.3  are used in the second region.The bistatic radar cross section (RCS) for H-H polarization is shown in Figure 5.It is obvious that the result of IE-DDM/MBFs agrees well with that of MIE series.

Efficiency of IE-DDM/
MBFs.In the second numerical example, EM scattering from a PEC circular cone with a tiny structure is investigated.The circular cone is illuminated by a plane wave of 1.2 GHz from the nose direction.The RWG basis functions with mesh size of 0.1  are used in the first subdomain, and the second-order HOHV basis functions with mesh size of 0.3  are used in the second subdomain.
The dimension and mesh of these two regions are shown in Figure 6.The computational cost statistics of IE-DDM/MBFs, traditional MoM with RWG basis functions (MoM-RWG), and IE-DDM with RWG basis functions (IE-DDM/RWGs) accelerated by the MLFMA are shown in Table 1.For this multiscale structure, the total number of unknowns required by RWG basis function is 31,950, and the one by IE-DDM/RWGs is 24,213.On the other hand, the total number of unknowns by IE-DDM/MBFs is only 8,712.The reduction of number of unknowns is based on the nonconformal DDM framework and the property of higher order basis functions.The computational advantage of the IE-DDM/MBFs over traditional MoM is also remarkable.As shown in Table 1, the total time by IE-DDM/MBFs is only 18 minutes, while MoM requires 25 times longer CPU time than IE-DDM/MBFs.It should be mentioned that, due to the MLFMA accelerating  the IE-DDM/RWGs, its CPU time is less than that of IE-DDM/MBFs.Figure 7 presents the bistatic RCS results calculated by three methods.The polarization is H-H and the frequency is 1.2 GHz.It is seen that the results of the proposed method agree well with the one of IE-DDM/RWGs, MoM-RWG methods.In the third numerical example, EM scattering from a simplified PEC ship model is investigated.The incident plane wave illuminates the model along the negative direction of -axis, and its frequency is 0.6 GHz.The dimension and geometry of the ship model is shown in Figure 8.This ship model is divided into 4 closed sub-domains, whose shapes and meshes are shown in Figure 9.
The second-order HOHV basis functions with the maximum mesh size of 0.5  are used in the first subdomain, and the RWG basis functions with the minimum mesh size of 0.01  are used in the second, third, and fourth subdomains.The numbers of unknowns for the first to the fourth subdomains are 9,712, 2,220, 3,801, and 5,688, respectively.Traditional higher order MoM can also be used to solve this problem, but the system matrix is ill-conditioned due to the multiscale property of the geometry.The maximum and minimum sizes of HOHV basis functions are 0.5  and 0.01 , respectively.For the meshes with size less than 0.1 , the first-order basis functions are used.The secondorder basis functions are used for the other meshes.Detailed computational cost statistics of IE-DDM/MBFs, MoM-RWG, and higher order MoM (HO-MoM) are shown in Table 2.For this complicated structure, the higher order MoM requires 2420 iterations to converge to relative residual error of 0.01.But the preconditioning IE-DDM only requires 9 outer iterations to converge.The iteration time by HO-MoM and IE-DDM/MBFs is 67 minutes and 8 minutes, respectively.Therefore, compared to HO-MoM, the proposed method achieves significant reduction of iteration time.Although it is possible to improve the iterative convergence by adopting 6 International Journal of Antennas and Propagation   preconditioner technique in higher order MoM, it is tedious and time-consuming to construct a good preconditioner for higher order MoM when modeling multiscale objects.For such multiscale structure, traditional MoM-RWG leads into a large number of unknowns.As shown in Table 2, the number of unknowns required by HO-MoM and IE-DDM/MBFs is 25,034 and 21,421, respectively, while MoM-RWG requires 89,739 unknowns.Obviously, MoM-RWG requires large memory and long CPU time.Figure 10 presents the H-H polarized bistatic RCS results of the simplified ship model from different methods.Good agreement is observed between the proposed method and other methods.

Conclusions
The integral equation domain decomposition method with mixed basis functions is developed to solve time harmonic multiscale EM scattering problems.Based on the nonconformal, nonoverlapping domain decomposition method, the original problem is decomposed into several closed subdomains in which different basis functions can be used according to the geometry characteristics of subdomains.By using IE-DDM/MBFs, the multiscale problems have been successfully solved.Although mixed higher order hierarchical vector basis functions and RWG basis functions are used here, the present method can be extended to incorporate other basis functions.

Figure 1 :
Figure 1: The notation of the decomposition of the original PEC object.

Figure 2 :
Figure 2: (a) An example of the nonconformal patches on the interface.(b) The description of evaluation of the cement matrix.

Figure 3 :
Figure 3: The CPU time of determining the common regions of meshes with octree algorithm and direct technique.

Figure 6 :
Figure 6: The dimension and the mesh of the circular cone.

Figure 7 :
Figure 7: The bistatic RCS of the circular cone.

Figure 8 :
Figure 8: The geometry and dimension of the simplified ship model.

Figure 9 :Figure 10 :
Figure 9: Domain partition by using IE-DDM/MBFs to solve the EM scattering of simplified ship model.

Table 1 :
Computational statistics of the circular cone.

Table 2 :
Computational statistics of the ship.