3-D Electromagnetic Scattering Computation in Free-Space With the FETI-FDP2 Method

The electromagnetic dual-primal finite element tearing and interconnecting (FETI-DPEM) method is a nonoverlapping domain decomposition method developed for the finite element analysis of large-scale electromagnetic problems, where the corner edges are globally numbered. This paper presents an extension of the FETI-DPEM2 method, named FETI-full dual primal (FETI-FDP2), where more flexible Robin-type boundary conditions are imposed, on the inner interfaces between subdomains as well as on the corner edges, leading to a new interface problem. Its capacities are tested in the framework of a three-dimensional (3-D) free-space scattering problem, with a scattered field formulation and a computational domain truncated by perfectly mathed layers (PML). First, we compare its accuracy with respect to other FETI-DPEM2 methods and to a complete resolution of the FEM problem, thanks to a direct sparse solver. We show that the convergence of iterative solvers is affected by the presence of the PML and can be accelerated by means of a more accurate approximation, between adjacent subdomains, of the Dirichlet-to-Neumann (DtN) operator. The effectiveness of the iterative solvers are also considered for different test cases. The advantages of the proposed FETI-FDP2 method combined with the associated DtN approximation is numerically demonstrated, regardless the chosen working frequency or the iterative solvers.


I. INTRODUCTION
The need of engineering simulations of large and complex structures is rapidly growing, requiring numerical methods which are more and more efficient. As such, the Finite Element Method (FEM) applied to the resolution of time harmonic electromagnetic wave scattering has become very popular over the past decades. Indeed, the resolution of the weak form of the Helmholtz equation has shown its potential and its versatility among different configurations and types of problems to be solved (complex media, periodic structures, anisotropic materials . . . [1]- [3]).
In this method, the unknown of interest (the electric or magnetic field) is expanded onto a set of basis functions. A linear system is then calculated by projecting the Helmholtz equation onto the same set of test functions, as advocated by the Galerkin method. Thus, the efficiency of the method is mainly dependent upon one's own ability to solve the resulting sparse linear system, which can be time and memory consuming, especially in three-dimensional (3D) configurations.
Among the different schemes proposed to solve large scale models and preserve the versatility of the FEM method, the Domain Decomposition Method (DDM) and its different evolutions are especially appealing [4], [5]. The closely related Finite Element Tearing and Interconnecting (FETI) method seems also very robust when one is dealing with arbitrarily mesh partitions. The general principle of FETI methods is first to divide the entire computational domain into smaller non-overlapping subdomains. In each of these subdomains, local linear systems can now more easily be inverted. Simultaneously, the different adjacent subdomains are glued at their common interfaces thanks to appropriate boundary conditions and constraints imposed to the field and its derivatives, leading to a so-called global interface problem. Once this interface problem is solved, the solution inside each subdomain can be evaluated independently by using the known mixed boundary conditions at the internal interfaces between subdomains. This method has been applied in many domains going from mechanics [6], [7], to acoustic wave propagation [8]- [10], and electromagnetism [11]- [17].
In order to improve the convergence and the scalability of the method, one can notice the existence of two techniques. The first one uses the plane wave spectrum operator [8], while the second one uses dual-primal techniques, denoted as FETI-DPEM when applied to electromagnetic cases, which can be seen as coarse grid corrections [9], [12], [14], [18]. Over the past years, the latter has been successfully applied to the simulation of large-scale electromagnetic problems [11], [19]. In this last method, the corner nodes in 2D or the corner edges in 3D (the "corners" corresponding to geometrical entities belonging to more than two adjacent subdomains) are extracted from each subdomain and are globally and uniquely numbered, in combination with the boundary conditions imposed at these corners.
The various boundary conditions can be easily imposed by means of appropriate Lagrange multipliers. In FETI-DPEM, a single Lagrange multiplier is introduced in order to impose a Neumann boundary conditions at each internal interface, as well as at the corner points. In FETI-DPEM2, each internal interface apart from the corners is equipped with two Lagrange multipliers, which is equivalent to imposing a Robin type boundary condition, thus avoiding the appearing of spurious solutions. This boundary condition can be also seen as a crude approximation of a transparency operator and many efforts have also been performed in order to optimize the coefficients arising in this boundary condition for plane interfaces between subdomains [10], [16], [17], [20].
In a previous work, performed in 2D [21], we have extended the FETI-DPEM2 method by imposing everywhere, even on the corner nodes, more flexible Robin-type boundary conditions. The resulting interface problem is completely rewritten and acts as a new global coarse preconditioner. In 2D, a series of numerical tests have shown the efficiency of this modified method.
In this paper, we follow the same ideas applied this time to a 3D configuration and for the vectorial Helmholtz equation. Moreover, we analyse a 3D free-space scattering problem, which leads us to implement a scattered field formulation within a domain truncated by Perfectly Matched Layers (PML) [22], [23]. After developing the linear system associated to the internal interfaces problems due to the different variations of the FETI-DPEM2 methods and to our proposed FETI-FDP2 method, we test the accuracy of each method with respect to a direct FEM solution computed thanks to a direct sparse solver. Then we focus our attention on the convergence behaviour of the iterative algorithm which must be employed to solve the interface problem. Indeed, iterative procedures appear unavoidable when dealing with large size problems. We show in this paper that the convergence of the iterative algorithm is strongly affected by the presence of anisotropic materials. We also show how the method can be strongly accelerated by partly adopting the Evanescent Modes Damping Algorithm (EMDA) [24] when PML are present. This EMDA method originally was proposed as a more accurate approximation of the Dirichlet-to-Neumann (DtN) operator for the nonoverlapping Schwarz DDM [24], instead of classical Sommerfield boundary conditions between adjacent domains. One can notice that second order transmission boundary conditions or high order transmission boundary conditions [10], [25]- [27] have been developped and have shown their efficiency, though we have focused our study on the implementation and optimization of the EMDA approach in the framework of our proposed method since it avoids the appearing of supplemental terms in the discrete finite element linear system. The paper is organized as follows. Section 2 describes the considered 3D free-space electromagnetic scattering configuration and the associated scattered field formulation. In Section 3, following the notations introduced in [13] [21], the mathematical formalism of the different FETI-DPEM2 and FETI-FDP2 methods are described. The constructions of the full and reduced linear systems obtained for the various interfaces unknowns, either the corner edges or the inner interface edges, are detailed. All the numerical results are gathered in Section 4. At first, in Section 4.A, a small 3D scattering problem is constructed. It enables to test, in Section 4.B, the accuracy of all the presented methods with respect to a classical resolution of the same FEM problem. In Section 4.C, we investigate the convergence behaviour of a robust Generalized Minimal Residual Method (GMRES) [28] iterative algorithm for solving the interface problem. In particular, the influence of the different methods as well as the presence of PML is analysed. The convergence improvement obtained thanks to a Mixed approach which exploits in part the EMDA method is also numerically shown. A similar analysis is performed in Section 4.D for a second type of iterative solver, that is the BiConjugate Gradient Stabilized method (BiCGStab) [28]. Finally, the behaviour of all these methods are presented in Section 4.E for large-scale problems. Conclusions follow in Section 5.

