035067-2-T ADAPTIVE INTEGRAL METHOD FOR HYBRID FE-BI MODELING OF 3D DOUBLY PERIODIC STRUCTURES T. F. Eibert J. L. Volakis Rockwell International Science Center 1049 Camino Dos Rios P.O. Box 1085 Thousand Oaks, CA 91360 July 1998 35067-2-T = RL-2483

PROJECT INFORMATION PROJECT TITLE: REPORT TITLE: U-M REPORT No.: CONTRACT START DATE: END DATE: DATE: SPONSOR: Hybrid Finite Element-AIM Methods for Antennas Adaptive Integral Method for Hybrid FE-BI Modeling of 3D Doubly Periodic Structures 035067-12 —T October 1998 March 1999 July 1998 Marek Bleszynski Rockwell International Science Center 1049 Camino Dos Rios P.O. Box 1085 Thousand Oaks, CA 91360 SPONSOR CONTRACT No.: U-M PRINCIPAL INVESTIGATOR: John L. Volakis EECS Dept. University of Michigan 1301 Beal Ave Ann Arbor, MI48109-2122 Phone: (313) 764-0500 FAX: (313) 747-2106 volakis @umich.edu http://www-personal.engin.umich.edu/-volakis/ CONTRIBUTORS TO THIS REPORT: T. F. Eibert and J. L. 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, AIM is derived using Taylor series expansions for the involved coupling integrals. Implementation issues are discussed as related to the periodic formulation and central processing unit (CPU) timings / 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 a unit cell representing the periodic structure, 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 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. Until 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 [7]. However, for planar BI terminations, AIM is the method of choice because it directly results in low 0(ns) storage and 0(9ns log ns) CPU time requirements for executing 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 log ns) complexity. Also, FMM cannot be directly applied to periodic problems because it is based on a multipole expansion of 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 use of Fast Fourier Transform (FFT) algorithms has been extensively applied in Conjugate Gradient (CG) - FFT like algorithms for many years [12, 13]. Originally, these methods were limited due to their inherent requirement for uniform meshes. To obtain more flexibility, 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 [14, 15], algorithms referred to as "Precorrected-FFT Methods" were presented and applied to electrostatics. The projections between the uniform and irregular meshes were achieved via multipole expansions of the potentials as described in [16] and 6-functions were used 1

