On the Use of a Localized Huygens’ Surface Scheme for the Adaptive H-Refinement of Multiscale Problems

This work proposes a domain decomposition method (DDM) based on Huygens’ equivalence principle to efficiently perform an adaptive h-refinement technique for the electromagnetic analysis of multiscale structures via surface integral equations (SIEs). The procedure starts with the discretization of the structure under analysis via an initial coarse mesh, divided into domains. Then, each domain is treated independently, and the coupling to the rest of the object is obtained through the electric and magnetic current densities on the equivalent Huygens’ surfaces (EHSs), surrounding each domain. From the initial solution, the error is estimated on the whole structure, and an adaptive h-refinement is applied accordingly. Both the error estimation and the adaptive h-refined solution are obtained through the defined EHSs, keeping the problem local. The adaptive h-refinement is obtained by a nonconformal submeshing, where multibranch Rao–Wilton–Glisson (MB-RWG) basis functions are defined. Numerical experiments of multiscale perfect-electric-conductor (PEC) structures in air, analyzed via the combined field integral equation, show the performance of the proposed approach.


I. INTRODUCTION
C OMPUTATIONAL electromagnetism problems require a good mesh quality for an accurate and efficient solution.
In general, the mesh density leads to a higher accuracy of the numerical solution, but the increase of unknowns requires high computational resources.A tradeoff between accuracy and computational cost is usually faced in large multiscale problems.One would like good accuracy in the whole object, independent of its smoothness or dimension while saving computational resources.
Several techniques, such as the multilevel fast multipole algorithm (MLFMA) [1], [2], the fast Fourier coupling matrix compression for MLFMA (MLFMA-FFT) [3], or the domain decomposition method (DDM) [4], [5], [6], [7], have been proposed as efficient ways to face the multiscale, large problems.These powerful techniques are usually applied to general-purpose meshes, which could have a high number of unknowns, maybe keeping a fine mesh in regions of the object where it is not needed.The regions where a fine mesh is needed depend on both the geometrical details and the behavior of induced currents, determined by the excitation.These dependencies lead the computational community to look for adaptive techniques that find an optimal mesh for a specific situation, reducing the number of degrees of freedom.These techniques should also be automatic to give optimal discretization without user intervention.Some methods focus mainly on the refinement due to the geometry, as in the case of scatterers with sharp edges [8], [9].Other techniques try to find indicators of the discretization quality in two different ways.The first approach is to check the accuracy of the solution via an error estimation residual: Bibby and Peterson [10] and Saeed and Peterson [11] used overdetermined systems via the electric field integral equation (EFIE) for h-and p-refinements, while, in [12] and [13], the jump of the charge and current at cell boundaries is evaluated.A second approach, known as goal-oriented adaptivity, checks the accuracy and efficiency of a quantity of interest [14], [15].
Tobon Vasquez et al. [16] proposed a local mesh refinement technique, for triangular meshes, in the method of moments (MoM) applied to surface integral equations (SIEs).The main goal was to obtain an automatic technique to discretize with a fine mesh only the regions that needed to be refined, after estimating the error of an initial coarse solution.The local refinement led to mesh nonconformity that was faced through half Rao-Wilton-Glisson (RWG) basis functions and discontinuous Galerkin (DG) [17], [18], [19].The complete problem was first solved for the initial coarse mesh, then a local refinement was implemented to improve the solution accuracy via an error estimation procedure, and, finally, the new discretization was used to compute the solution.Although the error estimation step dealt with a projection on a reference fine mesh, with no need for system matrix inversion, it still involved the filling of a large MoM matrix, coupling far regions with fine meshes.Moreover, in the final problem, filling a large matrix, related to the h-refined mesh, and its inversion were needed.Hence, although the mesh refinement was a localized operation, both the error estimation and final solution were performed globally, due to the coupled nature of the problem, with an evident additional computational cost.
In this work, the numerical inefficiencies of [16] are overcome by a DDM based on Huygens' equivalence principle to perform the proposed adaptive h-refinement procedure.The basic idea is to isolate the different parts of the structure under analysis and perform the error estimation, the mesh refinement, and the final solution, locally.To isolate each domain of the problem, we propose using Huygens' equivalence theorem [20], where properly defined equivalent surfaces will act as an intermediary between each domain and all the rest of the object.This proxy, generated at the initial coarse mesh, is used in both the error estimation and the final problem solution, with its main advantage being in the latter.Moreover, the cell-based nature of the error estimation and refinement leads to a mesh nonconformity that, here, is faced with the div-conforming multibranch RWG (MB-RWG) basis function, recently proposed in [21], with no need for extra penalty conditions on the SIE kernel as in [16].The analysis in this article is focused, without limiting the generality, on the combined field integral equation (CFIE) of perfect-electricconductor (PEC) structures in a homogenous medium (i.e., air).
This article is organized as follows.Section II presents the used SIE formulation.Section III describes all the steps of the proposed approach using Huygens' surfaces and the MB-RWG basis functions.Finally, Section IV shows and discusses some numerical results employing the proposed approach to 3-D geometrically complex structures.Preliminary numerical results have been presented in [22] and [23].