II. SCATTERED FIELD FORMALISM
We consider a three-dimensional electromagnetic scattering problem, where a known monochromatic incident electromagnetic wave is impinging on an inhomogeneous target D ⊂ Ω, whose relative permittivity varies with respect to the surrounding. The relative permittivity and permeability distribution in absence (resp. in presence) of the scatterer are denoted by ε inc r ( r) and µ inc r ( r) (resp. ε tot r ( r) and µ tot r ( r)) with From the linearity of Maxwell's equations, the total field (associated to the permittivity distribution ε tot r ) can be decomposed into the incident field E inc (associated to ε inc r ) and the scattered field E sc which satisfies a Helmholtz equation where k 0 is the vacuum wavenumber. The induced currents are defined as (3) The scattered field also satisfies a radiation condition (RBC) The computational domain can be bounded thanks to a Perfectly Matched Layer (PML) [22], [23] which can be interpreted as an anistropic absorber [29] where the tensors µ pml r and ε pml r are defined with respect to µ tot r and ε tot r along with stretching factors. This type of absorbing layer helps to prevent spurious reflections from the external boundaries which may arise when the radiation boundary condition is simply set.
When the sources and the receivers are not necessarily located in the vicinity of the target, the meshing of the model space is too memory-consuming. Indeed, in finite-element methods, for first-order vector basis elements as it is the case herein, the number of unknowns corresponds to the global number of edges in the mesh. In that case, the computation domain can be restricted to the scatterer's area and a Near-to-Far-Field transformation based on Huygen's principle [2] can be employed to extract the fields at the receivers locations.
Nonetheless, when the scatterer's external dimensions are large with respect to the wavelength, the size of the model domain, even if it is moulded to the shape of the target, requires heavy computational resources. It is therefore of interest to derive ad-hoc and optimized numerical schemes for solving the aforementioned Helmholtz equation.

