035067-12-T ADAPTIVE INTEGRAL METHOD FOR HYBRID FE-BI MODELING OF 3D DOUBLY PERIODIC STRUCTURES Sanders, A Lockheed Martin Co. 95 Canal Street NCA1-6268 P.O. Box 868 Nashua, NH 030601-0868 Thomas F. Eibert John L. Volakis July 1998 35067-12-T = RL-2492

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: U-M REPORT No.: CONTRACT START DATE: END DATE: DATE: SPONSOR: SPONSOR CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: Hybrid Finite Element Design Codes for the SERAT Array Adaptive Integral Method for Hybrid FE-BI Modeling of 3D Doubly Periodic Structures 035067-12-T October 1996 September 1998 July 28, 1998 Roland Gilbert SANDERS, INC, A Lockheed Martin Co. MER 24-1583 PO Box 868 Nashua, NH 030601-0868 Phone: (603) 885-5861 Email: RGILBERT@mailgw.sanders.lockheed.com P.O. QP2047 John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave Ann Arbor, MI 48109-2122 Phone: (313) 764-0500 FAX: (313) 747-2106 volakis @umich.edu http://www-personal.engin.umich.edu/-volakis/ CONTRIBUTORS TO THIS REPORT: T. Eibert, J. Volakis

Adaptive Integral Method for Hybrid FE/BI Modeling of 3D Doubly Periodic Structures Thomas F. Eibert and John L. Volakis Radiation Laboratory, EECS Department The University of Michigan, Ann Arbor, MI 48109-2122 Abstract The adaptive integral method (AIM) is applied in conjunction with a hybrid finite element (FE)/boundary integral (BI) approach for the electromagnetic analysis of three-dimensional (3D) doubly periodic structures. Starting from the conventional mixed potential integral equation (MPIE) formulation, the AIM is derived using Taylor series expansions for the involved coupling integrals. Implementation issues are discussed which are related to the periodic formulation. Central processing unit (CPU) timings and memory requirements are given and compared to the conventional BI implementation. 1 Introduction The application of hybrid FE/BI methods to infinite periodic structures (antennas or frequency selective surfaces) [1, 2, 3] is very attractive because it provides full 3D modeling flexibility. Basically, the FE method is used to model the dielectric section of the unit cell, whereas the BI provides for a rigorous mesh truncation at the upper and/or lower surfaces of the discretized unit cell. Most practical problems can be analyzed using planar BI surfaces and the corresponding half-space periodic Green's functions. However, for large unit cell apertures, the resulting fully populated BI matrix leads to CPU intensive solutions. To alleviate the 1

CPU and memory bottleneck, fast integral algorithms such as the fast multipole method (FMM) [4, 5, 6, 7, 8] or the adaptive integral method (AIM) [9, 10, 11] may be employed. Till now, fast integral methods were mostly used in conjunction with pure integral equation formulations. An application of the FMM to the FE/BI method was presented in [6] and also in [7]. However, for planar BI terminations, AIM is the method of choice because it directly results in low 0(ns) storage and 0(nslogns) CPU time requirements for the execution of the matrix-vector products in the iterative solver (ns: number of BI unknowns). In contrast, FMM must be implemented in its more involved multi-level formulation [8] to achieve the same 0(ns logns) complexity. Also, FMM cannot be directly applied to periodic problems because its concept is based on the free-space Green's function. The speed-up and memory savings of AIM are attained by taking advantage of the Toeplitz properties of the Green's functions when uniform grids are employed. The underlying concept which is based on the convolution theorem and the use of Fast Fourier Transform (FFT) algorithms has been extensively applied in Conjugate Grandient (CG) - FFT like algorithms for many years [12]. Originally these methods were limited due to the inherent necessity for uniform meshes. To obtain more flexible methods, it was necessary to project arbitrary irregular meshes onto a uniform mesh and vice verca. Such a procedure can produce accurate results only for far interactions and therefore near interactions are usually calculated in a conventional manner. The latter can be considered as a correction to the FFT-based fast algorithm. In [13, 14] algorithms, referred to as "Precorrected-FFT Methods", were presented and applied to electrostatics. The projections between the regular and irregular meshes were achieved based on multipole expansions of the potentials as described in [15] and d-functions were used to expand the charges on the regular grid. In [16], a full-wave method was published for the analysis of cavity backed antennas. Both, the uniform and the irregular BI meshes 2

