Method of Moments Based on Equivalent Periodic Problem and FFT with NURBS Surfaces for Analysis of Multilayer Periodic Structures

: In this paper, an efficient technique of computation of method of moments (MM) matrix entries for multilayer periodic structures with NURBS surface and Bézier patches modelling is proposed. An approximation in terms of constant pulses of generalized rooftop basis functions (BFs) defined on Bézier patches is proposed. This approximation leads discrete convolutions instead of usual continuous convolution between Green’s functions and BFs obtained by the direct mixed potential integral equation (MPIE) approach. An equivalent periodic problem (EPP) which contains the original problem is proposed to transform the discrete convolutions in discrete cyclic convolutions. The resultant discrete cyclic convolutions are computed by efficiently using the Fast Fourier Transform (FFT) procedure. The performance of the proposed method and direct computation of the MM entries are compared for phases of reflection coefficient. The proposed method is between 9 and 50 times faster than the direct computation for phase errors less than 1 deg. The proposed method exhibits a behaviour of CPU time consumption of O(N b Log10N b ) as the number N b of BFs increases. This behaviour provides significant CPU time savings with respect to the expected behaviour of O(N b2 ) provided by the direct computation of the MM matrix entries. This efficient pre ‐ processing computation show consumption of CPU time of 2.247 s for 2  256  256 integrals. Using these pre ‐ processing results, an analysis of a unit cell was carried out. Comparisons of consumption of CPU time by the direct MM and the proposed method was carried out and the comparisons show that the proposed method can be between 10 and 50 times faster than the direct method when the number of BFs increases from 138 to 772. Dependence of CPU time consumption with as a function of BFs was carried out. O (N b Log 10 (N b )) behaviour of CPU time consumption by the proposed method was found. This behaviour provides significant CPU time savings with respect to the expected behaviour of O(N b2 ) provided by the direct computation of the MM matrix entries. Finally, several validations of the proposed and direct methods were successfully carried out.