METHOD
Among the available domain decomposition-based methods, we have adopted the FETI method, and in particular the FETI-DPEM2 versions [13], [14] to solve our scattered field problem. While recalling the various variations of the FETI-DPEM2 method with the same notations as in [13], we will detail how we have extended it in order to impose everywhere more flexible Robin-type boundary conditions. This new approach will be denoted as the FETI-FDP2 (full dual-primal) method.

A. Geometrical feature extractions
According to the Domain Decomposition idea, we assume that the domain Ω (with external boundary ∂Ω) is divided into a set of N s non-overlapping subdomains Ω = For a given subdomain Ω i , the indices of all its adjacent subdomains are called neighbor(i). The internal boundary of the subdomain Ω i is denoted as Γ i . The internal boundary between two subdomains Ω i and Ω j , with j ∈ neighbor(i), is denoted as Γ ij .
Following the idea of the FETI aproach [8], we construct at first an augmented Lagrangian functional. Its associated Karush-Kuhn-Tucker conditions [11] lead to the following equation in each i-th subdomain K i is the sum of the stiffness matrix, the mass matrix and the external and internal boundary condition matrices where N i denotes a column vector containing the curlconforming vector basis function in the i-th subdomain. The role of α i is explained in Section III-C. E i is the unknown components of the vector field E sc in Ω i and f i represents the discretized version of the right-hand-side of (Eq. 2). D i is a Boolean matrix which extracts only the edges of Ω i which are on Γ i . As for the dual Lagrange multipliers λ i , they are also involved within the unknown boundary conditions which are imposed at the internal interfaces between adjacent subdomains.

B. Ordering with global corner edges
In each subdomain, we split the field unknowns in the following way where the notations E i r denote all the internal (E i V ) and interface (E i I ) Degrees of Freedom (DOFs) belonging to the ith subdomain except for the corner DOFs which are denoted by E i c . The main feature of the FETI-DPEM methods is to consider the corner edges in a different way from the other edges. Indeed, these corner edges are merged, then globally and uniquely numbered to provide a common vector E c . With such notations, we can now rewrite (Eq. 5) as follows We now introduce several Boolean matrices to gather the DOFs of interest. The matrix D i r selects from all the DOFs of Ω i the ones which are located on its interface, i.e., The projection matrix T i→j r extracts the DOFs associated to the interface Γ ij from the ones belonging to Γ i , A similar Boolean matrix T i→j c is here to select all the corner DOFs belonging to Γ i and associated to Γ ij , Two supplemental Boolean matrices Q i λr and Q i Ec are introduced such that From the first line of (Eq. 7), the electric field in the i-th subdomain can be easily derived From the second line of (Eq. 7), the corner DOFS can be estimated By assembling and summing all the subdomains contributions, the global corner points related system equation is obtained as follows with where the operator S Ec is defined by

