033256 -L- T A HYBRID FINITE ELEMENT-FAST MULTIPOLE TECHNIQUE FOR ELECTROMAGNETIC SCATTERING S.S. Bindiganavale and J.L. Volakis Naval Air Warfare Center Weapons Div. China Lake, CA 93555-6001 January 1996 33256-2-T = RL-2456

Report #033762-1-T Contract Title: Further development of FEMATS including prismatic meshes, graphical interface and new mesh truncation schemes Contract No.: Report Title: Report Authors: N68936-95-C-0349 Primary University Collaborator: Contract monitor: University Address: A hybrid FEM-FMM technique for electromagnetic scattering Sunil S. Bindiganavale and John L. Volakis John L. Volakis Dr. Helen Wang Code 455520D Bldg 31567 Naval Air Warfare Center Weapons Division China Lake, CA 93555 Radiation Laboratory Department of Electrical Engineering and Computer Science Ann Arbor, MI 48109-2122 Email: volakis@umich.edu Phone: (313)764-0500 Date: January 1996

A hybrid FEM-FMM technique for electromagnetic scattering Sunil S. Bindiganavale and John L. Volakis Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Abstract In this paper, we apply a version of the Fast Multipole Method (FMM) to reduce the storage and computational requirement of the boundary integral in the finite element-boundary integral method. By virtue of its 0(N1'5) operation count, the application of the singlestage FMM, results in substantial speed-up of the boundary integral portion of the code, independent of the shape of the BI contour. We will discuss the efficiency of the method and present the application of this technique to the computation of electromagnetic scattering from large grooves recessed in a ground plane. 1 Introduction Over the past few years different hybrid versions of the finite element method have been explored for application to scattering by composite structures. Among them, the finite element-boundary integral equation (FE-BI) and the finite element-absorber boundary condition (FE-ABC) methods have been quite popular and extensively applied to many applications. The FE-BI method [1],[2] employs the exact boundary integral equation which provides 1

an independent relation between the tangential E and H fields on the mesh outer boundary and is therefore an exact method. This is in contrast to the FE-ABC method [3] which employs an approximate truncation operator but leads to fully sparse systems. On the other hand, the FE-BI method, although "exact", leads to a partly full and partly sparse system (see figure 1) and is thus more computationally intensive. For special cases, where the boundary is rectangular or circular, the FFT can be used [4] to reduce the memory and CPU requirements to N log N. However, in general, the boundary integral is not convolutional and in that case the CPU requirements will be of O(Nb) where Nb denotes the unknowns on the boundary. In this paper, we apply a version of the Fast Multipole Method (FMM) to reduce the storage and computational requirement of the boundary integral when the size of the contour is large. By virtue of its 0(N1 5) operation count, the application of the FMM, results in substantial speed-up of the boundary integral portion of the code, independent of the shape of the BI contour. We consider the application of this technique (referred to hereon as the FEM-FMM method) to scattering from a material-filled groove in a conducting plane. The FEM is employed to formulate the fields within the cavity and establish a relationship with those at the aperture. The fields external to the groove are expressed as an integral over the aperture and a system of integral equations is then obtained by enforcing field continuity across the aperture. To reduce the storage requirement and speed up the computation of the boundary integral, an approximate version of the fast multipole method [5] is employed. 2 Formulation Consider a filled PEC groove of width w and depth d in an otherwise uniform ground plane as shown in figure 2(a). The material filling the groove has a permittivity cr and permeability,u The free space region exterior to the groove is denoted as region I while the groove itself is denoted as region II (see figure 2). The fields in the two regions are decoupled by closing the aperture with a perfect conductor and introducing an equivalent magnetic current based on the equivalence principle M1 = E1x (1) 2

where E1 is the electric field at the aperture. The ground plane can be removed by application of image theory and hence the field in region I can be expressed as the radiation caused by M1 and the external sources (Ji, Mi). The coupling of the fields in each region is achieved by requiring continuity of the tangential magnetic field across the aperture HIaJy = HI, l (2) Htan y=0 = Htan y=0 (2) The magnetic field, HII, inside the groove is formulated via the finite element method which has the inherent geometrical and material adaptability and low O(N) storage requirement. This procedure and the coupling in (2) is described in detail in [1]. HI is expressed as an integral of M1 using the free space Green's function. To this integral equation we apply a version of the fast multipole method [5] to reduce the operation count from O(Nb) to O(Nb'5) where Nb is the number of unknowns on the boundary; the next section describes this procedure for TE incidence. 2.1 Boundary integral for TE case The groove is illuminated by the plane wave Hi- = ejko(x sin o +y cos o) (3) where ko = 2x7/Ao is the free space wavenumber and qo is the angle of incidence. With a z-directed impressed magnetic field Hi, the scattered magnetic field will also be z-directed, and consequently the equivalent magnetic current M1 may be written as Ml = zM1(x) (4) From figure 2(b), the magnetic field in region I due to M1 is given by HI(r) = Hnc(x,y) + Hefl(x, y) - jkoYo 10 Ml(x')Go(p, p')dx' (5) where Go(P) = j{ ) (ko(x - 2 + (y - y) } (6) 3