Introduction
In the design and analysis of electromagnetic devices likes frequency selective surfaces (FSSs) [1], reflectarrays/transmitarrays [2], leaky wave antennas [3] and metasurfaces antennas [4], efficient electromagnetic analysis tools of multilayer periodic structures are required. Although reflectarrays/transmitarrays, leaky wave antennas and metasurfaces antennas are not strictly periodic structures, in the design of these antennas, it is common practice to assume that each element of the antenna is located in periodic environment. This is known as the local periodicity assumption [2]. Some popular of numerical tools of periodic structures are based on Finite Element (FE) [5], Finite Difference Time Domain (FDTD) [6] and MM [7], all methods under periodic boundary conditions. Although FE and FDTD are mature methods, these methods involve volumetric meshes to model the multilayer medium which hosts the unit cell of the periodic structure. This volumetric mesh involves a large amount of numerical computations which slow down the electromagnetic analysis. Thus, the MM is preferable because the multilayer medium is modelled by Green's functions and the volumetric mesh is avoided. The work shown in this paper will be focused in MM formulation for an analysis of multilayer periodic structures. There are recent papers that attempted sophisticated improvements on the MM [8][9][10], also for multilayer periodic structures [11][12][13]. The MM is usually used to solve electric field integral equations (EFIEs) where the unknowns are the induced current densities on the metallic layout. These EFIEs show hypersingular behaviour of the kernel (Green's functions) which cause difficulties in the solution procedure [14]. This hypersingular behaviour may be avoided if the electric fields are expressed in terms of vector and scalar potentials with weakly singular kernels. This led to the development of MPIEs [15]. In this latter approach, for the case of periodic structures, one has to face with the computation of multilayer periodic Green's functions consisting of slowly convergent double infinite summations. Fortunately, there are preprocessing techniques to interpolate efficiently the Green's functions in terms of Chebyshev polynomials after extracting known singular behaviour around the source point [11]. Therefore, this approach involves surface mesh to model the geometry layout hosted in the unit cell of the multilayer periodic structures with accurate computation of the multilayer periodic Green's functions. Once an accurate computation of multilayer periodic Green's functions is available, the induced current density on the metallic layout is expanded in terms on BFs weighted by unknown coefficients. When this expansion is taken into account in the MPIEs and a method of weighted residual is applied, a linear system of equations is obtained. The elements of the coefficient matrix of the resultant linear system of equations involve the computation of convolutions between the Green's function and the BFs used in the expansion of the density current. These convolutions involve the integration of singular behaviour of the integrand introduced by the Green's function when the observation point and the source point are close. Fortunately, there are direct methods of numerical integration of the singular behaviour of the integrand as double exponential (DE) formulas [16] or Ma-Rokhlin-Wandzura (MRW) quadrature rules [17,18].
On the other hand, in order to provide accurate electromagnetic models of complex geometries of the layout, NURBS surfaces are usually used by computer-aided geometry design (CAGD) tools [19]. NURBS surfaces are quite useful because they provide invariance under rotation, scaling, translation and perspective transformation of control points. Moreover, they allow complex shapes to be defined by means of a small number of NURBS. These NURBS surfaces are efficiently written in terms of piecewise Bézier patches as it is described in [19] using Cox-de Boor transformation algorithm [20]. Bézier patches are parametric surfaces defined by Bernstein polynomials [21,22] and they are suitable for numerical computation of parameters associated with the surface (curvature, derivatives, integrations) thanks to the properties of Bernstein polynomials. Thus, Bézier patches are suitable domains to define known subsectional BFs which approximates the induced density currents on the metallic surface. In fact, generalized subsectional rooftop BFs are defined on a pair of adjacent Bézier patches in [19] to approximate the surface density currents induced on metallic surfaces. Despite all these improvements, the direct computation of MM matrix entries leads to computational complexity of CPU time consumption as a function of the number of BFs Nb which is roughly O(Nb 2 ) [23,24]. This computational complexity is inherent to the direct computation of MM matrix entries because it requires the computation of a MM matrix with dimensions of NbNb element by element.
In [25] and [26], the Conjugate Gradient Fast Fourier Transform (CG-FFT) method is proposed to solve electromagnetic problems. This CG-FFT method uses as BFs rooftops defined on a regular rectangular mesh to expand the current density. The advantage of this approach is that the convolutions between Green's functions and the BFs are expressed as discrete cyclic convolutions which are computed by FFT algorithm. In this way, the consumption of CPU time involved in the computation of these cyclic convolutions are proportional to NbLog10(Nb). Thus, the behaviour of CPU time with respect to the number of BFs, Nb, involves important CPU time savings as Nb increases when it is compared to that provided by the direct computation of MM matrix entries. However, this approach has a drawback because it uses regular rectangular meshes. In order to model accurately complex geometries (e.g., geometries with a combination of thin strips and/or patches, patches with curve boundaries, etc.) it needs a very dense mesh. This requires a high number of rooftops to approximate the surface currents which leads to systems of equations with a high number of unknowns.
In other very general formulations for multilayered problems [27,28], Rao-Wilton-Glisson (RWG) BFs [29] defined on triangular surfaces are also approximated in an expansion of pulses. However, this approximation is accurate for sufficient separation of the source and the observation triangular surfaces. Under this condition, the MM matrix entries are computed using a conventional pre-corrected FFT algorithm [30].
In this work, we used generalized rooftop defined over pair of adjacent Bézier patches. These rooftops are approximated in an expansion of pulses with enough dense equi-spaced meshes in order to keep the high-order description of the geometry. Thus, the approach is suitable for complex geometries (geometries with a combination of thin strips and/or patches, patches with curve boundaries, etc.) when very dense equi-spaced meshes are used in the pulse expansion. This approximation leads to discrete convolutions instead of usual continuous convolution between Green's functions and BFs. Since the Green's functions and generalized rooftops involved in the discrete convolutions are not strictly periodic functions, an EPP which contains the original problem is proposed [26]. We show that this new approach leads to discrete cyclic convolutions which were be efficiently computed for a very dense equi-spaced mesh by means of an FFT procedure. Since the original problem is contained in the EPP, this approach does not require considering different formulations when the source and the observation points are close or far. Moreover, unlike the CG-FFT method, the increment of density of the equi-spaced mesh does not lead to an increment of number of unknowns.
This paper is organized as follows. Section II shows a description of the problem which is faceted by the direct MM using a generalized rooftop defined on a pair of adjacent Bézier patches and 'razorblade' as weighted functions. In this section, the formulation of direct numerical integration involved in the MM matrix entries is shown. Section III shows an efficient computation of these matrix entries using the pulse approximation of the generalized rooftop BFs and EPP approach. Section IV shows an efficient computation of integrals of the Green's functions in each pulse domain which are required in the formulation described in the previous section. Section V shows convergence results and CPU time consumption obtained by the proposed method and by the direct computation of the MM matrix entries. Validations of both numerical techniques are also shown with different periodic structures. Finally, the conclusions are provided in Section VI. Figure 1 shows a multilayer periodic structure. The patches are assumed to be PEC with negligible thickness. The multilayer substrate consists of N lossy dielectric layers with thickness dk and complex permittivity k=0r,k(1-jtank) (k=1,…,N). The lower limit of the multilayer substrate is a ground plane. The qth-interface hosts a periodical array of patches with arbitrary geometry. In order to simplify the notation, we show a formulation with unique interface which hosts a periodical array of patches. The formulation can be easily extended with more interfaces with periodical structures. The periodic structure of Figure 1 is illuminated by a linearly polarized plane wave with an arbitrary polarization direction and the incidence direction is given by the angular spherical coordinate inc and inc. A time dependence of the type e jt is assumed and this dependence is suppressed throughout. In order to determine the electric field scattered by the periodic structure of Figure 1, we need to determine the induced surface current densities J(x,y) hosted on the unit cell of the qth-interface of the periodic array of patches. These surface current densities can be obtained by solving MPIE [15]:

Description of the Problem
where Et exc (x,y,z=-hq) is the tangential electric field generated in the observation point (x,y,z=-hq) by the plane wave impinging on the multilayer substrate in the absence of the patches. G A xx and G  are the periodic Green's functions for the x-component of the vector potential and the scalar potential, respectively, of the multilayer substrate of Figure 1.  is the surface charge density induced on the surface S of the patch hosted in the unit cell of the periodic array of patches of the qth-interface. The induced surface charge density  is related with the induced surface current density J by the known continuity equation: Note that gradient operators work on prime variables x' and y' of the source point on the surfaces of patches. We would like to point out that in this work, the periodic Green's function of the multilayer substrate for the x-component of the vector potential G A xx and for the scalar potential G  were efficiently obtained by the pre-processing procedure of interpolation described in [11]. In this procedure, the periodic Green's functions of the multilayer substrate are judiciously interpolated in the spatial domain in terms of 2-D Chebyshev polynomials after extracting the singular behaviour of the Green's functions around the source points (which includes the source singularities plus the images through the closest layers). If the MM is used to solve the MPIE shown in Equation (1), J(x,y) would have to be expanded in terms of known BFs Jj(x,y) (j=1,…,Nb), as shown below: Since the surface charge densities and current densities are related by the continuity Equation (2), the BF for the surface charge densities j(x,y) are related with the BF for the current densities Jj(x,y) in similar way as shown in Equation (3). In this paper, the surface of the metallic layout hosted in the unit cell of the qth-interface is modelled by NURBS surfaces. These NURBS surfaces are efficiently written in terms of piecewise Bézier patches as is described in [19] using the Cox-de Boor transformation algorithm [20]. Thus, the BFs Jj(x,y) (j=1,…,Nb) used in this paper are generalized "rooftop" functions defined on pairs of adjacent Bézier patches that share a common boundary line, as described in [19]. According to [19], the BF j(x,y) is constant in each Bézier patch.
In order to solve the MPIE given in Equation (1), the method of weighted residual was carried out using "razor-blade" as the weighting function. These "razor-blades" are defined in the qthinterface over the isoparametric curved lines Ci. This curved line joins the centres of the pair of adjacent Bézier patches associated to each BF Jj(x,y) (see [19]). Thus, when Equation (3) is introduced into Equation (1) and "razor-blade" functions are used, the following system of linear equations for the unknown coefficients cj is obtained: where the coefficients ei of the system of the linear equations can be computed by the next line integral: The elements of the coefficient matrix of the system of linear equations of (4) can be broken down in inductive and capacitive contribution: where the inductive contribution of each element of the coefficient matrix is given in terms of the line integral of the bi-dimensional convolutions between periodic Green's function G A xx for the xcomponent of the vector potential and the BF Jj(x,y), as shown in Equation (7). These bi-dimensional convolutions involve the integration of singular behaviour of the integrand introduced by the periodic Green's function when the observation point (x,y,z=-hq) and source point (x',y',z'=-hq) are close.
The numerical integration of the singular behaviour of the integrand can be accurately computed by DE formulas [16] or MRW quadrature rules [17,18]. On the other hand, the capacitive contribution is given by the line integral of the gradient of the bi-dimensional convolution between periodic Green's function G  for scalar potential and the BF j(x,y). Since the capacitive contribution involves a line integral whose integrand is an exact differential, this line integral can be analytically expressed as the difference of bi-dimensional convolutions, as shown in Equation (8). According to [19], the observation points (x + i, y + i, -hq) and (xi, yi, -hq) are the centres of the pair of adjacent Bézier patches associated to BF j(x,y) (i.e., the extremes points of the isoparametric curved lines Ci).
Again, the bi-dimensional convolutions involve the integration of singular behaviour which can be accurately computed by DE formulas [16] or MRW quadrature rules [17,18].
The main drawback of the direct computation of the inductive and capacitive contributions of the coefficient matrix given in Equations (7) and (8) is that the CPU time consumption is proportional to the size of the coefficient matrix of the linear system of equations (i.e., the computational complexity of CPU time consumption as a function of the number of unknown Nb is roughly O(Nb 2 )) [23,24]. Thus, if a very high number Nb of the BFs is required for the approximation of the induced surface density current on the patches, the required CPU time consumption involved in the computation of all elements of coefficient matrix may be prohibitive. In the next section, we describe a procedure for the computation of the inductive and capacitive contributions which involves the FFT algorithm. This approach does not require considering the different formulations when the source and the observation points are close or far, as is required in [27,28]. This procedure provides significant CPU time saving and behaviour of CPU time consumption of O(NbLog10Nb) as Nb increases, as will be shown in the following sections.

Efficient Computation of MM Matrix Entries Using Pulses Expansions of BFs
In the previous section, the induced surface current J(x,y) and charge (x,y) density are expanded in terms of known BFs Jj(x,y) and j(x,y) (j=1,…,Nb) to solve the MPIE given in Equation (1). These BFs can be expanded in terms of pulses as shown below: where xm'=x0+m'x, yn'=y0+n'y are equi-spaced mesh points and P(•) is the pulse function defined as The expansions given in Equations (9) and (10) approximate the values of the BF Jj(x',y') and j(x',y') as the constant values Jj(xm',yn') and j(xm',yn'), respectively, in the intervals xm'-x/2<x'<xm'+x/2 and yn'-y/2<y'<yn'+y/2. Thus, the pulse approximation of the BF improves as the density of the equispaced mesh increases (i.e., the values of Nx and Ny of the upper limit of the expansion increase as the values of x and y decrease). We would like to point out that the BFs are defined on pair of adjacent Bézier patches associated with NURBS surfaces. Thus, the BFs are capable of taking into account a high-order description of the geometry of layout hosted by the unit cell. Since these same BFs are expanded in terms of pulses with very dense equi-spaced meshes, the high-order description of the geometry is kept.
When Equations (9) and (10) are introduced into Equations (7) and (8), the following inductive and capacitive contributions of the coefficients matrix are obtained: and The direct evaluation of the inductive and capacitive contributions given in Equations (12) and (13) involves the summation of double series Equations (14) and (15) whose addends involve double integrals given by Equations (16) and (17). Moreover, these double series have a high number of addends since a high density of the equi-spaced mesh is required in order to provide accurate approximation of the BFs in terms of pulses. Thus, this approach may seem inappropriate. Let us define the following discrete functions for the BF of the surface current and charge densities: In Equations (18) and (19) the m'n'th-element of matrixes of size (2Nx+1)x(2Ny+1) is defined as the values of the BFs in the equi-spaced mesh point (xm',yn'). Now, let us define the following discrete functions for the double integrals given in Equations (16) and (17) (14) and (15), expressions in terms of discrete convolution are obtained:   Second, we multiply, element by element, the elements of the resultant discrete functions computed in the previous step.  Finally, we compute the inverse 2D-FFT of the resultant discrete function provided by the multiplications, element by element, in the previous step. However, the discrete functions J d j,x/y and Gq d,A involved in Equation (22) and the discrete functions  d j and Gq d, involved in Equation (23) are not discrete periodic functions. Thus, in this case, an EPP has to be found [25,26]. In our case, we define the EPP as the discrete functions J d j,x/y[m',n'] and  d j[m',n'] are defined in Equations (18) and (19) (20) and (21) in the whole domain of the discrete variables m and n. Thus, the discrete functions involved in Equations (22) and (23) are periodic with periods 2Nx+1 and 2Ny+1. In this way, Equation (22) and (23) are discrete cyclic convolutions and can be efficiently evaluated by the previous 2D-FFTs procedure. Note that the EPP has a periodic cell with dimensions of 2a2b which are twice the original periodic cell. This period is the minimum period which guarantees that aliasing is avoided [26]. On the order hand, J d j,x/y[m',n'] and Gq d,A [m,n] (or  d j[m',n'] and Gq d, [m,n]) contains values of the original problem for 0<m'<Nx, 0<n'<Ny. This fact guarantees that the convolutions of the EPP contains the convolutions of the original problem for 0<m'<Nx, 0<n'<Ny [26]. We would like to point out that thanks to EPP, this approach does not require considering different formulations when the source and the observation points are close or far, as is required in [27,28].
Since the FFT algorithm is very efficient, the main computational cost of this procedure is the computation of the discrete functions defined in Equations (20) and (21). Since Green's functions G A xx and G  involved in Equation (16), (17), (20) and (21) are efficiently interpolated in terms of 2D-Chebyshev polynomial [11], the computational cost of the computation of Gq d,A [m,n] and Gq d, [m,n] is improved with respect to that obtained if the conventional computation of multilayer periodic Green's functions was made by means of slowly convergent double infinite summations. Despite this improvement, in this work, we make an additional effort to alleviate this computational cost in the next section.
Finally, the computation of the inductive contributions given in Equation (12) of the elements of the coefficient matrix is computed by conventional quadrature (for example Gauss-Legendre quadrature rules [31]) where the samples of the components of the vector function fj ind (x,y,z=-hq) are computed from the values of the discrete functions f d,ind j,q,x/y[m,n] by conventional bilinear interpolation [31]. A similar procedure is carried out for the computation of fj cap (x,y,z=-hq) from the samples f d,cap j,q[m,n]. Since a dense equi-spaced mesh is required to approximate efficiently the BFs by pulses expansions given in Equations (9) and (10), the bilinear interpolation provides enough accuracy for practical purposes. Thus, once discrete functions f d,ind j,q,x/y[m,n] and f d,cap j,q[m,n] are available for jth-BF, these discrete functions are reused for the computation Zij ind and Zij cap given in (12) and (13) for i=1,…,Nb. Thus, CPU time savings with respect to the direct method (Equations (7) and (8)) are expected as the number of BFs increases.

Efficient Computation of Gq d,A and Gq d,
In this section, we show an efficient procedure to compute efficiently the discrete functions Gq d,A [m,n] and Gq d, [m,n] which are required for the computation of the discrete cyclic convolution given by Equations (22) and (23). These elements have to be computed for the discrete values 0≤m≤2Nx+1, 0≤n≤2Ny+1. According to the previous section, the values of Gq d,A [m,n] and Gq d, [m,n] for 0≤m≤ Nx, 0≤n≤Ny are the values of the integrals given by Equations (16) and (17) in the points x=xm, y=yn and z=-hq when 0≤xm-xm'<a and 0≤yn-yn'<b. Note that, in similar way, the values of Gq d,A [m,n] and Gq d, [m,n] for the discrete values m> Nx and/or n>Ny are the values of the integrals given by Equations (16) and (17) in the points x=xm, y=yn and z=-hq when xm-xm'>a and yn-yn'>b. Let G0 represent any of Green's functions G A xx and G   The Floquet representation of G0see Equations (12) and (13) in [11] ensures the following periodic properties: where kx0=k0sin(inc)cos(inc), ky0=k0sin(inc)sin(inc) and k0=2π/ is the vacuum wavelength If we take into account these periodic properties in Equations (16) and (17) (16) and (17) (16) and (17) for the observation point (x, y, z)=(xm, yn, -hq) using MRW quadrature rules [17].

If the singularity behaviour of the Green functions is the dominant behaviour, then the values of Gq d,A [m,n] and Gq d, [m,n] can be computed by numerical integration of Equations
The solid line in Figure 2 shows the computed results of the normalized values of Gq d, [m,m] for 0≤m≤Nx=127 for a grounded dielectric multilayer (N=q=4 in Figure 1). The thickness and dielectric constant of each dielectric layer are the following: d4=0.7 mm, r4=2.1, d3=0.3 mm, r3=12.5, d2=0.5 mm, r2=9.8, and d1=0.3 mm, r1=8.6. This multilayer medium was selected to show the capability of the proposed method to analyse periodic structure on a complex multilayer medium in fast way. A similar complex multilayer medium was used in [32]. The periodic structure has a square period given by a=b=10 mm. The results are obtained under oblique incidence given by inc=20 deg and inc=0 deg at 10 GHz. The results shown as the solid line were obtained by means of the computation of Equations (16) and (17) using DE formulas with 203203 quadrature points (i.e., five levels of the quadrature rule [16]). These values are considered virtually exact. Relative error of the results obtained by MRW quadrature rules with respect to the results obtained by DE formulas are also shown as a dashed line. These results are obtained using 33 quadrature points when m≠0 and with 4040 quadrature points when m=0. We can see that the error level obtained by MRW quadrature rules is lower than 0.03% for the worst case (m=0). The relative error of the results obtained by the approximation given in Equation (28) with respect to the results obtained by DE formulas is also shown as the dotted line. These errors reach values lower than 0.16% when mx/or (mx-a)/ However,the error level increases as mx/or (mx-a)/decreases. Therefore, a mixed computation is implemented to obtain enough accuracy: computation of Gq d, [m,n] and Gq d, [m,n] by the approximations given by Equations (27) and (28) when mx/or (mx-a)/and computation by MRW quadrature rules otherwise. The relative error obtained by this mixed computation is also shown as the dash-dotted line.
We would like to point out that the total consumption of CPU time to compute all the values of Gq d,A [m,n] and Gq d, [m,n] for (2Nx+1)(2Ny+1)=256256 is 2.247 s, which was obtained in a laptop computer with a processor Intel Core i7-6700HQ at 2.6 GHz with 32 GB of RAM memory.
Finally, the computation of the inductive contributions given in Equation (12) of the elements of the coefficient matrix is computed by conventional quadrature (for example Gauss-Legendre quadrature rules [31]) where the samples of the components of the vector function fj ind (x,y,z=-hq) are computed from the values of the discrete functions f d,ind j,q,x/y[m,n] by conventional bilinear interpolation [31]. A similar procedure was carried out for the computation of fj cap (x,y,z=-hq) from the elements f d,cap j,q[m,n]. Since a dense equi-spaced mesh is required to approximate efficiently the BFs by pulses expansions given in Equations (9) and (10), the bilinear interpolation provides enough accuracy for practical purposes.

Numerical Results
In this section, we show the results of four cases of studies of surface geometry of layout: (1) two ellipse halves, (2) square loop/patch combination, (3) two concentric split rings and (4) two stacked crosses. Figure 3 shows the definition of the geometrical parameters of the unit cell of the periodic structure which is considered in this work. The unit cell hosts a metallic layout which consists of two ellipse halves on the grounded multilayer periodic medium used in the previous section. This surface geometry of the layout was selected to show the capabilities of the proposed method to analyse the resonant layout with a complex surface geometry. The semi-major and semi minor axis of each ellipse half are given by L1 and L2=2 parameters where L2=2L1/3. The ellipse halves are separated by a fixed gap of w=0.5 mm. The NURBS model of the geometry was made by a periodic structure module of newFASANT [33]. We would like to point out that this layout is modelled by only four NURBS surfaces.  Figure 4 shows equi-spaced mesh points x=xm, y=yn, z=-hq 0<m<Nx, 0<n<Ny that are inside of the periodic cell 1515 mm 2 . These points were used to define the discrete functions described in section 3 given by Equations (18)- (23). Values of Nx=Ny=127 and L1=4.1365 mm are considered. This length value of L1 is the value for which the periodic structure is resonant at 10 GHz under oblique incidence given by inc=20 deg and inc=0 deg. The blue points stands for the layout outer equi-spaced mesh points while that the red points stands for the layout inner equi-spaced mesh points. The centre points of the Bézier patches used in the definition of the BFs are also shown. We would like to point out that the equi-spaced mesh is dense enough to consider the elliptical curves and the aperture between the two ellipse halves.   (22) and (23) associated with BF defined on adjacent Bézier located in the central region of the upper ellipse half (i.e., when j=102). These results were obtained at 10 GHz for L1=4.1365 mm under oblique incidence given by inc=20 deg and inc=0 deg. In these computations, the values of Gq d,A [m,n] and Gq d, [m,n] for (2Nx+1)(2Ny+1)=256256 shown in Section IV were used. Note that the discrete functions f d,ind j,q,y[m,n] and f d,cap j,q[m,n] contain the values of continuous functions fj ind (x,y,z=hq) and fj cap (x,y,z=-hq) of the original problem for equi-spaced mesh points x=xm, y=yn, z=-hq 0<m,n<127 (this region is limited by the dotted grey lines in Figure 5). These results are shown to emphasize that the proposed EPP contains the results of the original problem for equi-spaced mesh points x=xm, y=yn, z=-hq 0<m<Nx, 0<n<Ny.  (22) and (23)  Since a dense equi-spaced mesh is required, samples of the continuous functions f ind j,x/y(x,y,z=-hq) and fj cap (x,y,z=-hq) for 0≤x≤a, 0≤y≤b are efficiently computed from the values of f d,ind j,q,x/y[m,n] and f d,cap j,q[m,n] by conventional bilinear interpolations with enough accuracy for practical purpose.

Two Ellipse Halves
The inductive and capacitive contributions of the elements of the coefficient matrix given in Equations (12) and (13) were efficiently computed from these samples. These inductive and capacitive contributions were computed when i=j for different values of the size of the equi-spaced mesh (2Nx+1)(2Ny+1) at 10 GHz under oblique incidence given by inc=20 deg and inc=0 deg for L1=4.1365 mm. These inductive and capacitive contributions were also computed by direct integration of Equations (7) and (8) using MRW quadrature rules. Figure 6 shows the relative error of sum of both contributions, Zii, with respect to the results obtained by means of MRW quadrature rules with 4040 quadrature points. In these computations, 182 BFs Jj(x',y') (i.e., 0<j<182) for the expansion of the surface density current of (3) were taken into account. We can see that the relative errors decrease as the size of the equi-spaced mesh (2Nx+1)(2Ny+1) increases. Since increasing of the density of the equi-spaced mesh produces a better pulse approximation, these results were expected. Table 1 shows the phase of the copolar reflection coefficient for transverse magnetic polarized incident wave when the unit cell was analysed by the two different methods: the proposed method where Zij was computed by Equations (6), (12) and (13) with Equations (12) and (13) computed by the procedure described in Section III, and the direct MM where Zij was computed by Equations (6)- (8) with the direct evaluation of Equations (7) and (8) by MRW quadrature rules. The result were obtained at 10 GHz under oblique incidence given by inc=20 deg and inc=0 deg for L1=4.1365 mm. Again 182 BFs Jj(x',y') were considered. The CPU time consumption is also shown. These CPU times were obtained with a laptop computer with processor Intel Core i7-6700HQ at 2.6 GHz with 32 GB of RAM memory. We can see that the computation of Zij by means of Equations (6), (12) and (13) produced a gain of CPU time of 12.9 times with respect to that produced by the direct MM when 5x5 MRW quadrature points were used. If we assume that the results provided by the direct MM with 40x40 MRW quadrature points are virtually exact, the phase of the reflection coefficient obtained by the proposed method provides a phase error less than 0.3 deg, while the phase obtained by the direct MM with 5x5 MRW quadrature points provides a phase error of roughly 2.8 deg. We would like to point out that the direct MM takes into account the generalized rooftops BFs defined on pairs of adjacent Bézier patches associated with NURBS surfaces. In this way, the direct MM uses a highorder description of the geometric layout. These results show that although the proposed method used an equi-spaced mesh model of geometry, the equi-spaced mesh is dense enough to provide accurate results with low CPU time consumption.  (6), (12) and (13) with respect to the direct method given by the direct computation of Equations (6)-(8) using MRW quadrature rules with 4040 quadrature points. The relative errors are shown for different values of the size of the equi-spaced mesh (2Nx+1)(2Ny+1). The results were obtained at 10 GHz for L1=4.1365 mm under oblique incidence given by inc=20 deg and inc=0 deg.  Table 2 shows the phase of the copolar reflection coefficient for transverse magnetic polarized incidence when the unit cell was analysed by the two different methods for different numbers Nb of the BFs Jj(x',y') in the approximation of the surface density current given in (3). The results are shown for Nx=Ny=127 in the proposed method and 5x5 MRW quadrature points in the direct MM.
The consumption of CPU time obtained by both methods is shown. Again, these CPU times were obtained with a laptop computer with a processor Intel Core i7-6700HQ at 2.6 GHz with 32 GB of RAM memory. We can see that the CPU time provided by the proposed method increased from 6.72 (with Nb=138) to 29.8 s (with Nb=772) while the CPU time provided by the direct MM shows a larger increase (from 61.4 to 1625 s). This fact means that the CPU time gain provided by the proposed method is between nine and 50 times faster than the direct method. This justifies the efforts carried out in Section III to compute the inductive and capacitive contributions of the MM matrix entries appearing in Equations (12) and (13) in efficient way. We analysed the dependence of CPU time consumption as a function of the number Nb of the BFs for both methods. In this way, we fitted the logarithm of CPU time provided by the direct MM versus Log10(Nb) to a straight line model by least square. The resultant straight line model is Log10(time)=1.955Log10(Nb)-2.389 with coefficient of determination R 2 =0.9962. Thus, we can see that CPU time consumption provided by the direct MM as a function of the number of BFs, Nb, is roughly O(Nb 2 ) as expected. On the other hand, we fitted the CPU time provided by the proposed method versus NbLog10(Nb) to a straight line model by least square. The resultant straight line model is time=0.01189NbLog10(Nb)-3.467 with coefficient of determination R 2 =0.9987. Thus, we can see that CPU time consumption provided by the proposed method as a function of the number of the BFs, Nb, is roughly O(NbLog10(Nb)). Figure 7 shows the CPU time versus Nb for the values of CPU time and Nb shown in Table 2 and the fitted least square models.   Table 2 for both direct and proposed methods. The fitted least square models are also shown.
In the following section, we will show several validations of different periodic structures. Figure  8 shows the phase curves of the copolar reflection coefficient for transverse magnetic polarized incident wave for the unit shown in Figure 3 at 10 GHz under oblique incidence inc=20 deg and inc=0 deg for variation of the L1 parameter. These phase curves were obtained by both methods, the proposed method and the direct method. We can see that excellent agreements were found between both sets of numerical results.  Figure 9 shows the phase curves of the of the copolar reflection coefficient for X-polarization of the Figure 3 of [34] provided by reflectarray element proposed by Zhou et al. under normal incidence. According to [34], the reflectarray element hosts a metallic layout which consists of a square loop/patch combination. The size of the side of square patch is given by L2 while the size of the side of the loop is given by L1. The width of the loop is given by w. As shown in [34], the fixed relations L2=0.69L1 and w=0.135L1 are considered. The periodic lengths are 13.33x13.33mm 2 . The results shown in Figure 3 of [34] were reproduced by the proposed method using (2Nx+1)=(2Ny+1)=256 and direct methods using 5x5 MRW quadrature points. In [34], shaped beam transmit-receive flat reflectarray was designed using this reflectarray element with stringent requirements and a high performance. We can see that good agreements were found between the three sets of numerical results.  Figure 10 shows the phase curves of the copolar reflection coefficient for right handed circular polarization (RHCP) and left handed circular polarization (LHCP) of Figure 7(a) of [13] at 29.75 GHz under oblique incidence inc=30 deg and inc=0 deg. These results are provided by a dual-band reflectarray element which consist of two concentric split rings with 5x5 mm 2 of period lengths. The results provided in [13] wree obtained with the novel sophisticated technique which provides very accurate results with very fast computations. However, this technique is only suitable when the surface geometry of the layout is split rings or single circular arcs. The results were reproduced by the proposed method using (2Nx+1)=(2Ny+1)=256 and direct methods using 5x5 MRW quadrature points. These reproduced results were obtained using the following fixed geometrical parameters 2=150.4 deg, 2=1.85 mm, 1=1.20 mm, w=0.2 mm. The geometrical variables 1 and used are given in Figure 7(b) of [13]. In [13] and [35], this reflectarray element was used to design dual circular polarized and dual band focused reflectarray. We can see that acceptable agreements were found between the three sets of numerical results. Figure 10. The phase curves shown in Figure 7(a) of [35] and the reproduced results by the proposed and direct methods.

Two Stacked Crosses
Finally, Figure 11(a) shows the square of magnitude of the transmission coefficient, |T| 2 , and reflection coefficient, |R| 2 , provided by two stacked crosses printed on a single dielectric layer studied in Figure 4(b) of [36]. The period of the structure is a=b=10 mm and the dielectric layer has a permittivity r=2.58 and thickness of d=2.362 mm. The results were obtained at normal incident for a for transverse electric polarized incident wave with arm lengths and width of the crosses of L=6.875 mm and w=0.625 mm, respectively. The results were reproduced by the proposed method using (2Nx+1)=(2Ny+1)=256 and direct methods using 5x5 MRW quadrature points. We can see that acceptable agreements were found between three sets of numerical results. Figure 11 (b) shows the phases of the transmission and reflection coefficients obtained by the proposed method and direct method. Again, we can see that acceptable agreements were found between the three sets of numerical results.
(a) (b) Figure 11. (a) Square of magnitude of the transmission coefficient, |T| 2 , and reflection coefficient, |R| 2 , of two stacked crosses printed on a single dielectric layer studied in Figure 4(b) of [36]. (b) shows the phases of the transmission and reflection coefficients.

Conclusions
In this work, several approaches to pulse expansion with dense equi-spaced mesh of generalized rooftop BFs were carried out, defined on a pair of adjacent Bézier patches and an EPP which contains the original problem. These approaches led discrete cyclic convolutions instead of usual continuous convolution between Green's functions and BFs obtained by the direct MPIE approach. These discrete cyclic convolutions involve periodic discrete functions of discrete values of integrals of Green's functions and BFs. These discrete cyclic convolutions are efficiently computed by FFT procedure for very dense equi-spaced mesh.
Prior to efficiently computing the discrete cyclic convolutions, a fast and efficient mixed preprocesing technique of computation of integrals of Green's functions was proposed. This efficient pre-processing computation show consumption of CPU time of 2.247 s for 2256256 integrals. Using these pre-processing results, an analysis of a unit cell was carried out. Comparisons of consumption of CPU time by the direct MM and the proposed method was carried out and the comparisons show that the proposed method can be between 10 and 50 times faster than the direct method when the number of BFs increases from 138 to 772. Dependence of CPU time consumption with as a function of BFs was carried out. O(NbLog10(Nb)) behaviour of CPU time consumption by the proposed method was found. This behaviour provides significant CPU time savings with respect to the expected behaviour of O(Nb 2 ) provided by the direct computation of the MM matrix entries. Finally, several validations of the proposed and direct methods were successfully carried out.