C. Boundary conditions at internal interfaces
A classical continuity boundary condition is first imposed on the tangential component of the electric field where i↔j corresponds to the jump boundary condition on Γ ij . This equation is rewritten as where E i→j (resp. E j→i ) is the subset of E i associated to the edges of Ω i (resp. Ω j ) which belong to the interface Γ ij . The introduction of dual Lagrangian multipliers enables to insert, at the internal interface Γ i , either an unknown Neumann-type boundary condition (when α = 0) or an unknown Robin-type boundary condition (when α = 0) It leads to the following set of equations, taking advantages of (Eq. 17), where the matrix W i↔j = M i→j + M j→i directly translates the boundary conditions. Following the aforementioned "rc"notations, this equation can be rewritten as follows: The FETI-DPEM2 and FETI-FDP2 methods differ in the way these boundary conditions are effectively introduced at the internal interfaces, and more particularly at the corner edges.
D. Full Interface Problem construction 1) FETI-DPEM2 (a) method: In [13], a Robin-type boundary condition is imposed for each internal interface, apart from the corner edges. Instead, for the DOFs related to the corners, the Neumann-type boundary condition is applied. Moreover, during the "gluing" process, no interactions are introduced between the corner and interface DOFs. In that case, (Eq. 20) is simplified into Thanks to (Eq. 12), we can eliminate E j→i r from the first line of (Eq. 21). By summing over all the subdomains, we obtain the first part of the FETI-DPEM2 (a) Full Interface Problem: where and the operator S λr is defined by Thanks to the definition of the operator S Ec , the second line of (Eq. 21) provides We can now assemble (Eq. 14), (Eq. 22) and (Eq. 24) to derive the FETI-DPEM2 (a) Full Interface Problem 2) FETI-DPEM2 (b) method: In a previous work, performed under a two-dimensional configuration [21], we showed that in order to limit computational errors, the connection between the corner and interface DOFs must be kept. In other words, the term W i↔j rc should not be neglected during the "gluing" process. This idea has also been adopted in the FETI-DPEM2 (b) method [26] where they transform the boundary conditions into With such a system, at the corner edges, a Neumann boundary condition is still applied while a Robin boundary condition is applied only for the rest of the interface DOFS. Taking this into account, the resulting equation (Eq. 22) is still valid [21], [26]. The only difference between the FETI-DPEM2 (a) and FETI-DPEM2 (b) methods relies on the definition of the matrix F λrEc , which now corresponds to 3) FETI-FDP2 method: In the framework of a twodimensional configuration [21], we not only showed that the term W i↔j rc should not be neglected, but we also showed that the method can be further enhanced by imposing a Robin-type boundary condition everywhere. In that case, (Eq. 20) should be kept as is.
Nevertheless, if we try to keep all the Lagrange multipliers λ c related to the corner points, we are facing singularity issues. To avoid it, in our previous 2D work [21], we suggested to look for a special set of Lagrange multipliers, such that which results in a matrix F λcλc = 2I.
In the current 3D configuration, we propose a more general framework. Instead of searching for λ i→j A new projection Boolean matrix Q i↔j λc is introduced to extract the double Lagrange multipliers λ i↔j c associated to Γ ij from the full list of λ c , i.e., λ i↔j c = Q i↔j λc λ c . A new operator S λc is also defined such that [11] By applying this new operator to the second line of (Eq. 28), we obtain where The matrix F Ecλc can not be introduced through the operator S Ec . It is a Boolean projection matrix which puts in correspondence all the global corner Lagrange multipliers λ i↔j c related to the global corner E k c for every k = 1...N Ec . The size of this matrix is equal to N Ec × N λc and its construction is quite simple. In each line (that corresponds to the global index of the corner edge) we need to put 1 if the given global corner Lagrange multiplier belongs to this corner and 0 if not.

E. Reduced Interface Problem construction
In practice, the Full Interface problem is not directly solved but replaced by a linear system with smaller dimension and better conditioning number [30]. Indeed, the Lagrange multipliers λ c are of not great interest and can be directly eliminated from the previous set of linear equations. Once these Lagrange multipliers are dropped out after few mathematical transformations, we arrive at the Reduced Interface Problem whereF We can now easily compute E c by solving the following system of equationŝ Finally, replacing E c and λ r in (Eq. 12) gives access to E i r in each subdomain Ω i . Note that, in this work, we only consider the Reduced Interface Problem.