II. BACKGROUND FORMULATION
An electromagnetic scattering problem is considered here for a PEC object in a homogeneous medium.We define J as the equivalent electric current density on the object surface and obtain the tangential EFIE (T-EFIE) and the normal magnetic field integral equation (N-MFIE) by applying the equivalence principle to the total electric and magnetic fields as where η is the intrinsic impedance of the background medium, n is the unit vector normal to the surface, and E i and H i are the incident electric and magnetic fields, respectively.The operators L and K are defined as where k is the wavenumber, r and r ′ are the observation and source points, respectively, ∇ ′ is the divergence in the source coordinates, and P V is the principal value of the integral in (4).The homogeneous Green's function g(r, r′) is defined as Combining ( 1) and ( 2), we obtain the CFIE where where α is a weight, between 0 and 1, that balances the EFIE and MFIE equations (α = 0.5 here).The MoM procedure is applied to (7), with the equivalent electric current densities approximated as a sum of N known vector basis functions where I n are the unknown coefficients.Applying the Galerkin testing procedure, we obtain the following dense matrix system: with where [Z EFIE ] and [Z MFIE ] are N × N matrices, formed by the coupling between all the basis and testing functions, via the EFIE and MFIE operators, respectively.[I ] is an N -column vector collecting the unknown coefficients of the current density expansion, and [V ] is the N -column excitation vector obtained from the excitation fields.
Equation (7) refers to the radiation of the combined (electric and magnetic) fields from the electric current density, J, in the case of a PEC object; hence, we can call it the combined operator for electric currents.However, the procedure described in Section III also requires the definition of the radiation of the combined fields from the magnetic current density defined over the equivalent Huygens' surfaces (EHSs).Then, we define the combined operator of the magnetic currents as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

III. HUYGEN'S-BASED DOMAIN DECOMPOSITION AND ADAPTIVE H-REFINEMENT
This section describes how to integrate a Huygens' surface DDM (HS-DDM) into the adaptive h-refinement technique, going toward a localized approach.In this way, both the error estimation procedure and the h-refined problem solution are performed as local operations.

A. Initial Mesh Solution
We start from an initial coarse mesh that describes the geometry of our problem of interest, with an average discretization size of λ/10.The solution of this initial mesh can be obtained by applying a DDM, as detailed in [6].A convenient, but not mandatory, choice is to keep the same decomposition in P domains throughout the h-refinement process.The initial problem can be expressed as where [Z 0 ] is the system matrix, [V 0 ] is the right-hand side (RHS) vector, and [I 0 ] is the vector collecting the unknown current density coefficients, as defined in ( 9) and (10), here applied to the initial mesh.
The current density coefficients vector, [I 0 ], will be the seed for the localized error estimation process (see Section III-C).It will also be used in the final adaptive h-refined mesh solution (see Section III-E) via the definition of Huygens' surfaces (see Section III-B) and the adaptive h-refinement (see Section III-D).Considering that each domain can be treated in parallel and independently, we will avoid the subscript p when talking of a generic domain.Moreover, the subdivision of the structure under analysis, defined here, will be used through all the following adaptive h-refinement processes.The solution of this initial system can also be obtained through other fast MoM methods, such as the MLFMA-FFT [3]; however, the subdivision into domains is still required for the following steps.