to expand the charges on the uniform grid. In [17], a full-wave method was published for the analysis of cavity backed antennas. Both the uniform and the irregular BI meshes were based on triangular surface elements 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 (-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 independent 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 CPU time and memory requirements of the BI portion of a hybrid FE/BI approach pertaining to 3D doubly periodic structures. The paper focuses especially on issues related to the implementation of 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 provides CPU speed-up and memory reduction curves afforded by the AIM implementation. 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 Acros A oss E 0 0 0 int fint A: A~ 6 \ Ey \^ 0 Ztop 0 0 0 Eint i ={ \n 1 Across 0 Abozd Ebound u0 Z bonp Ebound fbound botnAc 0 A bound E dbound A2,bot A bot I Ebot ) L Z bot bottom where only a unit cell of the array (see Fig. 1) is discretized. The A-matrices are sparse (20 to 40 non-zero elements per row) 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 employ triangular prismatic elements to generate the volume mesh giving rise to triangular surface meshes with Rao-Wilton-Glisson basis functions [18] for the magnetic currents on the planar BI surfaces. For arrays with arbitrary scan angles, all matrices in (1) may be nonsymmetric due to the PBCs and the periodic Green's function. For the solution of the system, 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) [19]. Due to the non-symmetric 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 2

to the storage of the search vectors within the restart loop. Also, additional computational effort is required for the orthogonalization of the search vectors [19]. The computational complexity per matrix-vector product of the total A-matrix is of O(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 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 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 method of moments (MoM) coupling integrals. Therefore, this subsection reviews 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 decompose the Z-matrices as [Z] = [z]near + [Z]far, (2) where [Z]near contains the elements of [Z] which are "near" the self-cell. Correspondingly, [Z]far 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 [Z]far [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 [13]. 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 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 overlays 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-1 J-1 Maux = E E [M + My ] 6(x-xo-i Ax) 6(y-yo-jAy) i=0 j=O N I-1J-1 -= E[An + An ] b( _x-0o-i Ax) (y- yo-j Ay), (3) n=l i=0 j=O 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. M', M.. are the expansion coefficients on the uniform grid, Mn are the expansion coefficients on the original mesh, and Ain, Ayj are the equivalent expansion coefficients on the uniform grid for each single expansion function on the original mesh. That is, summation over the expansion coefficients Mn of the original mesh multiplied by the equivalent expansion coefficients Ai n, Aijn gives the expansion coefficients MX, My. of the total surface current densities on the uniform mesh. Since we work with a MPIE 3

formulation for the BI implementation, we introduce the additional expansion I-1 J —1 Vs Ma = 6E Qi S( - o-iAx) 6(y - o - j y) i=O j=:O N I-1J-1 Qn AQn 6( - x - i Ax) 6(y - Y - j y) (4) n=l i=O j=O for the surface divergence of the surface currents. To find an appropriate relationship between the uniform expansion (3), (4) and the basis functions on the original possibly irregular triangular mesh for evaluation of the far interactions, we start with the common MPIE expression [18] Zmn = ) + (2) =- 2k 2 /JJGp(rrs) fm(r) fr(r) dS dS Sm S~ + 2 / (r l rs) V f m (r) Vs 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 fn and fm are the source and testing basis functions on the triangles with the areas Sn and Sm, respectively. Also, Gp(rlrs) is the scalar periodic Green's function of the array problem (for instance see [1, 2, 3]), ko is the free space wave number, and Vs' denotes the surface divergence. With the definitions qm(r) = Vs f(r), qn(rs) = Vs fn(rs) (6) and g(r)= JGp(rlrs) qn(rs) dSs, (7) Sn Zmn can be written as Z( = 2 Jqm(r)g(r) dS. (8) Sm Without loss of generality, the evaluation of the coupling integrals can be performed in the plane z = 0. Based on this observation and introducing the Taylor series expansion 00 g(r) = = aq (x -- xl)q (y - yl)q2 (9) q=ql+q2=0 around the expansion point r = (xl, Yi) in the plane z = 0, Z(2) can be written in the form 00oo Z(2) =2 aq qrnm(r)(x -l)q (y - yl)qd, (10) q=ql+q2-= Sn which is equivalent to summing up the moments of the testing function qm multiplied by the Taylor coefficients aq. In principle, the expansion point (xi,yl) can be selected arbitrarily 4

but it is numerically advantageous to put it at the center of the edge for which qm is defined [10]. The Taylor series converges if the Green's function singularity is outside the integration domain which is guaranteed for non-overlapping source and test subdomains. Based on (10), Z(2) can be calculated exactly by employing the equivalent basis functions on the uniform grid if the moments of the equivalent basis functions are equal to the moments of qm. So, in principle the equivalent amplitudes A9Qm for the mth basis function on the uniform grid can be obtained by enforcing the equality I-1 J-1 //qm(r)(x- Xl)ql (y - yl)q2 dS - AQ'm (xo + i Ax - Xl)ql (Yo +j Ay - yl)q2, (11) Sm i=O j=O q=ql + q2 =0,...,oo, where the filter property of the d-functions was utilized. 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 3 x 3 = 9 moments, so that each basis function qm must only be related to nine basis functions of the uniform grid. That is, (11) is evaluated for indices i = (im - 1),..., (irn + 1) and j = (jm - 1)...(j Um + 1), where im and jm are the indices of the uniform grid point which is closest to the center of edge m. The resulting 9 x 9 linear algebraic system can then be solved for the nine A9'm coefficients. Note that the moments of the Rao-Wilton-Glisson basis functions on the irregular mesh can be evaluated analytically as for instance discussed in [11]. The approximate expression for (m2) can be finally written as (im+l) (jm+'1) Z(2) 2 g(xo + iAx,yo + jAy) AQm, (12) i-(im-l) j=(jm-l) where again the filter property of the d-functions in (4) was utilized. A similar procedure 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 point r, for a fixed observation point r. The obtained approximate expression for g is (kn+l) (n+ 1) g(r) GP(rI(xo +kAx,yo+ lAy)) A (13) k=(kn-1) l=(ln-1) where kn, In are the indices of the closest uniform grid point with respect to the center of edge n and AQ n 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) (jm+l 1) (kn+l) (/n+l) Z(2) 2 E AQ n Gp((i Ax, j Ay)I(k Ax, lAy)) AQ'm (14) i-=(im-l) j-=(jm-l) k=(kn-1) I=(ln-l) in which xO, yo were omnitted because of the shift invariance of Gp. In matrix form, we can write [Z](2) 2 [A]Q[Gp][A]. (15) 5

Expression (15) represents a matrix product of the two sparse A-matrices and the fully populated Green's function matrix [Gp] of Toeplitz form. In some way, the A-matrices can be interpreted as mapping operators from the irregular mesh onto the uniform grid and back. Usually reverse mapping for a given operator is obtained by taking the inverse of the operator. 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, in the case of AIM it is important to note that the A-matrices or "mapping operators" between the two grids are constructed in a way which avoids the generation of inverse matrix. Z in (5) is evaluated in a similar way as (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]AIM -2ko ([A] [Gp][A] + [A]y[Gp][A]) + 2 [A]Q[Gp][A], (16) where [A]x and [A]y are sparse matrices containing the equivalent current amplitudes introduced in (3). 2.3 AIM Implementation for Periodic Arrays For the implementation of AIM, it is assumed that [Z]AIM is sufficiently accurate for [Z]far but not so for the near-zone elements. Therefore, we decompose the AIM matrices as [Z]AIM = [Z]AIM + [Z]AiM (17) and when [Z]far in (2) is approximated by [Z]QalM, we can rewrite the original Z-matrices as [Z]approx = ([Z]near - [Z] ) [Z]AIM (18) With this approximate 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(ns) memory and 0(ns logns) CPU time. In the final numerical implementation, a near-zone threshold is defined so that [Z]a;prox 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 periodic structures, the test functions in the near-zone around the source function must be extended by test functions whose supports are close to an image source in one of the neighboring periodic cells (see Fig. 2). For the calculation of the matrix-vector products in the iterative solver, [Z]AfaM is not computed explicitly. After mapping the actual source distributions onto the uniform grid through the A-matrices, the pertinent matrix-vector products are performed in the DFT domain using a FFT algorithm for the corresponding transformation (see for instance [13]). 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]AIM not only contains the intended far —zone contributions but also near-zone contributions, even though these near-zone contributions are not a good approximation to the original near-zone 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]IAea, it is essential that the coupling contributions on the uniform grid are collected 6

in the spatial domain (performing the convolutions) before the Toeplitz Green's function is transformed to the DFT domain. In the DFT domain, for each basis function the whole FFT procedure must be performed to obtain the elements of one column of [Z]Xf. Therfore, computation of [Z]"/' in the DFT domain would require O(n2 log ns) CPU time. 3 Results For validation of the AIM accelerated array analysis FE/BI code, we first consider a simple strip dipole FSS structure as presented in [21]. 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 20 x 20 32 x 32, and 44 x 44 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 on average approximately 3.5 uniform grid points were placed per triangle side. The near-zone threshold was set so that the original matrix elements were used within 15 uniform grid samples in the x and y directions around the source element center. That is, the near zone was a square of 30 x 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 FE/BI implementation. Keeping in mind that the evaluation of the far interactions in a conventional FE/BI formulation is usually done with a low-level numerical integration (one point sampling), AIM might even provide more accurate representations of the far-zone integrals. Increasing the mesh density for the unit cell in Fig. 3 from 20 x 20 to 32 x 32 subdivisions, shows a shift of the resonance curves to higher frequencies. However, the 44 x 44 mesh resonance curve coincides with the 32 x 32 curve. The data are compared to results obtained with a MoM formulation employing a multilayered Green's function. Two MoM curves with 3 x 8 and 9 x 18 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/AIM 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 t)y 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 dependent on the number of BI unknowns is depicted as a measure of the memory requirements of the algorithm (12 bytes per matrix element for single precision and sparse matrix storage). In the conventional BI formulation, the number of matrix elements increases with complexity O(n2) whereas AIM results in an optimal complexity of (0(rns). 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. 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 7

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(n2) to about 0(ns) (constant number of terms in the series representation of the Green's function assumed for both values). The complexity of the CPU time per iteration for the conventional and AIM accelerated BI approaches is depicted in Fig. 6. Results are given for the BiCG and GMRES solvers. In all cases, the GMRES solver (restarted every 50 iterations) needs approximately half the CPU time required by the BiCG solver since the GMRES solver needs only one matrix-vector product per iteration whereas the BiCG solver needs two matrix-vector products. As illustrated in Fig. 6, the AIM acceleration reduces the CPU time complexity per matrix-vector product from O(nr2) to O(ns logns). However, the necessary AIM overhead leads to increased CPU time requirements for very small numbers of unknowns. The total CPU time for the iterative solvers 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. 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 an intuitive mathematical derivation based on Taylor series expansions of the coupling integrals was given. Also, implementation issues relating to the array problem were addressed. For the considered planar BI surfaces, the method leads to low 0(rns) 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(n log n). With properly selected AIM parameters, results without any compromise in accuracy are obtained compared to the conventional BI implementation. Acknowledgements The first author acknowledges the fellowship support from the "German Academic Exchange Service (DAAD)". 8

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. [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] N. Bojarski, k-Space Formulation of the Electromagnetic Scattering Problem, Techn. Report, AFAL-TR-71-5, Mar. 1971. [13] 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. 9

[14] 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. [15] 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. [16] C. L. Berman, "Grid-Multipole Calculations", SIAM J. Sci. Comput., vol. 16, no. 5, pp. 1082-1091, Sep. 1995. [17] 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. [18] 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. [19] Y. Saad, Iterative Methods for Sparse Linear Systems, Boston: PWS Pub. Co., 1996. [20] 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. [21] 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. 10

Figure Captions Figure 1: Unit cell of array Figure 2: Triangular 13I mesh with periodic image sources Figure 3: Resonance curves for strip dipole FSS as sugested in [21] 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 11

BI PBC PBC BI Figure 1: - - image sources. - - - \ \ I I\ II \ \ I / / \ \ I / / / primary source Figure 2:

Number of matrix elements 0~ Power Reflection Coefficient C) C: ON:,MI CZ.e rw 5 cr m *1 0 tz 0* 0 r 0 0 4 rA C'D 0 N 1CQ 0 -

104 W ~ IIII/I a~l AIM acceer t 103 103 Number of RI unknowns Figure 5: 101 E100 10-1 103 Number of RI unknowns Figure 6: 104 104 14

103 o." 102 101 Io, *T 00, La A*qr IL RT MI acamwe BiCG --- GMRES M-1 MI y tirelliallA 103 104 Number of RI unknowns Figure 7: 15