IV. NUMERICAL RESULTS
We now analyze the numerical behavior of the proposed approaches in the framework of a 3D free-space scattering configuration.
The factorization of each sub-matrices is performed thanks to a direct sparse solver [31] and stored during the resolution of the interface problem. The factorisation of the matrixF EcEc is performed using Intel Math Kernel Library and kept during all the iterative processes. All computations have been performed on an Intel(R) Xeon(R) CPU X5570 @ 2.93GHz, with 48 GB of RAM, with no parallel programming specificities. tivities ε 1 r , ε 2 r , ε 3 r and ε 4 r are respectively equal to 1.5, 2.0, 2.5 and 3.0. A plane wave illuminates the set of scatterers, with an incidence angle of 45 • and a linear polarization along E θ .

A. Configuration description
The computational domain Ω is a parallelepiped whose size depends on the wavelength. For example, at 4GHz, the investigation domain Ω is a [0.32 × 0.34 × 0.22] m 3 box (≈ 61λ 3 , where λ is the background wavelength). The domain is discretized thanks to a free unstructured mesh generator Gmsh [32]. This global mesh is inserted into a home-made mesh partitioner which provides, for each subdomain, the structure of the local mesh as well as the boundaries and the corner edges lists between subdomains. The partitioning is definitely irregular as shown in Figure 2. This partitioner intensively uses subroutines provided by Metis [33].

B. Direct frontal solver
At first, the computational domain domain is surrounded with PML, with a width w PML of one wavelength. The scattering problem is solved thanks to a classical Finite Element Method without any partitioning, thanks to the frontal solver Mumps [31] based on a LU decomposition. This provides the scattered field E sc FEM . The same frontal solver is used to solve the Reduced Interface Problems associated to each of the three aforementionned methods, with α = jk 0 . A L 2 -norm error is introduced to compare the various scattered fields and the associated errors are provided in Table I. It is obvious that the relative error obtained with the FETI-DPEM2 (a) method is significantly higher than the other ones. This comforts us in the fact that, as proposed in [21], the term W i↔j rc of (Eq. 20) cannot be neglected. From now on, we will thus stop considering the FETI-DPEM2 (a) method. The second obvious conclusion is that the FETI-DPEM2 (b) and FETI-FDP2 scattered fields are identical (up to the computer precision) to the field computed with the classical FEM method. This fully validates the proposed approaches and their implementations.