were based on triangular basis functions and the projection operators were generated by simple interpolations between the basis functions on the different meshes. The AIM concept for radiation and scattering calculations was first suggested in [9] and described at some detail in [10]. It employs 5-functions on the uniform grid and generates the projection operators between the uniform and irregular meshes by equating a finite number of multipole moments of the basis functions (not multipole expansions of the Green's function as in the pre-corrected FFT). Even though the authors mention in [10] that other projection methods might be superior to the described matching of moments, we consider this projection algorithm as the characteristic feature of AIM. Its most important advantage as compared to the pre-corrected FFT is that it can be applied directly to any integral equation irrespective of the Green's function as long as the convolution property is retained. Therefore, AIM can be applied to periodic problems in a straightforward manner. In this paper, we employ AIM to reduce the CPU time and memory demand of the BI portion of a hybrid FE/BI method as pertains to 3D doubly periodic structures. The paper focuses especially on issues related to the implementation of the AIM in connection with infinite periodic structures. It begins by presenting an intuitive mathematical derivation of AIM using a Taylor series expansion for the evaluation of the coupling integrals. The latter part of the paper presents CPU speed-up and memory reduction achieved by the AIM implementation when compared to the more conventional simulation of periodic arrays. 3

2 Formulation 2.1 Hybrid FE/BI System The conventional implementation of the hybrid FE/BI method for an infinite periodic structure (ejwt time convention) results in a linear algebraic system of the form Aint Across Across Eint ] 0 Eint fint 1,top JXI b E 0 0 0 E i l i t Across Abd Ebound + 0 Ztop O Ebound boud 2,bot bot bot (bo J bottom to 40 non-zero elements) and are associated with the FE portion of the hybrid method. The BI matrices Z are fully populated and are associated with the boundary edges on the top and bottom boundaries of the discretized periodic cell. The right-hand side vector elements f account for excitations in the FE and BI portions of the method. Mesh termination at the side walls of the unit cell is achieved by imposing the appropriate phase boundary condition (PBC) according to the Floquet theorem. In our implementation, we employed triangular prismatic elements to generate the volume mesh giving rise to triangular surface meshes with Rao-Wilton-Glisson basis functions [17] for the magnetic currents on the planar BI surfaces. For arrays with arbitrary scan angles, all matrices in (1) may be nonsymmetric and this is due to the PBCs and the periodic Green's function. For the solution of the system, the biconjugate gradient (BiCG) or generalized minimal residual (GMRES) solvers are used, both of which rely on an efficient evaluation of the matrix-vector products associated with (1) [18]. Due to the nonsymmetric properties of the system matrix, the employed BiCG solver requires two matrix-vector products for the generation of one search vector whereas the GMRES solver needs only one. However, the GMRES solver has a considerably higher memory demand due 4