B. Domain Decomposition Through Huygens' Surfaces
After the initial solution and the definition of the domains, an EHS is defined around each domain.The EHS is based on Huygens' equivalence principle, where a radiation problem is substituted by an equivalent problem using an imaginary closed surface and fictitious electric and magnetic surface current densities on it, as graphically shown in Fig. 1.This surface acts as a proxy between the enclosed domain and the other regions of the object, isolating the domain and allowing its independent treatment.
Each EHS is discretized with triangular flat patches with an average cell size around λ/10 [20].In this work, we used spherical surfaces, but any convenient surface can be used as EHS.Fig. 2 depicts an example of a discretized EHS around a domain of interest.As it is evident from Fig. 2, the EHS intersects the object.From the intersection between the EHS and the object, we define four regions, shown in different colors in Fig. 2. The first region is the domain itself, shown in blue.The second region, called "inner buffer" and shown Fig. 1.Huygens' equivalence principle scheme; J S and M S are the electric and magnetic equivalent current densities, respectively.Fig. 2. Example of a discretized EHS and object under analysis, highlighting the considered domain, the buffer regions, and the far-field region associated with that domain.in pink, is the region adjacent to the domain, not belonging to it but still inside its EHS.The third region is called "outer buffer" and is shown in yellow: it includes the cells outside the EHS, but yet near to it.The fourth region is called the "far-field region" (shown in pale green), and it is the zone out of the EHS and far from it.In the following, both buffer regions are excluded from the generation of the equivalent current densities on the EHS as they are in their near field.This leads to a size of the buffer regions of, at least, λ/2 to guarantee that the radiated fields from the surface current in the far-field region via the Huygens equivalent surfaces to the subdomain under consideration are dominated by the traveling waves [20].
The equivalent current densities are obtained on each EHS from the electric and magnetic fields, E S and H S , respectively, radiated from the far-field region of the corresponding domain p, as where np is the inner normal to the EHS around the domain p.These current densities are discretized as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
with f S p,i being the ith RWG basis function, for i = 1, . . ., K p , defined on the EHS around the pth domain.C J p,i and C M p,i are the coefficients obtained from the projection of the continuous current densities on the basis function space.These coefficients are collected in an array of dimension 2K p defined as and they will be used in the following for the error estimation and adaptive h-refined problem solution to represent the coupling of the domain p with its far-field region.

C. Error Estimation Using HS-DDM
The goal of the error estimation is to find a metric that indicates where the initial mesh and, from there, the basis functions defined on it cannot accurately describe the solution current variations of the object under analysis.As in [16], the error estimation used here is residual-based, as it compares the residual of the initial solution with that of a different basis function space, defined on a finer mesh.The Calderon identity [24] shows that the error in the approximate solution is bounded above and below by the residual error of the discretized SIE.Then, we can improve the accuracy of the solution by refining the mesh in the areas where the residual is higher.
The residual is defined as where C J is the CFIE operator as defined in (7) and J 0 is the initial solution current density.The testing of (18) on the same basis functions' space used for obtaining J 0 will lead to zero, due to the orthogonality that was exploited to obtain that solution through Galerkin.Instead, if we use a different set of basis functions as testing, the error can be estimated through the norm of the residual.As we are going toward an h-refinement, a possible choice is to select a set of basis functions that gives information on a more refined mesh.Hence, multiple dyadic divisions of the initial mesh are performed, leading to a finer discretization where the new basis functions are defined.This fine mesh is called the following "reference mesh."The choice of which basis functions can be defined in the reference mesh is arbitrary.However, it is convenient to have basis functions whose linear combination can describe exactly the basis functions defined on the initial mesh.According to this consideration, we chose to test RWG basis functions defined on the reference mesh.
The norm of the residual, R, is estimated through its normalized projection on the RWG basis functions, f L , defined in the reference mesh (ℓ = L), with L being the number of dyadic subdivision from the initial mesh (ℓ = 0) to arrive to the reference one.Hence, in each mth cell of the reference mesh, the residual norm is estimated as where f L m(i) , with i = 1, 2 and 3, are the three RWG basis functions defined on the edges of the selected triangle m of the reference mesh.Then, the cell m estimated residual error is equal to the highest estimated error among the three RWG basis functions defined on the cell m edges.The denominator in (19) corresponds to the norm of each f L m(i) basis function and can be evaluated in closed form.
Substituting (18) into (19) and removing the index m(i) to simplify the notation, the inner product ⟨f L , R⟩ in the numerator of ( 19) can be expressed as To evaluate (20), the current density, J 0 , obtained with the initial mesh, is projected on the f L j basis functions defined on the reference mesh, obtaining where I 0,L j are the current density coefficients obtained with the initial mesh (ℓ = 0) and projected on the N L basis functions defined on the reference mesh (ℓ = L), following the procedure described in [25].We collect them in the vector Substituting ( 21) into (20), we obtain the following matrix system: where [R L ] is the projected residual N L × 1 vector defined as [A L ] is the usual CFIE system matrix of dimension N L × N L , and ] is evaluated as with i, j = 1, . . ., N L .We remark that, if it is needed to compare the error maps obtained for different excitations and geometries, the residual in (24) should be normalized by, for example, the norm of the corresponding RHS (26) The dimension N L reflects the number of levels L used to generate the reference mesh.This number tends to increase rapidly as each dyadic division of the triangle means multiplying the number of mesh cells by 4, leading to the computation of a large matrix.Although (23) does not need a matrix inversion, the matrix filling and the matrix-vector product could be made more efficiently via the HS-DDM technique presented in Section III-B, transforming a global problem into several local ones through the EHSs as detailed in the following.
We consider the initial division into P domains of the object under analysis (see Section III-A).Then, for each domain p, with p = 1, . . ., P, the radiation of the current densities on the object via ( 23) is separated into two components: one due to the local interaction of the current densities inside the domain and in its buffers regions (see Fig. 2), and the second due to the far-field current densities through the EHS.
The interaction within the domain and its buffer regions is obtained through the matrix [ A L p ], with each element equal to where f L p, j is the jth source basis function among N L p defined in the domain together with its buffers, while f L p,i is the ith test basis function among N L p defined in the domain only.Hence, the matrix [ A L p ] is not square as N L p > N L p .Instead, the interaction with the far-field region of the object is obtained via the matrix [Z S,L p ], where each element is equal to and i = 1, . . ., N L p .Finally, combining ( 27) and ( 28) with ( 23), we obtain where [B L p ] is the RHS on the domain p, being a subset of ], and [ I 0,L p ] collects the initial current density coefficients, projected in the reference level, for the basis functions defined on the domain and its buffers.To generate the error estimation associated with all the cells of the reference mesh, the described procedure has to be repeated for all the domains p, with p = 1, . . ., P.