C. GMRES iterative solver
As direct solvers require a large memory storage and as we want to tackle large-scale electromagnetic scattering problems, we investigate the influence of iterative solvers on the various methods. At first, we employ a robust iterative solver, a crude Generalized Minimal Residual Method (GMRES) [28], [34]. We set m = 15 for the size of the Krylov subspace and 10 −6 for the stopping criterion. It is well-known that methods based on Krylov-subspaces, such as GMRES, are well adapted for solving linear systems, even if they might poorly converge when the systems are not necessarily positive definite. Unfortunately, neither the Full Interface Problem nor the Reduced Interface Problem are positive definite systems. As finding a good preconditioned for solving a given sparse linear system is often viewed as a combination of art and science, we propose to work without any preconditioner and to focus on the various transmission boundary conditions influences.
1) PML influence: We divide the investigation domain into 42 subdomains and follow the convergence process of the iterative GMRES solver in presence of a λ-width PML layer. Indeed, several numerical tests have pointed out that the PML layers might have a drastic influence on the convergence behaviour of the iterative solver. Figure 3 shows the convergence history of the modulus of Lagrange multipliers associated to the subdomain Ω 15 . The choice of this specific subdomain is due to the fact that, among all of them, the relative L 2 -error in Ω 15 is the highest one. It turns out that this domain contains the biggest number of Lagrange multipliers located within a PML zone. As one can see from Figure 3, after 10 GMRES iterations, there is a high correlation between the numerical and exact solutions in the areas not associated to PML, while there is little correspondence between these values in the PML region.
2) RBC influence: Instead of surrounding the computational domain with PML, we employ a radiation boundary condition (RBC). The corresponding convergence results are presented in Figure 4, where the residue corresponds here to the error committed at each iteration of the iterative process between the left-hand-side and the right-hand-side of the Reduced Interface Problem. The significant difference in terms of number of iterations confirms the fact that the PML badly impact the convergence behavior, whatever the chosen FETI method. However, this is not the case when RBC are applied.
Indeed, as the Lagrange multipliers are approximated within the PML, once the local problems are solved, they are all affected by these errors. Thus, the misfit between the "true" FEM solution and the FETI solution is large. On the contrary, when there is only an absorbing boundary condition, the Lagrange multiplier error is small everywhere and does not affect as much the local problems. The misfit is thus reduced.
3) Internal boundary conditions influence: Internal boundary conditions are also known to impact the convergence of the domain decomposition algorithms. Indeed, when α is simply set to jk 0 as in the Despres approach [35], the approximation of the Dirichlet-to-Neumann operator of (Eq. 18) does not treat efficiently the evanescent modes [10]. Recent works are mainly devoted to the implementation of second or higherorder transmission conditions and have shown their efficiency [10], [25]- [27]. Here, we prefer to focus on simpler first-order transmission conditions and see how we can combine them in order to improve the convergence behaviour.
The role and the impact of α on the convergence of the domain decomposition algorithms has already been studied in the framework of DDM and FETI-methods but for plane interfaces [17] [15] [25] [36]. Since we are dealing with an arbitrary mesh partitioner yielding strongly irregular surfaces with varying curvatures, we have decided to only select the simplest choice for this parameter.
In particular, we exploit the Evanescent Modes Damping algorithm (EMDA) proposed in [24] which enables to extend the transmission boundary conditions to the evanescent modes. We only consider here the EMDA case where α = jk 0 (1+jχ), χ being a real-valued positive coefficient. The optimal value for χ depends in particular from the mean curvature on the interface [37]. As the partitioning is performed in an irregular fashion, we simply set χ = 0.5 in the following computations as it is favorized in [37].
We have played with the various types of transmission coefficients, and more specifically with their locations ( Figure 5). We thus compare: 1) the Despres approach, which states α = jk 0 everywhere, 2) the EMDA approach, which states α = jk 0 (1 + jχ) everywhere, 3) a new Mixed approach, with α = jk 0 (1 + jχ) in the PML areas and α = jk 0 everywhere else. The Reduced Interface Problem is first constructed with RBC and no PML. The Mixed (or equivalently Despres) approach and the EMDA approach are thus applicable. The associated convergence histories, obtained for the FETI-DPEM2 (b) and FETI-FDP2 methods are presented in Figure 6. It is clearly visible that the convergence is decreased by the use of the Mixed approach, whatever the employed FETI method.
The Reduced Interface Problem is now constructed with PML. The convergence results obtained with the FETI-FDP2 method are presented in Figure 7. The associated eigenspectra of the Reduced Interface Problems are also visible in Figure 8.  is equal to 501, 411 and 292 for the Despres, EMDA and Mixed approaches respectively. In presence of anisotropic media, the EMDA algorithm improves the conditioning of the Interface problem, and as a consequence, the convergence results. Nevertheless, it seems more interesting to follow the Mixed approach and to impose the EMDA only out of the PML zone in order to further improve the convergence process.  A similar case is considered, but this time with smaller PML (w PML = 0.4λ). The convergence behaviour of the GMRES method is presented in Figure 9. Several conclusions can be drawn. Firstly, the Mixed approach definitely improves the convergence process. Secondly, whatever the chosen transmission coefficients approximation, the FETI-FDP2 method converges as fast as the FETI-DPEM2 (b) method.

D. BiCGstab iterative solver
Various research works take advantage of other iterative methods, such as the BiCGStab method. We thus apply this iterative method (with a stopping criterion equal to 10 −3 ) to solve the Reduced Interface Problem in presence of PML. Unfortunately, for the three types of transmission conditions, this method only converges for one of them, that is the Mixed approach ( Figure 10). The Despres and EMDA approaches failed, whatever the employed FETI method.
The same problem is solved this time without PML but with RBC. The corresponding results are given in Table II. It clearly shows that the EMDA approach deteriorates the convergence and can even prevent the iterative solver from converging.
All these numerical tests lead to the conclusion that the Mixed approach should be favorized regardless the configuration and the choice of the iterative solver. Moreover, it seems that the combination of the FETI-FDP2 method and the Mixed  approach is more stable and should be promoted.

