Abstract
Linear embedding via Green’s operators (LEGO) is a domain decomposition method particularly well suited for the solution of scattering and radiation problems comprised of many objects. The latter are enclosed in simple-shaped subdomains (electromagnetic bricks) which are in turn described by means of scattering operators. In this paper we outline the extension of the LEGO approach to the case of penetrable objects with dyadic permittivity or permeability. Since a volume integral equation is only required to solve the scattering problem inside a brick and the scattering operators are inherently surface operators, the LEGO procedure per se can afford a reduction of the number of unknowns in the numerical solution with the Method of Moments and subsectional basis functions. Further substantial reduction is achieved with the eigencurrents expansion method (EEM) which employs the eigenvectors of the scattering operator as local entire-domain basis functions over a brick’s surface. Through a few selected numerical examples we discuss the validation and the efficiency of the LEGO-EEM technique applied to clusters of anisotropic bodies.
1. Introduction
It is straightforward to incorporate anisotropic materials into macroscopic Maxwell’s equations by allowing for constitutive relationships with dyadic permittivity or permeability [1]. Nevertheless solving an electromagnetic (EM) scattering or radiation problem which involves objects with dyadic constitutive parameters demands the application of some numerical method.
For instance, the direct solution of Maxwell’s equations in differential form with the finite-element method (FEM) [2] is tempting, because the latter procedure leads to highly sparse matrices. On the other hand, radiation boundary conditions are naturally taken into account by means of volume integral equations (VIEs) [3–6], but their numerical treatment with the Method of Moments (MoM) produces dense matrices (e.g., [7, 8]). In either case, when the structure of concern is electrically large, one could apply an iterative technique to solve the resulting algebraic linear system. However, the convergence may be slow, unless the matrix is well conditioned; poor conditioning seems to especially plague problems that involve objects with large contrast with respect to the background medium [9].
Solution strategies based on domain decomposition methods (DDMs) may help overcome the foresaid drawbacks, because rather than tackling a given EM problem as a whole, DDMs handle the local strong interactions and the distant weaker ones in a large composite structure separately. This is accomplished by enclosing “small” parts of a large structure into subdomains that are treated independently of one another. The multiple scattering that occurs among the various subdomains is then described by means of translation or transfer operators. In this way, the usage of VIEs is restricted to the solution of the EM problem within the subdomains and hence, the resulting “smaller” matrices are likely less prone to ill conditioning. Finally, a surface integral equation (SIE) for the unknown fields or equivalent currents on the subdomains’ boundaries is stated and solved numerically. Among the methods that follow this paradigm we mention the equivalence principle algorithm (EPA) [10] and the tangential equivalence principle algorithm (T-EPA) [11], the generalized surface integral equation (GSIE) method [12], and the linear embedding via Green’s operator (LEGO) technique [13, 14].
The domain decomposition in LEGO is effected by enclosing the objects into 3D simple-shaped domains (Figure 1) which constitute the building EM “bricks” of the model. Firstly, each brick is described through a scattering operator of equivalent surface current densities defined over a brick’s surface. Secondly, the multiple scattering occurring between any two bricks is captured by means of transfer operators. Since the mathematical framework for the separation is provided by the Surface Equivalence Principle (SEP) [1, 15], LEGO yields a system of as many coupled SIEs for the unknown surface current densities as the number of bricks in the model [14, 16].
Now, handling anisotropic media with LEGO is straightforward, in that there is no limitation as to the constitutive parameters of the objects inside a brick. Needless to say, the calculation of the relevant scattering operator requires the solution of the EM scattering problem inside the brick by means of a VIE, as anticipated. And yet, once the scattering operators of all the bricks in the model have been computed, the usual equations of LEGO can be stated and solved as before [14, 16, 17].
The important point is that a VIE has to be solved only for one object at a time, if many different objects form the cluster, or just once, in case the objects are all equal to each other. Furthermore, in the context of the numerical solution with the MoM, the (surface) basis functions on a brick’s boundary are usually far fewer than the (volume) basis functions within the object. Therefore, it is expected that the LEGO procedure can reduce the computational burden as compared to the workload that would be required to solve a single VIE for all the objects simultaneously. Additionally, substantial reduction of the degrees of freedom is achieved thanks to the eigencurrents expansion method (EEM) in which the eigenvectors of the scattering operator are employed as local entire-domain basis functions over a brick’s surface [14, 18, 19].
The rest of this paper is organized as follows. In Section 2 we outline the derivation of the scattering operator of a brick that embeds an anisotropic object, and we give detailed expressions for the special case of dyadic permittivity and scalar permeability. Besides, we review the LEGO equations involving the total inverse scattering operator of a composite structure. In Section 3 we describe the numerical solution of the LEGO equations with the MoM and the order reduction through the EEM. Lastly, in Section 4 we present an example of validation, a study of the convergence properties of the EEM, and a numerical example that involves a complex multiscale structure comprised of a magnetized plasma. Efficiency in terms of matrix sizes and computational times is briefly addressed as well.
A time dependence for fields and sources in the form is assumed and suppressed throughout.
2. Formulation with LEGO
2.1. Generalities
Simply put, the solution of a scattering problem which involves an aggregate of anisotropic bodies consists of three main steps [14, 16, 17]: (i) definition of the EM bricks , and calculation of the scattering operators thereof (Figure 1); (ii) calculation of the transfer operators , , , and , which enter the total inverse scattering operator of the structure; and (iii) numerical inversion of resulting surface integral equation.
In this section, we elaborate on steps (i) and (ii), whereas we will deal with step (iii) in Section 3. At this stage, we make no hypothesis with regard to the position, the shape, and the constitutive parameters of the bodies. Nevertheless, in the formulation to be developed hereinafter we consider bricks of identical shape. The latter assumption—which should not constitute an issue especially for arbitrary displacements of objects and bricks—greatly facilitates the numerical solution, as will be clear in Section 4.
2.2. EM Bricks with Anisotropic Media
Sketched in Figure 1 is a quite general example of EM brick immersed in an infinite background homogeneous medium (labelled ①). The brick comprises a simple-shaped domain filled with an isotropic homogeneous host medium (tagged ②) which in turn wraps the anisotropic body (tagged ③). For the sake of completeness, we point out that in [14] we considered the case of PEC or penetrable isotropic objects and media ① and ② with the same constitutive parameters, whereas in [16] we introduced the notion of scattering operator of a material interface (i.e., a brick’s boundary, ) which allows handling the instance of different media ① and ② efficiently.
To keep the exposition lucid, with the aid of Figure 2 we first outline the general procedure for deriving the scattering operator in the case of background and host medium having the same EM properties. Then, we specialize the equations and give the explicit form of the operators involved for the instance of medium ③ with dyadic permittivity and scalar permeability . (The case of medium ③ with dyadic permeability and scalar permittivity is similar and ensues in a straightforward manner.) Moreover, the scattering operator of a material interface along with the cascading procedure described in [16] can be employed to obtain for the general situation of Figure 1, though we will not examine this configuration in what follows.
(a)
(b)
(c)
(d)
With these preliminaries we are now ready to elaborate on the derivation of .
The starting point is a solitary brick illuminated by incident fields , , radiated by some external sources. An application of the SEP [1] allows replacing the effect of the external sources with equivalent electric () and magnetic () surface current densities flowing on [14]. The latter are collected in the abstract vector for ease of manipulation; namely, with being the intrinsic impedance in medium ①. The superscript “” serves to remind us that the currents in question reproduce the incident fields towards the inside of the brick (see Figure 2(a)). We notice in passing that and can be thought of as being positioned on either side of , because the brick’s boundary is not a material interface by hypothesis (cf. [16]).
To proceed we observe that the object inside is illuminated by the EM field generated by . As suggested in Figure 2(a), such field can be related to its sources symbolically by means of a linear operator; that is, where is the abstract vector and the propagator is a matrix of dyadic surface integral operators whose kernel is the dyadic Green’s function of media ① and ②.
Now, by invoking the Volume Equivalence Principle (VEP) [1] we replace the object with volume current densities [20]. In the general case of both electrically and magnetically penetrable objects, said currents are related to the total flux densities and (e.g., [7]). The latter can be obtained by solving a set of VIEs; namely, where and is a matrix of volume integral operators over the region occupied by medium ③.
Next, we observe that the volume “currents” generate scattered fields on the boundary , as shown in Figure 2(c); in symbols, this reads where is the abstract vector and the propagator is a abstract matrix of dyadic volume integral operators involving dyadic Green’s function of media ① and ②; besides, the subscript “” stands for “tangential.” The superscript “” (short for “scattered”) signifies that the fields of concern here propagate away from the object and outwards the brick.
Lastly, as is graphically explained in Figure 2(d), a second application of the SEP over allows substituting the scattered fields with equivalent surface electric () and magnetic () current densities [14]. In symbols, with being an abstract vector defined in a similar fashion to (1). The propagator is exactly the same as in [14, Table 1], since it depends only on the shape of and not on the media inside .
By formally solving (2), (4), (6), and (8) for as a function of we find (this expression differs from [14, Equation ] for a minus sign, which is due to different definitions adopted for the operator ) which provides us with a suitable recipe for computing the scattering operator .
In the notable instance of an anisotropic dielectric object with , the operators , , and reduce to , , and abstract matrices whose explicit expressions read with being the wavenumber in medium ①, being the distance between source () and observation () points, and and denotes the surface unit dyadic on . The null vector in (7) and the null dyadic appearing in in (12) are a consequence of the definition of the propagator which relies only on the magnetic field to yield the corresponding equivalent surface current densities on (cf. [14], Table 1).
We conclude by observing that the propagators , and the operator are dimensionless, as can be ascertained by inspection of (10), (12), and (11). In contrast, current and field vectors , , , and are all normalized so as to carry the physical dimension of .
2.3. Total Inverse Scattering Operator and LEGO Equations
When two or more bricks are considered, (9) must be modified to account for the multiple scattering that occurs between the bricks. This is readily accomplished upon recognizing that the equivalent scattered currents on produce fields that reach the boundary of alongside the incident fields due to the independent sources of the EM problem. By invoking the SEP, the extra incident fields can be turned into equivalent surface current densities on . In symbols, where the propagators , are the same as in [14, Table 1], as they depend only on the shape of . The formal inversion of (14) yields where the dimensionless transfer operator maps scattered currents on onto incident currents on [14].
In the light of the previous discussion and thanks to (15), the extension of (9) to the case of interacting bricks reads [14] which represents a system of coupled surface integral equations. Upon grouping all the scattered and incident currents column wise into and , (16) can be expressed in a compact form as with where is the total inverse scattering operator of the structure. Equation (17) along with (19) constitutes the formulation of the EM problem.
Notice that (19) is a formal result not only because ordinarily and its inverse cannot be obtained in closed form, but also because may not exist per se when happens to be singular. In spite of these complications, the numerical calculation of through (19) is viable and stable, as long as we carefully handle the possible null or small eigenvalues of [14, 18], as we show in Section 3.2.
3. Numerical Solution with MoM and EEM
The numerical solution of (17) with the Method of Moments (MoM) in Galerkin’s form and the eigencurrents expansion method (EEM) was extensively presented in [14, 18]. For ease of reference we briefly recall the key points of the procedure herein below.
3.1. Reduction to a Weak Form with MoM
We employ the MoM with subsectional basis functions and a symmetric inner product repeatedly in order to (a) turn (2) and (6) into a weak form; and (b) formally invert (4), (8), and (14). In particular, we model the bricks’ boundaries and the volume of the anisotropic objects with 3D triangular-faceted tessellations and tetrahedral meshes, respectively. Then, we associate () surface (volume) divergence-conforming linear vector functions [21, 22] with the surface (volume) mesh on (inside) each brick in order to expand the surface currents and (the flux density or ).
As a result, we can write the algebraic counterparts of and as where and have size and the rank of is .
Upon using (20) and (21) as building blocks, the algebraic counterpart of (17) can be written formally as where the vectors store the MoM expansion coefficients of . As anticipated, though, explicitly building and then numerically inverting the linear system (22) are far from being trivial for two reasons.
First and foremost, the inverse scattering operators may not be defined in some cases. Specifically, if the volume basis functions in the tetrahedral mesh are fewer than the number of surface basis functions on , then, as is clear from (20), the scattering operator is singular, with its rank being . At any rate, even though , may still be nearly singular due to the presence of a few very small eigenvalues. On the whole, inverting is not only inconvenient, but also a potentially ill-conditioned problem.
Secondly, difficulties posed by the inversion of aside, the size of the total inverse scattering operator, may be quite large, if one considers a problem modelled with a moderate to large number of EM bricks. As a consequence, inversion with direct solvers such as LU factorization [23] may not be feasible. Then, iterative solvers seem an obvious alternative, but the convergence rate is expected to be poor, since is ill conditioned, in the light of the previous observations.
To circumvent the hurdles above we apply the EEM [14].
3.2. Compression and Inversion with Eigencurrents
In the EEM we employ the eigenvectors of as local entire-domain basis functions to expand the scattered currents on . Since these eigenvectors in tandem with the low-level basis functions define surface currents over , we refer to the columns of as to the eigencurrents of .
The eigencurrents can be separated into two subsets: coupled and uncoupled. The former—associated with the larger eigenvalues of —produce substantial radiation and thus contribute to the multiple scattering that occurs among the bricks. Conversely, the uncoupled eigencurrents—which correspond with the smaller or null eigenvalues of —yield scattered fields that are practically confined around a brick and hence serve to describe the scattered currents on a brick. We denote by the number of coupled eigencurrents.
The notion of coupled-uncoupled eigencurrents allows us to partition the matrix (expressed in the basis of the eigencurrents) into four subblocks [14, Equation ]; namely, where the subscripts “” (“”) stand for coupled (uncoupled).
In the spectral reordered representation (23) of the total inverse scattering operator the subblock —relevant to the interaction between any two uncoupled eigencurrents from either two distinct bricks or the same brick—can be approximated with its diagonal, that is, the reciprocal of the eigenvalues associated with all the uncoupled eigencurrents, because the latter interact only with themselves. Likewise, the off-diagonal subblocks , can be zeroed altogether, as they account for the negligible interaction between a pair of coupled and uncoupled eigencurrents. Lastly, the diagonal subblock —relative to the interaction of any two coupled eigencurrents from either two distinct bricks or the same brick—has to be retained in full.
In the end, the inverse of is approximately given by where is a diagonal matrix that contains the eigenvalues germane to the uncoupled eigencurrents of the problem. Since is known the moment the eigenvalues of have been determined, evidently, the occurrence of small or null eigenvalues does not pose any problem at all in (24). Furthermore, since the number of coupled eigencurrents contributed by is usually far smaller than the number of low-level basis functions over , the calculation of requires inverting a relatively small matrix—which can be carried out with direct solvers.
In conclusion, the EEM facilitates the inversion of , as anticipated. Once the expansion coefficients in the eigencurrents basis have been computed, can be retrieved through .
4. Results
4.1. Example of Validation
An extended numerical code has been developed for the calculation of and and the subsequent solution of (22) with the EEM in the instance of objects with either dyadic permittivity or dyadic permeability.
Among many other tests, as a sanity check we have examined the problem of a solitary dielectric sphere with scalar permittivity and embedded it in a cubic brick . We have computed the scattering operator with the newly developed code and the old one [14], which employs surface integral equations (e.g., the Poggio-Miller-Chang-Harrinton-Wu-Tai formulation (PMCHWT) [24]) at the interface between media ② and ③. With the two codes we have obtained the equivalent scattered currents on and the scattered fields in the Fraunhofer region when the sphere is illuminated with a plane wave. The number of surface basis functions over is = 1152; the number of basis functions for the solution of the scattering problem inside is for the PMCHWT approach and for the approach based on (4) and (11). (Geometrical and physical data are given in the caption of Figure 3.)
| (a) LEGO/SIE |
| (b) LEGO/VIE |
| (c) LEGO/SIE |
| (d) LEGO/VIE |
The total equivalent surface current densities and yielded by the two codes are plotted in Figure 3; as can be seen, the current distributions are in excellent agreement. This notable result is a consequence of the algebraic scattering operator being theoretically independent of the actual formulation of the EM problem inside . Besides, if the scattered currents are correct, then we may expect also the far-field quantities (which from are derived) to be in excellent agreement. The latter is confirmed by the plot of the bistatic radar cross-section (RCS) given in Figure 4.
As a final remark, we notice that the solution with VIE will usually require more basis functions than what is necessary to solve the SIE (when the latter makes sense). This is obviously due to the need for distributing basis functions inside a volume rather than a surface. In fact, our numerical tests show that in order to attain the same level of accuracy of LEGO/SIE and LEGO/VIE, the size of the triangular facets over the sphere’s surface (see insets of Figure 4) should be comparable. If that criterion is met, then the requirement for having a uniform mesh in somehow constrains the number of tetrahedra and hence basis functions in the object.
4.2. Convergence with Number of Coupled Eigencurrents
To test the efficiency of the EEM [14] in compressing the algebraic system (22) we have considered the plane-wave scattering from a linear arrangement of anisotropic spheres. The EM problem has been modelled with cubic bricks each containing one sphere, as shown in Figure 5. The number of surface basis functions over is = 1152, whereas the number of volume basis functions within the sphere is . (Geometrical data and constitutive parameters are given in the caption of Figure 5.)
Notice that in the LEGO setup for this problem, only needs to be computed, as the four bricks are identical. Furthermore, thanks to the limited translational symmetry of the structure in Figure 5 only transfer operators out of are independent and must be computed. This would not be the case, if the bricks had different shapes.
We have obtained a reference solution by solving (22) with no reduction at all, that is, with coupled eigencurrents. Subsequently, we have applied the EEM to its fullest, with coupled eigencurrents, and have solved the problem by using the compressed scattering operator (24). The RCSs computed with are plotted in Figure 6 for comparison.
It can be seen that coupled eigencurrents are already sufficient for the solution to capture the trend and most details of the RCS versus the elevation angle ; in particular, position of nulls and local maxima are predicted, even though the levels thereof depart from the reference values. However, with eigencurrents the RCS () is practically indistinguishable from the reference curve (−). From a numerical viewpoint this means an accurate solution can be achieved by inverting a matrix of rank as compared to the rank of the original uncompressed system (22).
We elaborate a bit further on this numerical experiment as regards the assumptions that constitute the groundwork of the reduction scheme based on the notion of coupled-uncoupled eigencurrents. In Section 3.2 we posit that the coupled (uncoupled) eigencurrents correspond with the larger (smaller) eigenvalues , , of the (algebraic) scattering operator . As an example, for the problem of Figure 5 the spectrum of is plotted in Figure 7. It can be ascertained that the eigenvalues decrease exponentially to zero (i.e., the threshold of numerical noise in double-precision floating point arithmetics). In fact, and differ by approximately one order of magnitude. In the light of the conspicuous, fast decay of the spectrum, it can be understood why the eigencurrents associated with the eigenvalues , contribute naught to the fields scattered by and hence can rightly be termed uncoupled.
In [18] we demonstrated that for isotropic objects the spectrum of is mostly affected by the constitutive parameters of medium ③ (if media ① and ② are the same), the size and the position of the object with respect to ; a weaker dependence on the working frequency was observed. For the time being, we mention that it seems reasonable to expect similar trends also in the case of dyadic constitutive parameters of medium ③, although a thorough and extensive numerical campaign is required to support this supposition. (Also see Section 4.3.)
4.3. Example of Application: Arrays of Microplasma Rods
We now apply LEGO and EEM to solve a relatively large and multiscale problem, namely, the scattering from arrays of microplasma rods; these structures appear to have found application as reconfigurable metamaterials (e.g., [25]).
We consider four arrays comprised of , , , and -aligned cylindrical rods arranged in a Cartesian lattice in the plane; shown in Figure 8 are the tetrahedral model of a rod and the 3D surface mesh of the cuboidal embedding brick. The lattice constants in both and directions are the same as a brick’s width. (See the caption of Figure 8 for applicable geometrical and mesh data.) Modelling the small cylindrical rods requires a reasonably fine mesh capable of accurately reproducing the small geometrical features. By contrast, the surface mesh on a brick’s boundary need not be as fine because the brick has got a simple shape. Thereby, we anticipate that the bare introduction of LEGO bricks can already afford reduction of unknowns.
(a)
(b)
The plasma is a collisional magnetized gas of electrons with density , -aligned static confining magnetic intensity mT, and collision frequency , where denotes the electron (angular) plasma frequency. With these values and for frequencies GHz, the plasma dyadic permittivity in Cartesian coordinates reads where the Stix parameters , , and [26] are given explicitly in Table 1. Finally, as the plasma rods are immersed in free space, then .
By applying the EEM to (22) with coupled eigencurrents we have computed the equivalent scattered currents and afterwards the forward and backward RCS for the four arrays of plasma rods. The results are plotted in Figure 9 as a function of the number of parallel rows along the -axis for an -directed -polarized incident plane wave. (See the caption of Figure 9 for relevant data.) The minimum of reflection appears to occur for the array of rods at GHz, which is well below the electron plasma frequency GHz.
The spectrum of for the brick in Figure 8 at the three frequencies of interest is plotted in Figure 10; the vertical dashed line (at ) separates the regions of eigenvalues corresponding with the coupled (left) and the uncoupled (right) eigencurrents. Based on the convergence properties examined in Section 4.2 and the decay pattern which is not dissimilar from the one shown in Figure 7, we can assume that coupled eigencurrents are sufficient to determine accurate far-field quantities.
Besides, the plot of Figure 10 seems to confirm that the spectrum of a scattering operator is only weakly dependent on frequency [18] even in the case of embedded anisotropic objects. In particular, the first eigenvalues are seen to become smaller in magnitude as the frequency increases. This behavior can be explained intuitively by observing that on average the electrical distance between a rod and its surrounding brick increases as the frequency does so. Now, since the scattered field decays away from the object that has generated it and in the basis of the eigencurrents the scattered field (for a single brick) is given by , then it follows that must be smaller if the distance from the brick’s surface and the included object is larger, everything else unchanged.
4.4. Efficiency
We conclude this section by discussing the time requirements of LEGO-EEM with the aid of the problems of Figures 5 and 8. For the sake of completeness we mention that the simulations for the spheres (Figure 5) have been run on a PC equipped with an Intel Core i7-3770 CPU 3.40-GHz processor and 8-GB RAM, whereas the simulations for the microplasma rods (Figure 8) have been run on a PC equipped with an Intel Core Quad CPU 2.66-GHz processor and 4-GB RAM.
In the extension of LEGO to anisotropic bodies the calculation of the scattering operator (20) appears to be relatively time consuming, owing to the need for computing volume integrals in tetrahedra and (formally) inverting , which is usually larger than that in case of SIEs. Therefore, the time taken to obtain appears to be roughly proportional to .
For the problems at hand, is about 49 min for the anisotropic sphere case and about 59 min for the plasma rod case. Nevertheless, the calculation of (and its spectral decomposition) has to be performed only once, if the structure can be modelled by means of identical bricks, as we have shown. This is especially useful if one plans on reutilizing the bricks over again in other simulations, as long as the frequency is not altered, of course.
When this expedient is applied to the case of the anisotropic spheres with different numbers of coupled eigencurrents, the overall calculation time actually becomes , which drops drastically, as can be ascertained from Table 2.
Likewise, we have to compute only once () for the four arrays of microplasma rods. The calculation times are listed in Table 3. Interestingly, grows only linearly with the number of bricks , as shown by the plot in Figure 11. This result is a consequence of the limited translational symmetry of the arrays; in practice, the transfer operators and come in clusters of identical specimens and hence it is sufficient to compute only one element of each cluster [14, 17]. Then, as long as the number of clusters of transfer operators grows linearly with , so does the time . It is worthwhile noticing, though, that if one examines different EM problems involving a fixed number of objects with increasingly larger size (i.e., ever more unknowns ), the complexity—being dominated by the calculation of —will rather grow quadratically.
Finally, the matrix sizes relative to the problem of arrays of microplasma rods are listed in Table 3: is the size of algebraic system that one would need to invert if the problem was formulated directly with a single VIE; is the size of the LEGO system (22); is the size of in (23). It is apparent that the LEGO approach is already competitive over the baseline VIE, whereas further and substantial reduction is afforded by the EEM.
5. Conclusion
We have shown that the extension of the LEGO procedure for solving the EM scattering from clusters of anisotropic bodies basically requires modifying the way the scattering operators are computed. In fact, as long as the intrinsic modularity of LEGO is properly exploited, then incorporating new types of EM bricks into an existing code is in principle straightforward. Besides, modularity implies that when calculations are suitably organized, the scattering operators can be reutilized in subsequent runs, thus enabling, for example, the efficient analysis of increasingly larger structures.
Along the same lines, the EEM can be applied as before because its implementation is independent of the specific types of EM bricks being considered. In this regard, we have demonstrated that when anisotropic objects are enclosed in EM bricks, the EEM performs as efficiently as it does in the instance of conducting or penetrable isotropic media [18].
All in all, the features above make the LEGO method an attractive technique for the numerical simulations of moderately large and yet finite-size structures comprised of many separate objects.
We conclude by observing that although we have presented results for homogeneous anisotropic media, the procedure described in Section 2.2 applies as well to the case of spatially dependent permittivity or permeability without any modification [20].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors would like to thank the reviewers for their suggestions.