D. Adaptive H-Refinement
The error estimation ( 19) is associated with each of the cells in the reference mesh.Then, as in [16], the adaptive h-refinement starts from the initial mesh through a dyadic subdivision.If the considered cell contains at least one reference mesh cell with an associated error above the chosen threshold, it is divided by four.If, instead, it does not contain any reference cell above the threshold, it is kept unmodified in the adaptive process.The dyadic subdivision is applied where required, recursively, up to the level L of the reference level.
Fig. 3(a) depicts an example of the initial mesh (thick lines), the corresponding reference mesh, and one triangle with an associated error above the threshold (highlighted in red).The adaptive h-refinement process will ensure the presence of that triangle above the threshold in the final h-refined mesh.It keeps the shared edges of the initial discretization, generating a nonconformal mesh due to the localized refinement, where one cell is divided but not its neighbor, as shown in Fig. 3(b).
The MB-RWG basis functions [21] are defined on a triangle, called positive, and a group of smaller triangles, called negative, with an edge in common with the positive triangle, where the edges from the negative side form exactly the edge in the positive side.By definition, the MB-RWG basis functions guarantee normal current continuity across common edges, meaning that no line charges exist and the functions are div-conforming.The nature of the dyadic h-refinement in Fig. 3(b) makes the MB-RWG basis functions very suitable for their use in the proposed approach.Fig. 3(c) shows the vectorial behavior of an MB-RWG basis function, in a mesh resulting from the performed adaptive h-refinement.Hence, the MB-RWG basis functions are defined, together with the RWG basis functions, on the obtained h-refined mesh.This choice avoids DG, which would imply implementing an extra penalty condition in the kernel of the SIE, as done in the formulation proposed in, for example, [18] and [19].