E. Large scale problems
As all these methods are meant for tackling large-scale electromagnetic problems, we increase the frequency to enlarge the number of unknowns. The stopping criterion for all the iterative solvers is set to 10 −3 , which is a relatively high number, but provides nevertheless a good indication of the behaviour of the algorithms. Figure 11 shows for example the repartition of the amplitude of the scattered field, when the frequency is set to 20 GHz. Table III summarizes the results obtained with the various methods for frequencies ranging from 4 GHz to 20 GHz, in presence of RBC and when the Mixed approach is used. In particular, in the last column, the factorization time of all the matrices K i rr is provided as well as the construction time and the factorization time required for the matrixF EcEc . In the FETI-DPEM2(b) method, only the computation of F EcEc is required. The addition of the new set of Lagrange multipliers λ c in the FETI-FDP2 method introduces 3 supplemental matrix-matrix multiplications, in order to obtainF λrEc ,F EcEc andd Ec which are fully performed at the pre-initialisation level. It is worth mentioning that the computational time associated to these operations is relatively small, as the biggest matrix F λcλr is very sparse and stored in a special format which represents a set of full matrices. For the largest problem considered in this work, this time varied from few dozens to few hundreds of seconds depending on the domain partitioning. The time required for one GMRES iteration (for m+1 multiplications of the Interface Problem) is also given and is the same for the FETI-DPEM2(b) and FETI-FDP2 method, as the factorization of the various matrices is performed beforehand and the number of unknowns is the same for the two Reduced Interface Problems. This table enables to appraise the computational burden associated to each of the proposed methods and their various combinations. When the frequency increases, the number of unknowns is getting too large and prevents us from inverting the linear system by means of a direct frontal solver. In those cases, the L 2error is thus not available. It is of interest to note that the L 2error, when computable, is of the same order as the stopping criterion associated to the iterative solvers. This indicates that the supplemental error due to the iterative algorithm must be added to classical numerical errors linked for example to mesh discretizations. For all the cases, the GMRES algorithm converges faster than the BiCGStab. Similarly, the FETI-FDP2 method converges as fast as or even faster than the FETI-FDPEM2 (b) method.

V. CONCLUSION
In this paper, a variation of a Domain-Decomposition method, named the FETI-FDP2 method, has been suggested. Thanks to Lagrange multipliers, this method applies Robintype boundary conditions everywhere, even at the corner points. The principle of this method has been detailed and its similarities and differences with the already existing FETI-DPEM2 methods has been provided.
With such a method, an effective code has been developed in order to solve 3D large-scale scattering electromagnetic problems which can contain inhomogeneities and anisotropy. This code is based on an arbitrary mesh partitioning, which is a complicate issue in 3D. Special care has thus been taken to handle internal interfaces which are not necessarily planar without introducing a major drawback in terms of computational accuracy. The accuracy as well as the implementation of this new method has been fully checked by comparing its results with the ones obtained with classical FEM and FETI methods, where direct solvers have been used.
Due to the increasing size of the underlying linear systems, classical iterative methods have also been investigated to solve the interface problem. The numerical results presented here have shown that the convergence speed of the iterative methods are seriously affected by the presence of PML and the condition between the internal interfaces. A more accurate approximation of the Dirichlet-to-Neumann operator has thus been investigated. We have been able to numerically show that the use of a Mixed approach, which enables to vary the transmission condition in and out of the anisotropic areas, strongly enhances the convergence process and even enables it to no longer diverge.
So far, we have investigated a sequential implementation of the FETI-FDP2 method. Future work will naturally concern the algorithm's parallelization. In order to further improve the convergence rates, more advanced boundary conditions [37] should also be investigated. At the same time, thanks to our constant progress in microwave scattering experiments, we will analyse numerically and experimentally the scattering behaviour of objects whose overall dimensions are large with respect to the wavelength [38] [39]. This is of course an intermediate step before tackling complicated and large-scale inverse scattering problems [40] [41], where an efficient and accurate forward solver is a key feature.