to storing the search vectors within the restart loop. Also, additional computational effort is required for the orthogonalization of the search vectors [18]. The computational complexity per matrix-vector product of the total A-matrix is of ((nv), where nV denotes the number of volume edges. However, the complexity in executing the matrix-vector products [Z]{x} is O(n2) and storage requirements are also of the same order. For a more efficient implementation of the method, it is therefore crucial to reduce the complexity of the BI matrix-vector products. We propose to accomplish this by using AIM. 2.2 Derivaton of the AIM Using Taylor Series All elements of AIM together with a detailed error analysis are discussed in [10]. However, for a more intuitive and straightforward derivation of AIM, we propose the introduction of Taylor series expansions for parts of the integrands involved in the MoM coupling integrals. Therefore, in this subsection we first review AIM based on this new method before we discuss its implementation into the BI portion of our hybrid FE/BI approach for periodic structures. The starting point of AIM is to split the Z-matrices as [Z] = [znear] + [Zfar], (2) where [znear] contains the elements of [Z] which are near the self-cell. Correspondingly [Zfar] contains the remaining "far-zone" elements of [Z]. AIM reduces the CPU time and memory demand of the iterative solver by exploiting the convolutional properties of the Green's function for the evaluation of the matrix-vector products associated with the mostly full matrix [Zfar] [10]. That is, the far-zone matrices are not explicitly generated and the matrix-vector products are performed in the discrete Fourier domain (DFT) utilizing appropriate 2D FFT algorithms [12]. However, the convolutional properties can only be exploited on a uniform grid and therefore, the basis functions on the original and possibly irregular triangular mesh must first 5

be replaced by an equivalent expansion using unknowns placed on a uniform grid. To do so, we introduce auxiliary basis functions on a uniform rectangular grid which is coincident with the original triangular grid. The auxiliary expansion for the magnetic surface currents (equivalent to the tangential electric field intensities on the top and bottom BI surfaces of the unit cell) is given by I J Maux = E [M + MY y] 8(x-xo-i AX) 6(y-yo -jAy) i=1 j=l N I J = Mn 1E [A n+ Ajn y] ((x-xo-i x iA) (6(y -y -j Ay), (3) n=l i=lj=l where (xo, Yo) is the lower left corner of the uniform grid and Ax, Ay are the sample distances in the x and y directions, respectively. Mij, My. are the expansion coefficients on the uniform grid, Mn are the expansion coefficients on the original mesh, and Ajn, AYn are the equivalent expansion coefficients on the uniform grid for each single expansion function on the original mesh. Since we work with a MPIE formulation for the BI implementation, we introduce the additional expansion I J V. Ms = E Qij (x-xo-iAx) 6(y-yo- Ay) i=l j=1 N I J Z Qn Aqin (x - xo-iAx) 6(y- o- j y) (4) n=l i=1 j=l for the divergence of the surface currents. To find an appropriate relationship between the uniform expansion (3), (4) and the expansion functions on the original possibly irregular triangular mesh for evaluation of the far interactions, we start with the common MPIE expression [17] Zmn = Z(1) +Z(2) ZmG = -2k2o Jf TGp(rlrs) f m(r) f n(r,) dS, dS Smn Sn 6

+ 2 JJ// Gp(r, r,) V, f m(r) V f n(rs) dSs dS (5) Sm Sn for the calculation of an arbitrary coupling element for the original expansion. Here r, and r are the source and observation points, whereas f and fm are the source and testing basis functions on the triangles with the areas Sn and Sm, respectively. Also, Gp(r, r,) is the scalar periodic Green's function of the array problem (for instance see [1, 3, 2]) and ko is the wave number of free space. 'With qm(r) = V fm(r), qn(rs) = Vs fn(rs), (6) and g(r)= f Gp(r, r,) qn(r) dSs (7) ((2) Ss Zmn can be written as Z() qm(r)(r))dS. (8) S Without loss of generality, the evaluation of the coupling integrals can be performed within the plane z = 0. Based on this observation and introducing the Taylor series expansion 00 g(r)= - aq (x _-xl)ql (y-y1)q2 (9) q=ql+q2=1 around the expansion point r = (x1, Y1, 0) within the plane of integration, Z(2 can be written in the form oo Z( = aq ]Jqm(r)(x - )q (y - y)q2 dS, (10) mn (10) q=ql+q2=1 S which is equivalent to summing up the moments of the testing function qm multiplied by the Taylor coefficients aq. In principle, the expansion point (xl, yl) can be selected arbitrarily but it is numerically adlvantageous to put it at the center of the edge for which qm is defined. 7

The Taylor series is converging if the Green's function singularity is outside the integration domain which is guranteed for non-overlapping source and test subdomains. Based on (10), Z(2) can exactly be calculated employing the equivalent basis functions on the uniform grid if the moments of the equivalent basis functions are equal to the moments of q,. So, in principle the equivalent amplitudes Atom for the mth bais function on the uniform grid can be obtained by evaluating I J fqm(r)(x -x y)q (y - yl)q2 dS = E Aqm (xo + i Ax - x)ql (o + Ay - yl2, (11) S i=1 j=1 q=q +q2 =,...,oo, where the filter property of the a-functions was utilized. Of course, in a numerical implementation, only a finite number of moments can be enforced to be equal. However, if the distance between source and test subdomains is not too small, the function g(r) is well-behaved and an accurate evaluation of the coupling integrals can be achieved with a small number of the lowest order moments. We restrict ourselves to 3x3 = 9 moments, so that each basis function qm must only be related to 9 basis functions of the uniform grid. That is, (11) is evaluated for indices i = (i,1 - 1),..., (im + 1) and j = (jm - 1),..., (jm + 1), where im and jm are the indices of the uniform grid point which is closest to the center of edge m. The resulting 9x9 linear algebraic system can then be solved for the 9 Aq.m. Note, that the moments of the Rao-Wilton-Glisson basis functions on the irregular mesh can be evaluated analytically as for (im +-1) (jm-+1) Z(2) g(o + i zx, yo + j Ay) (12) i=(im-1) j=(jm-l1) The same procedure as above can be employed for the source integral in (7). For this purpose, a Taylor series expansion of the Green's function Gp is introduced with respect to the source 8

point rS for a fixed observation point r. The obtained approximate expression for g is (k, +l) (l,+1) g(r) E E 2Gp(r|xo+AA k Ax, Yo +l )A, (13) 5(13) k=(kn-l) =(ln-1) where kn In are the indices of the closest uniform grid point with respect to the center of edge n and Aqk are the equivalent amplitudes for the respective basis functions on the uniform grid for basis function n on the original irregular mesh. The combination of (12) and (13) results in (im+l) (im+l 1) (kn+l) (/n+l) Zn) E E E Ak 2 Gp(i Ax, j AykAx, 1 Ay) (14) i=(im-1) j=(jm-L) k=(kn-1) l=(ln-1) in which xO, yo were omitted because of the shift invariance of Gp. In matrix form, we can write [Z]() 2 [A]q[Gp][A]Tq. (15) Expression (15) represents a matrix product of the two sparse A-matrices and the fully populated Green's function matrix [Gp] which is due to the evaluation on a uniform grid of Toeplitz form. In some way, the A-matrices can be interpreted as mapping operators from the irregular mesh onto the uniform grid and vice versa. However, it should be noted that usually for the application of a mapping operator in one direction, the corresponding inverse mapping operater must be applied in the reverse direction. The inverse matrix of a sparse mapping matrix would in general give a fully populated matrix resulting in an inefficient and possibly useless method. Therefore, a key aspect of the AIM concept is that the A-matrices or mapping operators between the two grids are constructed in a way which avoids the inverse mapping matrix. As pointed out in [10], the evaluation of the coupling elements as a sum over the uniform grid can be interpreted as a numerical quadrature with the weights given by the elements of the A-matrices. 9

Z(1) in (5) is evaluated in the same way as z(2) except that the x and y components of the magnetic currents must be considered. Therefore, the resulting representation for Zmn in matrix form is [Z] [Z] M = -2ko ([A]x[Gp][A]T + [A]y[Gp][Ay) + 2 [A]q[Gp][A]. (16) 2.3 AIM Implementation for Periodic Arrays For the implementation of AIM, it is assumed that [Z]tal is sufficiently accurate for [Zfar] but not so for the near-zone elements. Therefore, we decompose the AIM matrices as [otal [Z]near + [Z]far (17) AIM = [ZAIM + L-JAIM ( and when this is combined with (2), we can rewrite the original Z-matrices as [Z]approx = ([Z]near - [Z]near) + [Z]IMo (18) With this representation of [Z] the near-zone elements are evaluated without compromise in accuracy. However, since the great majority of [Z]approx includes the Toeplitz kernel [Gp], the associated matrix-vector products can be performed using only 0(rin) memory and 0((ns logns) CPU time. In the final numerical implementation, a near-zone threshold is defined so that [Z]approx is a sufficiently accurate representation of [Z]. This threshold is mostly affected by the quasi-static singularities of the Green's function and can be reduced if the density of the uniform grid is increased. In the case of an infinite periodic Green's function, we must also account for the singularities associated with the elements which neighbor the boundaries of the unit cell since the latter are adjacent to the next periodic cell (see Fig. 2). Therefore, matrix elements associated with test basis functions whose supports are close to these singularties must also be treated as those in the near-zone of the source element. 10

For the calculation of the matrix-vector products in the iterative solver, [Z]fI is not computed explicitly. After mapping the actual source distributions onto the regular grid through the A-matrices, the pertinent matrix-vector products are performed in the DFT domain using a FFT algorithm for the corresponding transformation (for instance see [12]). After transformation back into the spatial domain, the fields on the original mesh are obtained by reverse mapping between the auxiliary unknowns and the original grid unknowns. As indicated by (18), [Z]tA4aj not only contains the intended far-zone contributions but also near-zone contributions, even though these near-zone contributions are not a good approximation for the original elements. These AIM near-zone contributions are calculated prior to starting the iterative solver and subtracted from the original near-zone coupling elements to obtain the exact near-zone terms after performing the matrix-vector products. For the computation of [Z]nAef% it is essential that the coupling contributions on the uniform grid are collected in the spatial domain before the Toeplitz Green's function is transformed into the DFT domain. 3 Results For validation of the AIM accelerated array analysis FE/BI code, we first considered a simple strip dipole FSS structure presented in [20]. The geometry of the structure is illustrated in Fig. 3 together with resonance curves of the power reflection coefficients for different discretization models. The prismatic FE volume meshes were expanded from triangular surface meshes using 20x20, 32x32, and 44x44 subdivisions, respectively. The density of the uniform AIM grid on the BI surfaces in the top and bottom of the unit cell mesh was chosen so that in average approximately 3.5 uniform grid points were placed per triangle side. The nearzone threshold was set so that the original matrix elements were used within 15 uniform grid 11

samples in the x and y directions around the source element center. That is, the near zone was a square of 30 uniform sample distances. The near-zone threshold was kept constant for all frequencies. Based on the above AIM parameters, the accuracy of the obtained results is equivalent to results obtained with a conventional BI implementation. Keeping in mind that the evaluation of the far interactions in a conventional BI formulation is usually done with a low-level numerical integration (one point sampling), in many cases, AIM might even provide more accurate representations of the far-zone integrals. Increasing the mesh density for the unit cell in Fig. 3 from 20x20 to 32x32 subdivisions, shows a shift of the resonance curves to higher frequencies. However, the 44x44 mesh resonance curve lies exactly on the 32x32 curve. The data are compared to results obtained with a method of moments (MoM) formulation employing a multilayered Green's function. Two MoM curves with 3x8 and 9x18 subdivisons for the current expansion on the strip are depicted in the figure. The higher current sampling on the strip gives a slight shift of the resonance curve to lower frequencies. However, the converged FE/BI and MoM results still have a small frequency shift of about 1%. To illustrate the memory and CPU time savings of the AIM acceleration, we analyzed different mesh configurations of a microstrip dipole array. All meshes consisted of one volume prism along the depth of the array and a metallic backing on the bottom surface. Different meshes with increasing numbers of BI unknowns were generated by grouping several array elements into the discretized unit cell. The largest mesh was a four by four array and had 18208 BI unknowns (24448 volume mesh unknowns). The operation frequency was close to the half-wave resonance of the dipole elements. The uniform AIM grid was constructed as described in the preceding validation problem giving results without compromise in accuracy. In Fig. 4, the number of matrix elements in the system matrix is depicted as a measure of 12

the memory requirement of the algorithm. In the conventional BI formulation, the number of matrix elements increases with complexity O(n2) whereas AIM results in an optimal complexity of 0(ns). The most time consuming portion of our FE/BI approach was the BI fill time which is given in Fig. 5 together with the total solution time. However, for the given problem, the BI fill time was larger than that for most practical problems. Basically, due to the large unit cell having a sidelength of more than two wavelengths, the evaluation of the periodic Green's function needed relatively many terms of the applied series representation. To obtain consistent timing results, we kept the same number of terms even for the smaller problems. From Fig. 5 it can be seen that the AIM acceleration reduces the complexity of the BI fill from about O(nr) to about 0(rin). For larger problems, the CPU time requirements will finally be dominated by the iterative solver. The complexity of the CPU time per iteration for the conventional and AIM accelerated BI is depicted in Fig. 6. Results are given for the BiCG and GMRES solvers. In all cases, the GMRES solver (restarted every 50 iterations) needed approximately half the CPU time required by the BiCG solver. This is due to that the GMRES solver needed only one matrix-vector product per iteration whereas the BiCG solver needed two matrix-vector products. As illustrated in Fig. 6, the AIM acceleration reduces the CPU time complexity per matrix-vector product from O(n2) to O(ns log ns). However, the necessary AIM overhead leads to increased CPU time requirements for very small numbers of unknowns. The total CPU time for the iterative solver is given in Fig. 7. Due to the varying numbers of iterations for the different problems, the curves show a more irregular behavior than the timing curves per iteration. However, the trend of the reduced complexity of the AIM accelerated formulation is still obvious. 13

4 Conclusions In this paper, we demonstrated the application of the AIM acceleration scheme whithin a hybrid FE/BI approach for modeling doubly periodic structures. The AIM concept was reviewed and a intuitve mathematical derivation based on Taylor series expansions of the coupling integrals was given. Also, implementation issues related to the array problem were addressed. For the considered planar BI surfaces, the method leads to low 0(ns) complexity for the memory requirements within the BI portion of the approach and the CPU time complexity of the matrix-vector products in the applied iterative solver is reduced to 0(nr log ns). With properly selected AIM parameters, results without any compromise in accuracy are obtained compared to the conventional BI implementation. Acknowledgements This work was supported in part by the Office of Naval Research (ONR-313) through the Naval Air Warfare Center Aircraft Division (NAWCAD 4.5.5.5) contract with Sanders, A Lockheed Martin Corporation. The first author acknowledges the fellowship support from the "German Academic Exchange Service (DAAD)". 14

References [1] D. T. McGrath and V. P. Pyati, "Phased Array Antenna Analysis with the Hybrid Finite Element Method," IEEE Trans. Antennas Propagat., vol. 42, no. 12, pp. 1625-1630, Dec. 1994. [2] E. W. Lucas and T. W. Fontana, "A 3-D Hybrid Finite Element/Boundary Element Method for the Unified Radiation and Scattering Analysis of General Infinite Periodic Arrays," IEEE Trans. Antennas Propagat., vol. 43, no. 2, pp. 145-153, Feb. 1995. [3] D. T. McGrath and V. P. Pyati, "Periodic Structure Analysis Using a Hybrid Finite Element Method,' Radio Science, vol. 31, no. 5, pp. 1173-1179, Sep/Oct. 1996. [4] N. Engheta, W. D. Murphy, V. Rohklin, and M. S. Vassiliou, "The Fast Multipole Method (FMM) for Electromagnetic Scattering Problems," IEEE Trans. Antennas Propagat., vol. 40, no. 6, pp. 634-641, Jun. 1992. [5] R. Coifman, V. Rohklin, and S. Wandzura, "The Fast Multipole Method for the Wave Equation: A Pedestrian Prescription," IEEE Antennas Propagat. Mag., vol. 35, pp. 7-12, Jun. 1993. [6] N. Lu and J.-M. Jin, "Application of the Fast Multipole Method to Finite-Element Boundary-Integral Solution of Scattering Problems," IEEE Trans. Antennas Propagat., vol. 44, no. 6, pp. 781-786, Jun. 1996. [7] S. S. Bindiganavale and J. L. Volakis, "Comparison of Three FMM Techniques for Solving Hybrid FE-BI Systems," IEEE Antennas Propagat. Mag., vol. 39, no. 4, pp. 47-60, Aug. 97. 15

[8] J. M. Song, C. C. Lu, and W. C. Chew, "Multilevel Fast Multipole Algorithm for Electromagnetic Scattering by Large Complex Objects," IEEE Trans. Antennas Propagat., vol. 45, no. 10, pp. 1488-1493, Oct. 1997. [9] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, "A Fast Integral Equation Solver for Electromagnetic Scattering Problems," IEEE AP-S, pp. 416-419, 1994. [10] E. Bleszynski, M. Bleszynski, and T. Jaroszewicz, "AIM: Adaptive Integral Method Compression Algorithm for Solving Large-Scale Electromagnetic Scattering and Radiation Problems," Radio Science, vol. 31, no. 5, pp. 1225-1251, Sep./Oct. 1996. [11] H. T. Anastassiu, M. Smelyanskiy, S. Bindiganavale, and J. L. Volakis, "Scattering from Relatively Flat Surfaces Using the Adaptive Integral Method," Radio Science, vol. 33, no. 1, pp. 7-16, Jan. —Feb. 1998. [12] J. L. Volakis and K. Barkeshli, "Application of the Conjugate Gradient FFT Method to Radiaton and Scattering," Appl. of It. Meth. to Electromag. and Signal Proc., ed. T. Sarkar, Elsevier Pub. Co., pp. 159-240, 1991. [13] J. R. Phillips and J. K. White, "Efficient Capacitance Computation of 3-D Structures Using Generalized Pre-Corrected FFT Methods", Proc. 3rd Topical Meet. Elect. Perfomance Elect. Packaging, Monterey, CA, pp. 253-256, Nov. 2-4, 1994. [14] S. D. Senturia, N. Aluru, and J. White, "Simulating the Behavior of MEMS Devices: Computational Methods and Needs", IEEE Computational Science & Engineering, pp. 30-43, Jan.-Mar. 1997. [15] C. L. Berman, "Grid-Multipole Calculations", SIAM J. Sci. Comput., vol. 16, no. 5, pp. 1082-1091, Sep. 1995. 16

[16] J. Gong, J. L. Volakis, A. Woo, and H. Wang, "A Hybrid Finite Element-Boundary Integral Method for the Analysis of Cavity-Backed Antennas of Arbitrary Shape," IEEE Trans. Antennas Propagat., vol. 42, no. 9, pp. 1233-1242, Sep. 1994. [17] S. M. Rao, D. R. Wilton, and A. W. Glisson, "Electromagnetic Scattering by Surfaces of Arbitrary Shape," IEEE Trans. Antennas Propagat., vol. 30, no. 3, pp. 409-418, May 1982. [18] Y. Saad, Iterative Methods for Sparse Linear Systems, Boston: PWS Pub. Co., 1996. [19] T. Ozdemir and J. L. Volakis, "Triangular Prisms for Edge-Based Vector Finite Element Analysis of Conformal Antennas," IEEE Trans. Antennas Propagat., vol. 45, no. 5, pp. 788-797, May 1997. [20] R. Mittra, C. H. Chan, and T. Cwik, "Techniques for Analyzing Frequency Selective Surfaces - A Review," Proc. IEEE, vol. 76, no. 12, pp. 1593-1615, Dec. 1988. 17

Figure Captions Figure 1: Unit cell of array Figure 2: Triangular BI mesh with periodic image sources Figure 3: Resonance curves for strip dipole FSS as sugested in [20] Figure 4: Number of matrix elements in the system matrix (FE and BI near-zone) Figure 5: BI fill time and total solution time for test problem of microstrip dipole array Figure 6: CPU time per iteration for test problem of microstrip dipole array Figure 7: Total solver CPU time for test problem of microstrip dipole array 18

BI PBC PBC BI Figure 1

..: ---. r mage ~-.......'. ~ ~.-.-.-.-........... ~.:.:.:.:.:.: ~. -.:. I'. '.''.'.'.'.'.' I''''''''' '''''''''1 I......-....-.-...-. ~........ sources, --->. /.'.'.'..' ""A!::' ~ ~:':': ':': I <..-......-.....-. ~.:.:.:.:.:.:.:.:.:.:/ ~:.::.::.::.::.: ';,:i.i.i.i... i I'", -.-;,:::::::::;,:, primary source Figure 2

Power Reflection Coefficient 0 0 0 0 0 0b \ b o b a S a * x o -,*\ \\, * * * g W 00 4t.)O (1 -tU00 -^ ^ 0.2 cm 1cm 00 " F " fl0.5cm N

- W 106 x *Ea *4-m o 105 105.......................I4.-.............:::nti: A:l::::........... *................................................... i 4..... 4.... 4 b-' —''....<..................................................................................................................,................ 5.;X..... 4............................................... C.... C.... I...4........................ 4 —. —. ---.-.-v..... f —. ---.......;.......................:..................................................................................... 103 104 Number of BI unknowns Figure 4

104 ci W W 0 0 pog W 5 II'= 44 RI fill time.. total time........ onven'1 "U"BI...... -AIM aiccelerated ElI 103 103 104 Number of RI unknowns Figure 5

103 102 101 I_ _ _ I I _ _......I.....4.....I.............................................9.........I.............4....9...i.....4........................................................................................................................ c.....................................j.....................................................:::::: ' i ' ': e........................................................................................................................................................................::: i i i:: ':::: ' ~::::?:,,::: '::;:::: - - i................... —............................... ---- -i. r ~ohv.efitldnal BI. | g -:; $ 4 -...........................................~...............................................................:...................... ". 4.. 4... 'l"'.. S -..-...:..................................................................................:.................... - -:................................................................................ i................. *'...............................................................................................................................................................................................i...''..................... i..................................................................i......................:::::::::::::::1::::::::::::: 103 104 Number of BI unknowns Figure 6

101 *c v) V 0 0 PM CD Number of BI unknowns Figure 7