E. Adaptive H-Refined Solution
Once each domain p reaches an adaptive h-refined mesh, the complete problem can be solved again but, unlike [16], here the solution is efficiently achieved using the HS-DDM scheme.In the following, we will use the subscript H to indicate the adaptive h-refined mesh level.Similar to the error estimation part explained in Section III-C, the global coupling in the adaptive h-refined solution is separated into different parts.In this case, we find the unknown current density in each domain p, considering three different contributions: 1) the h-refined extended domain, defined by the domain itself Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
and its inner buffer, interacting with itself; 2) the outer buffer current densities radiating to the extended h-refined domain; and 3) the far-field region radiating through the EHS to the extended h-refined domain.
The interaction between the basis functions inside the domain and its inner buffer is represented by the matrix where f H p is a basis function defined in the extended h-refined domain, and i, j = 1, . . ., N H p .The radiation from the outer buffer to the h-refined domain of interest is done through the matrix [ ǍH p ].This matrix is rectangular as it goes from the Ň 0 p basis functions of the initial mesh in the outer buffer to the N H p basis functions in the h-refined extended domain.A generic element of the matrix is equal to where f0 p, j is the jth basis function defined in the outer buffer discretized via the initial mesh, with j = 1, . . ., Ň 0 p .The test basis functions are the same as in (30).
The last contribution is that from the far-field region and is made through the EHS.This is represented by the matrix [Z S,H p ], where the current densities on the EHS, [I S p ], radiate to the h-refined extended domain p.This matrix is of the same nature as the one in (28), but changing the testing functions.The matrix [Z S,H p ] dimensions are N H p , that is, the number of testing functions in the h-refined extended domain, by 2K p , that is, two times the basis functions on the EHS, as in (28).Then, each element of [Z S,H p ] is equal to where i = 1, . . ., N H p , and f H p,i is the ith testing function in the h-refined extended domain.
Collecting (30)-(32), we can set for each domain p the matrix system where [ B H p ] is the RHS vector in the h-refined extended domain obtained as in (26).The array [ Ǐ 0 p ] collects the coefficients of the initial current density in the outer buffer zone.The unknown in (33) is the coefficients array [ I H p ] representing the final current density in the adaptive h-refined extended domain.Hence, (33) can be written as where all the terms on the RHS are known and represent the excitation of the h-refined extended domain p.The solution [ I H p ] can be obtained through a direct inversion of the matrix, or via a fast MoM solver, such as the MLFMA, if the dimensions of the local problem require it.The system (33) is defined on the extended domain, hence, after obtaining the solution, we discard the coefficients in the inner buffer (that will be obtained via another domain), keeping those in the domain itself, only.This solution procedure is then repeated for each domain.

F. Workflow Summary
A summary of the different steps needed to obtain an adaptive h-refined solution is summarized here, following Fig. 4.
We start from an initial coarse mesh, where we obtain an initial solution and define the domains that will be used in the following steps [see Fig. 4(a)] as detailed in Section III-A.For each domain, an EHS is defined [see Fig. 4(b)] and its equivalent current densities are obtained (see Section III-B), as well as the localized error estimation with the definition of the reference mesh (see Section III-C).Via the error estimation, the adaptive h-refined mesh is generated (see Section III-D), selecting the reference mesh cells above the chosen error threshold [see Fig. 4(c)].
Fig. 4(d) depicts the collection of all the domains after the adaptive h-refinement, where a solution will be obtained.In Fig. 4(e), the HS-DDM is applied to each domain to obtain the adaptive h-refined solution (see Section III-E), shown in Fig. 4(f).

IV. NUMERICAL RESULTS
In this section, we illustrate the effectiveness of the proposed HS-DDM implementation in solving two different EM examples with real-life interest, considering both a scattering and a radiation problem.
All simulations shown in this article were performed in a cluster with 2 × AMD EPYC 7H12@2.6GHz (2 × 64 cores = 128 cores) and 2 TB of RAM memory.