is the free space Green's function. Applying Galerkin's technique to (5) yields the matrix equation [c]{m} = {nc} - [Pp]{~} (7) where {q} is the column vector with Nb + 1 elements representing the nodal magnetic field on C1, and Nb is the number of segments employed for the discretization of C1. {b} is the column vector with Nb elements representing the magnetic current on the boundary. The column vector {4nC} has Nb elements given by in' = (HtIC(Xm, ) + Hz;efl(x 0)) Am, m = 1, 2,..., Nb (8) with AXm being the segment length of the mth segment on C1 and Xm being its midpoint. The product [C]{+q} in (7) is of not much concern since [C] is a sparse matrix assembled from the interaction of the magnetic field and current on the boundary segments. The product [P]{} is of concern because of the 0(Nb2) storage required by the full matrix [P] whose elements are given by Pmn = 2 {H02 (ko IXm - a) AmAXn} m n (9) To employ the approximate version of the FMM in solving (7), the unknowns on the boundary are divided into groups with Mb unknowns in each group and thus the number of groups will then be -Nb. For large source to observation Mb e distances, the kernel in (9) is approximated by using the large argument expansion as [5] 2 2, F Ce Ax31'1 H(2)(k0 Ixm - x I) - e-jkorxlm -e-koxnl (10) where xltl is the distance between the center of the test group I and the center of the source group 1/; Xni' is the distance between the nth source element and its group center and xlm is the distance between the mth test element and its group center. The decoupling of test-source element interactions in the kernel as in (10) enables the computation of the matrix-vector product for far-field groups with a reduced operation count as detailed in the following sequence. 4

1. For each test group, the aggregation of source elements in a single source group involves Mb operations, corresponding to the number of elements in the source group. The aggregation operation corresponds to Mb bll= Mjeikoxn (11) j=l 2. Since the above aggregation operation needs to be done for all source groups the operation count becomes O(NbMb) ~ O(Nb), where -N Mb represents the total number of groups. Also this operation, being dependent only on the test group rather than the test element, needs to be repeated for Nb_ test groups leading to a total operation count of Mb o(:) for aggregation. 3. The next step would be a translation operation corresponding to -2jkzoxj cI =,-b (12) V7 \/kox Since this needs to be done only at the group level, it involves 0( -) operations for all possible test and source group combinations and is the least computationally intensive step. 4. The final step in the sequence would be the process of disaggregation corresponding to the operation Nb/Mb Iln= CLle-jkoixm (13) 1/=l Conceptually, this process is the converse of aggregation. Since this operation involves only the source group instead of the source element it needs to be done for each source group thus implying an 0(-b) operation to generate a single row of the matrix-vector product. To generate Mb rows corresponding to a test group the operation count would be O(Nb). With -b test groups, the operation count would be of ()Mb 5

5. The near field operation count being of O(NbMb) and the far field being N2 O( -) gives a total operation count of Mb Op.count - C1Nb Mb + C2 N (14) Typically, we can set the elements in each group, Mb = v/Nb and as a result the total operation count is O0 N15. Once the currents are computed in this manner, the far-fields are computed in th usual manner. Thus the maximum computational requirement in this hybrid algorithm is O(N1'5) unlike the usual FE-BI which results in O(N2) if the boundary has N unknowns. 3 Results A computer code based on the above formulation was implemented and executed on a HP 9000/750 workstation with a peak flop rate of 23.7 MFLOPS. The geometry considered was the rectangular groove shown in figure 1. Figure 3 shows the bistatic RCS at normal incidence for a groove 20A wide and 0.25A deep with a material filling of,r = 4 and /,r = 1. The discretization of this geometry results in 2106 unknowns with a boundary integral of 300 unknowns. To obtain a measure of the speed-up realized by the application of the FEM-FMM to this problem, we refer to the CPU time comparisons also shown in figure 3. The compiled data on the FE-BI and FEM-FMM methods include the CPU time, storage requirement and average error. Data are also give for the FE-BI implementation using the special CG-FFT solver but it should be noted that the application of the CG-FFT is only suited for planar apertures. Clearly, the compiled run data show that the FEM-FMM is about three times faster than the traditional FE-BI implementation and requires only one sixth of the memory. However, the FEM-FMM has some deterioration in accuracy (average error of 2 dB for this example [6]) and this should be considered depending on the application and usage of the results. Similar data shown in figure 4 for plane wave incidence at an off-normal incidence angle (450) seem to verify the previous conclusions. The groove analyzed in figure 4 is 30A wide and is 0.15A deep with a material filling of Or = 6 and fi = 1. The discretization of this geometry results in 1803 unknowns with a 6

boundary integral of 450 unknowns. The error in the FEM-FMM solution is lower but this could be attributed to the use of a larger near-group window (4.3A instead of 2.45A) which results in an increase in storage and execution time of the boundary integral. Clearly the CG-FFT solver is substantially faster and arguably uncontested by the FEM-FMM method. This simply gives credence to modifications of the boundary integral surface so that the CG-FFT can be applied whenever possible even though the boundary surface may not be completely planar. References [1] J.M. Jin and J.L. Volakis, "TE scattering by an inhomogeneously filled aperture in a thick conducting plane," IEEE Trans. Antennas Propagat., vol.38, pp.1280-1286, Aug. 1990. [2] J.M. Jin and J.L. Volakis, "TM scattering by an inhomogeneously filled aperture in a thick conducting plane," IEE Proc., pt. H, vol.137, pp. 153-159, June 1990. [3] A. Chatterjee, J.M. Jin, and J.L. Volakis, "Edge-based finite elements and vector ABCs applied to 3D scattering," IEEE Trans. Antennas Propagat., vol. 41, pp.221-226, Feb 1993. [4] J.D. Collins, J.L. Volakis, and J.M. Jin, "A combined finite elementboundary element formulation for solution of two-dimensional scattering problems via CGFFT," IEEE Trans. Antennas Propagat., vol. 38, pp. 1852-1858, 1990. [5] C.C. Lu and W.C. Chew, Fast far field approximation for calculating the RCS of large objects, Micro. Opt. Tech. Lett., 8(5):238-241, April 1995. [6] S.S. Bindiganavale and J.L. Volakis, Guidelines for using the fast multipole method to calculate the RCS of large objects, Micro. Opt. Tech. Lett., March 1996. 7

A y BI A x FEM n. 20 40 3. p4:. 'k'4, = th %. ' ' '. '11. ':.oo - *.....::. '%:. ' *:': %~ *%. '*i 1:.. ~'~. * *-*.~ *- *. H~ ~s~ s.:S.s. *s.ss ~~ s hSSs. *ss *sss **~ *. I. ~e. BI System 1. o 60 80 100 l:..U IL 'tL IL IL IL -1 OnI I/l. 0 20 40 60 column 80 100 120 Figure 1: The FE-BI system matrix arising from the scattering/radiation problem of a groove in a ground plane. Unknowns 88-105 are the boundary integral unknowns. 8

I A y Region 1 ) C A x w Region 2 (a) J. M., M 1- V Region 1 y=0 ground plane (b) () r Region2 1 ( 6r 9 ^t d Region 2 1 (c) Figure 2: (a) Geometry of the groove recessed in a ground plane (b) Equivalent problem for region 1, and (c) region 2 9

-b c) cr 0 o rn 0 10 20 30 40 50 60 70 80 90 Observation angle (deg.) I,, L - Time (for BI only) FEM-BI (CG) (Minutes,seconds) 15 mins 14 sees FEM-FMM 4 mins 48 sees - FEM-BI (CGFFT) 1 min 40 sees Memory (for BI only) LU FEM-FMM (KB) 680 116 - Average (dB) - error FEM-FMM 2.016 FEM-BI (CGFFT) 0.98 Figure 3: Bistatic scattering from a rectangular groove 20A wide and 0.35A deep at normal incidence 10

co -2 10 0 20 40 60 80 100 120 140 160 180 Observation angle (deg.) Time (for BI only) FEM-BI (CG) | FEM-FMM FEM-BI (CGFFT) I [(Minutes,seconds) 30 mins 36 secs 9 mins 1 sec 2 min 15 secs Memory (for BI only) LU FEM-FMM (KB) 1544 436 Average error FEM-FMM FEM-BI (CGFFT) (dB) 1.553 1.0267 Figure 4: Bistatic scatterin ro a ng from a rectangular groove 30 wide and 0.15 deep at 450 incidence 11