A. Scattering Problem: Rafale Aircraft
A first numerical example is introduced to examine the efficiency of the proposed method to obtain an accurate solution for an adaptive h-refined mesh.A morphed version of a Rafale aircraft, shown in Fig. 5, is considered as the structure to analyze.This example has small and thin regions around the wings and the nose and some sharp edges around all the fuselage.The excitation of the problem is an oblique (θ = 30 • , φ = 30 • ) impinging plane wave at 137 MHz.The dimensions of the aircraft are approximately 15.3 m in length, 4.6 m in wingspan, and 3.9 m in height (7λ × 4.6λ × 1.8λ at the working wavelength, λ).
We start obtaining the solution of an initial mesh of the problem through a fast MoM method, that was modeled using 7986 unknowns corresponding to RWG basis functions.For this example, since the geometry does not exhibit multiscale features, we use an MLFMA-FFT [3] method to obtain the initial solution.The domain decomposition of the geometry is  then applied, splitting the Rafale aircraft into eight domains, shown with different colors in Fig. 5 together with the initial mesh.Next, each domain is encapsulated in an EHS, and we evaluate the error of the initial solution with respect to a reference one, obtained with a mesh that has two levels of cells dyadic subdivision (L = 2) from the initial one, as described in Section III-C; the reference mesh is formed by 127 776 RWG basis functions.The obtained error estimation is shown in Fig. 6 that represents the error values assigned to each triangle of the reference mesh.In the resulting error map, we can identify the regions with a high estimated error, concentrated mainly in the sharp parts (edges and borders).Then, the automatic h-refinement procedure, described in Section III-D, is applied to the initial mesh, comparing the error map with an error threshold equal to −3.5 (in logarithmic scale) [16] and obtaining an adaptive h-refined mesh, shown in Fig. 7, that provides a total of 54 405 unknowns (corresponding to 50 073 RWG and 4332 MB-RWG basis functions).
Finally, to find the solution in the adaptive h-refined mesh, we apply the proposed HS-DDM scheme, as described in Section III-E; in the following, we will refer to this approach as "adaptive h-refined (local)."To characterize and compare the performance of the HS-DDM method, also an MLFMA-FFT Fig. 7. Adaptive h-refined mesh of the Rafale aircraft.In gray, the plane region where the near field is evaluated.

TABLE I SIMULATION TIMES (IN SECONDS) FOR THE RAFALE EXAMPLE
solver is applied to the whole structure to obtain the solution of the adaptive h-refined mesh; the results obtained with the MLFMA-FFT solver are marked in the following with "adaptive h-refined (global)." Fig. 8(a) shows the convergence in terms of the iteration number under an iterative Krylov resolution of the matrix system (we are using GMRES [26]), comparing different solutions to better understand the improvement provided by the proposed approach.In Fig. 8, the green line corresponds to the solution with the initial mesh using the MLFMA-FFT, the black line shows the convergence of the MLFMA-FFT using the reference mesh, the red line represents the convergence of the adaptive h-refined mesh applying the MLFMA-FFT to the whole structure, while blue lines illustrate the convergence of the adaptive h-refined mesh problem using the HS-DDM method, where each blue line refers to the solution of one domain (obtained via a local MLFMA-FFT).We can observe that the convergence of all domains is faster than the rest of the solutions.Table I lists the times needed for the setup and solution of the initial problem (using the initial mesh), the time to generate the error estimation map and the currents on the EHSs, and the setup and solution times for the solution using the final refined meshes.The time for the local refinement is not reported because it is negligible.In columns two and three, Table I reports the required times using the adaptive h-refined global and local approaches, respectively; it is evident that the proposed local approach (i.e., the HS-DDM) outperforms the global one.The last column reports the setup and solution times using the reference mesh.In this case, the reference solution requires higher times than the adaptive h-refined solutions (both global and local).
However, considering that the different convergences shown in Fig. 8(a) corresponds to problems with a different number of unknowns, the number of iterations cannot be reliable to compare the performance of the different results, such as the times per iteration are different.A better figure of merit is shown in Fig. 8(b), where the above iterative convergences are expressed in terms of the total time to solve each problem.In Fig. 8(b), comparing the convergence of the proposed HS-DDM scheme (blue lines) with the MLFMA-FFT in the case of the adaptive h-refined problem (red line), it can be observed a speedup factor greater than ×7 applying a parallel implementation of the HS-DDM, which can be easily integrated in the method such as the solution of one domain is totally independent of the rest of domains.But, even without applying a parallel scheme in the solution of the different  domains, the HS-DDM solver is more than ×2 faster than the MLFMA-FFT one.Finally, comparing the adaptive h-refined solutions with respect to the reference one, there is a speedup factor of ×15 using the HS-DDM (parallelizing the domains solutions) and ×9 using the MLFMM-FFT (parallelized with eight processors); in both cases (HS-DDM and MLFMM-FFT), we consider also the time for the generation of the initial solution that is present in both approaches.
To verify the solution accuracy of the different methods, first, we examine the radiated electric near-field on a plane at λ/20 to the upper wing, as shown in Fig. 7.The obtained near fields are reported in Fig. 9  evaluated with the initial mesh is quite different with respect to the reference one with a maximum error of −6.76 dB.
These findings are coherent with the obtained electric current density over the aircraft surface, shown in Fig. 10, where there is a good agreement between the adaptive h-refined mesh solutions [see Fig.

B. Antenna Problem: Modern Vessel
As a second example, the radiation of a realistic and multiscale structure is considered to evaluate the accuracy and efficiency of the proposed method.The structure is a modern vessel with four integrated patch antennas embedded into the main mast (see Fig. 11).The approximate dimensions of the ship are 140 m in length, 20 m in width, and 40 m in height (23.3 λ × 3.3 λ × 6.7 λ in the used working wavelength corresponding to 50 MHz).These antennas are designed to work together as an omnidirectional VHF communication system.The excitation is a delta gap in the feeding points of the antennas.
In this case, the geometry of the problem exhibits very deep multiscale features and a localized excitation, so, instead of the MLFMA-FFT solver (that has a slow convergence), the DDM of [6] is applied to obtain the solution of the initial mesh, the global solution of the adaptive h-refined mesh and the one of the reference mesh.For this structure, the initial mesh is composed of 76 578 RWG basis functions (see Fig. 12), and the reference one by 306 900 unknowns, that correspond to one level of dyadic cells subdivisions (L = 1) from the initial mesh.The entire problem is partitioned into 14 domains, shown in Fig. 11: ten large domains belonging to the vessel superstructure, and four electrically small domains corresponding to the four patch antennas.

TABLE II SIMULATION TIMES (IN SECONDS) FOR THE VESSEL EXAMPLE
After obtaining the electric current density on the surfaces of the initial mesh, the error estimator algorithm is applied using the proposed local scheme to generate the error map shown in Fig. 13.We can see that, in this case, the maxima of the error are concentrated in the proximity of the antennas and the illuminated parts of the structure.Then, after comparing that error map to an error threshold of −3.5 (in logarithmic scale), the technique is applied, resulting in the adaptive h-refined mesh shown in Fig. 14, providing a total of 17 7752 unknowns (176 010 RWG and 1742 MB-RWG basis functions).We can observe that the obtained adaptive h-refined mesh concentrates the refinement close to the radiating systems; the vessel parts, far from the antennas, maintain the initial mesh (except on the edges), while the discretization of the masts and the antennas is refined to improve the accuracy of the solution.15, where, for the HS-DDM method (blue line), the solution of the domain with the worse convergence in terms of time is represented only, to avoid the superposition of too many lines.The performance improvement in the solution of the adaptive h-refined problem via the HS-DDM (blue line) is also evident for this radiation case, with a speedup factor of ×10 with respect to its solution via the global DDM (red line) [6].As in the previous case, Table II summarizes, for the vessel case, the required times for the initial setup and solution, for the generation of the error map and the currents on the EHSs, the final mesh setup  and solutions, and the reference mesh setup and solution times.In this case, the proposed local approach (i.e., the HS-DDM) outperforms the global one and there is a gain in time with respect to the reference solution (last column).It is important to notice that the reference mesh, and consequently the adaptive h-refined mesh, has only one level of refinement below the initial coarse mesh, and, in the adaptive h-refined mesh, the refinement is concentrated mainly on the mast.As in the previous example, we evaluate the radiated electric near-fields on a plane at λ/10 from one of the antennas, shown in Fig. 12. Fig. 16(a)-(d) reports the near fields obtained from the solution of the initial mesh, the reference one, and the adaptive h-refined mesh, both with the DDM and the HS-DDM one.We can observe a good agreement between the reference field and the ones obtained with the adaptive h-refined mesh, with a maximum error of −34.89 dB for the global DDM and −27.71 dB for the HS-DDM, while with the initial mesh, the maximum error is equal to −10.75 dB.Finally, the electric current density is shown in Fig. 17 for all the approaches.To better show the differences among the solutions, in Fig. 18, there is a detailed view close to one excited antenna.In the evaluation of the current density, the maximum error of the adaptive h-refined solutions, with respect to the reference one, is equal to −33.22 dB for the global DDM and −27.82 dB for the HS-DDM, while it is −14.62 dB using the initial mesh.

V. CONCLUSION AND PERSPECTIVES
In this work, we presented a new HS-DDM for obtaining an efficient and accurate adaptive h-refined problem solution for multiscale realistic PEC structures.The error estimation and solution of the adaptive h-refined problem are, in general, globally coupled, requiring the handling of large matrices.In this work, we proposed dividing the structure into domains and considering the effects of external (far-field) regions through equivalent surfaces at half-wavelength from those regions and based on Huygens' equivalence principle, going toward a localized problem.The initial solution, obtained via an average λ/10 mesh discretization, is the seed for all the proposed technique steps and defines the equivalent currents that link the domains in all the processes.Error estimation and adaptive h-refined solution are performed locally through the use of Huygens' surfaces.Moreover, in the adaptive h-refined mesh, which can be nonconformal, we used MB-RWG and RWG basis functions, avoiding the extra interior penalties of the DG MoM needed, instead, if half-RWG basis functions are applied.As the numerical results showed, the HS-DDM allows a significant speed-up in the solution process, keeping a good accuracy in the evaluated current density and near-fields.
Two possible future developments of the approach presented here are the extension to objects composed of multiple materials, and the implementation of the multiresolution preconditioner [27], [28] for the solution of the local domains to further improve the overall solution time.Moreover, the number levels in the reference mesh could be different for each region or triangle to better adapt the scheme to the initial mesh size and the problem excitation.

Fig. 4 .
Fig. 4. Workflow of the HS-DDM procedure to generate an adaptive h-refined mesh and to find the corresponding current densities.(a) Initial mesh and domain decomposition.(b) Initial mesh domain and its equivalent Huygens' surface (EHS).(c) Refined mesh for a domain.(d) Total refined mesh and domain decomposition.(e) Refined mesh domain and its equivalent Huygens' surface (EHS).(f) Solution current.

Fig. 5 .
Fig. 5. Partition into domains and initial mesh of the Rafale aircraft.

Fig. 6 .
Fig. 6.Rafale aircraft error map associated with the initial mesh.

Fig. 8 .
Fig. 8. Rafale aircraft: iterative convergence in terms of (a) number of iterations, and (b) solution time.

Fig. 9 .
Fig. 9. (a) Reference radiated electric near field.(b) Difference between initial and reference electric near fields.(c) Difference between adaptive h-refined mesh (global solution) and reference electric near fields.(d) Difference between adaptive h-refined mesh (local solution) and reference electric near fields.All the reported data are in dBµV/m.The plane where the near fields are evaluated is shown in Fig. 7 (inset).

Fig. 10 .
Fig. 10.Rafale aircraft surface current density and differences (dBµA/m).(a) Reference current density, (b) difference between initial and reference current densities, (c) difference between adaptive h-refined mesh (global solution) and reference current densities, and (d) difference between adaptive h-refined mesh (local solution) and reference current densities.
(a)-(d) for the solution of the initial mesh, the reference mesh, and the adaptive h-refined mesh, both with the global MLFMA-FFT solver and the HS-DDM one.The near-field obtained with the adaptive h-refined mesh matches very well with the reference one, with a maximum error of −25.96 dB for the global MLFMA-FFT case and −23.65 dB for the HS-DDM one.Instead, the near-field

Fig. 11 .
Fig. 11.Geometry and partition into domains of the vessel.

10
(c) and (d)], with respect to the reference one [see Fig. 10(b)], with maximum errors equal to −18.51 dB for the MLFMA-FFT solver and −16.82 dB for the DDM-HS one, while the initial mesh current density is quite different [see Fig. 10(a)], with a maximum error of −7.36 dB.

Fig. 12 .
Fig. 12. Vessel initial mesh.Plane where the near field is evaluated.

Fig. 13 .
Fig. 13.Vessel error map associated with the initial mesh.

Fig. 15 .
Fig. 15.Iterative convergence in terms of the solution time for the vessel.

Fig. 16 .
Fig. 16.(a) Reference radiated electric near-field.(b) Difference between initial and reference electric near-fields.(c) Difference between adaptive h-refined mesh (global solution) and reference electric near-fields.(d) Difference between adaptive h-refined mesh (local solution) and reference electric near-fields.All the reported data are in dBµV/m.The plane where the near-fields are evaluated is shown in Fig. 12 (inset).

Fig. 18 .
Fig. 18.Detailed view of the vessel surface current density and differences (dBµA/m).(a) Reference current density.(b) Difference between initial and reference current densities.(c) Difference between adaptive h-refined mesh (global solution) and reference current densities.(d) Difference between adaptive h-refined mesh (local solution) and reference current densities.