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

Ivan Voznyuk, Hervé Tortel, Amelie Litman

► To cite this version:

Ivan Voznyuk, Hervé Tortel, Amelie Litman. 3-D Electromagnetic Scattering Computation in Free-Space With the FETI-FDP2 Method. IEEE Transactions on Antennas and Propagation, 2015, 63 (6), pp.2604 - 2613. 10.1109/TAP.2015.2417977. hal-01279552

HAL Id: hal-01279552

https://amu.hal.science/hal-01279552

Submitted on 5 Nov 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# 3D electromagnetic scattering computation in free-space with the FETI-FDP2 method

Ivan Voznyuk, Herve Tortel, and Amelie Litman, member, IEEE

**Abstract**—The electromagnetic dual-primal finite element tearing and interconnecting (FETI-DPEM) method is a non-overlapping 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 3D 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 are 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.

**Index Terms**—Finite Element Method (FEM), Domain decomposition method, FETI-DPEM method, non-conformal mesh, arbitrary partitioning, Krylov subspace iterative method, electromagnetic propagation and scattering, three-dimensional configuration, Perfectly Matched Layer (PML)

## 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 condition 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 authors are with Aix-Marseille Université/CNRS/Ecole Centrale Marseille - Institut Fresnel, UMR 7249 - Avenue Escadrille Normandie-Niemen - 13397 Marseille Cedex - France. e-mail: herve.tortel@fresnel.fr.
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 non-overlapping 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 developed 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 \subset \Omega$, 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 $\varepsilon_r^{\text{inc}}(\vec{r})$ and $\mu_r^{\text{inc}}(\vec{r})$ (resp. $\varepsilon_r^{\text{tot}}(\vec{r})$ and $\mu_r^{\text{tot}}(\vec{r})$) with

$$ \begin{aligned} \varepsilon_r^{\text{tot}}(\vec{r}) &\neq \varepsilon_r^{\text{inc}}(\vec{r}) & \forall \vec{r} \in D \\ \varepsilon_r^{\text{tot}}(\vec{r}) &= \varepsilon_r^{\text{inc}}(\vec{r}) & \forall \vec{r} \in \Omega \setminus D \end{aligned} \quad (1) $$

From the linearity of Maxwell's equations, the total field (associated to the permittivity distribution $\varepsilon_r^{\text{tot}}$) can be decomposed into the incident field $\vec{\mathcal{E}}^{\text{inc}}$ (associated to $\varepsilon_r^{\text{inc}}$) and the scattered field $\mathbf{E}^{\text{sc}}$ which satisfies a Helmholtz equation

$$ \vec{\nabla} \times \left( \frac{1}{\mu_r^{\text{tot}}} \vec{\nabla} \times \vec{\mathcal{E}}^{\text{sc}} \right) - k_0^2 \varepsilon_r^{\text{tot}} \vec{\mathcal{E}}^{\text{sc}} = \vec{\mathcal{J}}^{\text{sc}} \quad \text{in } \Omega \quad (2) $$

where $k_0$ is the vacuum wavenumber. The induced currents are defined as

$$ \vec{\mathcal{J}}^{\text{sc}} = -\vec{\nabla} \times \left( \left[ \frac{1}{\mu_r^{\text{tot}}} - \frac{1}{\mu_r^{\text{inc}}} \right] \vec{\nabla} \times \vec{\mathcal{E}}^{\text{inc}} \right) + k_0^2 \left[ \varepsilon_r^{\text{tot}} - \varepsilon_r^{\text{inc}} \right] \vec{\mathcal{E}}^{\text{inc}} \quad (3) $$

The scattered field also satisfies a radiation condition (RBC)

$$ \vec{n} \times \left( \frac{1}{\mu_r^{\text{tot}}} \vec{\nabla} \times \vec{\mathcal{E}}^{\text{sc}} \right) + j k_0 \vec{n} \times \vec{n} \times \vec{\mathcal{E}}^{\text{sc}} = 0 \quad \text{on } \partial\Omega \quad (4) $$

The computational domain can be bounded thanks to a Perfectly Matched Layer (PML) [22], [23] which can be interpreted as an anisotropic absorber [29] where the tensors $\mu_r^{\text{pml}}$ and $\varepsilon_r^{\text{pml}}$ are defined with respect to $\mu_r^{\text{tot}}$ and $\varepsilon_r^{\text{tot}}$ 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.
### III. FINITE ELEMENT TEARING AND INTERCONNECTING 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 $\Omega$ (with external boundary $\partial\Omega$) is divided into a set of $N_s$ non-overlapping subdomains $\Omega = \cup_{i=1}^{N_s} \Omega^i$, with $\Omega^i \cap \Omega^j = \emptyset, \forall i \neq j$. For a given subdomain $\Omega^i$, the indices of all its adjacent subdomains are called neighbor(i). The internal boundary of the subdomain $\Omega^i$ is denoted as $\Gamma^i$. The internal boundary between two subdomains $\Omega^i$ and $\Omega^j$, with $j \in$ neighbor(i), is denoted as $\Gamma^{ij}$.

Following the idea of the FETI approach [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 E^i = f^i - D^i \lambda^i \quad (5)$$

$K^i$ is the sum of the stiffness matrix, the mass matrix and the external and internal boundary condition matrices

$$K^i = \int_{\Omega^i} \left[ \frac{1}{\mu_r} (\nabla \times \vec{N}^i) \cdot (\nabla \times \vec{N}^i)^T - k_0^2 \epsilon_r^{\text{tot}} \vec{N}^i \cdot \vec{N}^i \right] dV \\ + \alpha^i \int_{\Gamma^i} (\vec{n} \times \vec{N}^i) \cdot (\vec{n} \times \vec{N}^i)^T dS \\ + jk_0 \int_{\partial\Omega^i \cap \partial\Omega} (\vec{n} \times \vec{N}^i) \cdot (\vec{n} \times \vec{N}^i)^T dS$$

where $\vec{N}^i$ denotes a column vector containing the curl-conforming vector basis function in the i-th subdomain. The role of $\alpha^i$ is explained in Section III-C. $E^i$ is the unknown components of the vector field $\vec{\mathcal{E}}^{\text{sc}}$ in $\Omega^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 $\Omega^i$ which are on $\Gamma^i$. As for the dual Lagrange multipliers $\lambda^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

$$E^i = \begin{bmatrix} E_r^i \\ E_c^i \end{bmatrix} \quad (6)$$

where the notations $E_r^i$ denote all the internal ($E_r^i$) and interface ($E_r^i$) Degrees of Freedom (DOFs) belonging to the i-th subdomain except for the corner DOFs which are denoted

by $E_c^i$. 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

$$\begin{bmatrix} K_{rr}^i & K_{rc}^i \\ K_{cr}^i & K_{cc}^i \end{bmatrix} \begin{bmatrix} E_r^i \\ E_c^i \end{bmatrix} = \begin{bmatrix} f_r^i \\ f_c^i \end{bmatrix} - \begin{bmatrix} D_r^i T \lambda_r^i \\ \lambda_c^i \end{bmatrix} \quad (7)$$

We now introduce several Boolean matrices to gather the DOFs of interest. The matrix $D_r^i$ selects from all the DOFs of $\Omega^i$ the ones which are located on its interface, i.e.,

$$D_r^i E_r^i = E_I^i \quad (8)$$

The projection matrix $T_r^{i \to j}$ extracts the DOFs associated to the interface $\Gamma^{ij}$ from the ones belonging to $\Gamma^i$,

$$T_r^{i \to j} \lambda_r^i = \lambda_r^{i \to j} \quad \text{and} \quad T_r^{i \to j} E_I^i = E_r^{i \to j} \quad (9)$$

A similar Boolean matrix $T_c^{i \to j}$ is here to select all the corner DOFs belonging to $\Gamma^i$ and associated to $\Gamma^{ij}$,

$$T_c^{i \to j} E_c^i = E_c^{i \to j} \quad (10)$$

Two supplemental Boolean matrices $Q_{\lambda_r}^i$ and $Q_{E_c}^i$ are introduced such that

$$Q_{\lambda_r}^i \lambda_r = \lambda_r^i \quad \text{and} \quad Q_{E_c}^i E_c = E_c^i \quad (11)$$

From the first line of (Eq. 7), the electric field in the i-th subdomain can be easily derived

$$E_r^i = K_{rr}^{i-1} (f_r^i - D_r^i T \lambda_r^i - K_{rc}^i E_c^i) \quad (12)$$

From the second line of (Eq. 7), the corner DOFS can be estimated

$$(K_{cc}^i - K_{cr}^i K_{rr}^{i-1} K_{rc}^i) Q_{E_c}^i E_c = f_c^i - \lambda_c^i - K_{cr}^i K_{rr}^{i-1} f_r^i + K_{cr}^i K_{rr}^{i-1} D_r^i T \lambda_r^i \quad (13)$$

By assembling and summing all the subdomains contributions, the global corner points related system equation is obtained as follows

$$F_{E_c E_c} E_c - F_{E_c \lambda_r} \lambda_r + F_{E_c \lambda_c} \lambda_c = d_{E_c} \quad (14)$$

with

$$F_{E_c E_c} = S_{E_c} (K_{cc}^i - K_{cr}^i K_{rr}^{i-1} K_{rc}^i) Q_{E_c}^i$$

$$F_{E_c \lambda_r} = S_{E_c} (K_{cr}^i K_{rr}^{i-1} D_r^i) Q_{\lambda_r}^i$$

$$d_{E_c} = S_{E_c} (f_c^i - K_{cr}^i K_{rr}^{i-1} f_r^i)$$

where the operator $S_{E_c}$ is defined by

$$S_{E_c} u = \sum_{i=1}^{N_s} Q_{E_c}^i u^i \quad (15)$$
### C. Boundary conditions at internal interfaces

A classical continuity boundary condition is first imposed on the tangential component of the electric field

$$ \llbracket \vec{n} \times \vec{\mathcal{E}}^{\text{sc}} \rrbracket_{i \leftrightarrow j} = 0 \text{ on } \Gamma^{ij} \quad (16) $$

where $\llbracket \phantom{\times} \rrbracket_{i \leftrightarrow j}$ corresponds to the jump boundary condition on $\Gamma^{ij}$. This equation is rewritten as

$$ E^{j \to i} = E^{i \to j} \quad \forall i, j \text{ s.t. } \Gamma^{ij} \neq \emptyset \quad (17) $$

where $E^{i \to j}$ (resp. $E^{j \to i}$) is the subset of $E^i$ associated to the edges of $\Omega^i$ (resp. $\Omega^j$) which belong to the interface $\Gamma^{ij}$. The introduction of dual Lagrangian multipliers enables to insert, at the internal interface $\Gamma^i$, either an unknown Neumann-type boundary condition (when $\alpha = 0$) or an unknown Robin-type boundary condition (when $\alpha \neq 0$)

$$ \vec{n} \times \left( \frac{1}{\mu_r^{\text{tot}}} \vec{\nabla} \times \vec{\mathcal{E}}^{\text{sc}} \right) + \alpha \vec{n} \times \vec{n} \times \vec{\mathcal{E}}^{\text{sc}} = \Lambda^i \text{ on } \Gamma^i \quad (18) $$

It leads to the following set of equations, taking advantages of (Eq. 17),

$$ \lambda^{i \to j} + \lambda^{j \to i} = -W^{i \leftrightarrow j} E^{j \to i} = -W^{i \leftrightarrow j} E^{i \to j} \quad (19) $$

where the matrix $W^{i \leftrightarrow j} = M^{i \to j} + M^{j \to i}$ directly translates the boundary conditions. Following the aforementioned "rc"- notations, this equation can be rewritten as follows:

$$ \begin{bmatrix} \lambda_r^{i \to j} \\ \lambda_c^{i \to j} \end{bmatrix} + \begin{bmatrix} \lambda_r^{j \to i} \\ \lambda_c^{j \to i} \end{bmatrix} = - \begin{bmatrix} W_{rr}^{i \leftrightarrow j} & W_{rc}^{i \leftrightarrow j} \\ W_{cr}^{i \leftrightarrow j} & W_{cc}^{i \leftrightarrow j} \end{bmatrix} \begin{bmatrix} E_r^{j \to i} \\ E_c^{j \to i} \end{bmatrix} \quad (20) $$

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

$$ \begin{bmatrix} \lambda_r^{i \to j} \\ \lambda_c^{i \to j} \end{bmatrix} + \begin{bmatrix} \lambda_r^{j \to i} \\ \lambda_c^{j \to i} \end{bmatrix} = - \begin{bmatrix} W_{rr}^{i \leftrightarrow j} & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} E_r^{j \to i} \\ E_c^{j \to i} \end{bmatrix} \quad (21) $$

Thanks to (Eq. 12), we can eliminate $E_r^{j \to i}$ 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:

$$ F_{\lambda_r \lambda_r} \lambda_r - F_{\lambda_r \mathbf{E}_c} \mathbf{E}_c = -d_{\lambda_r} \quad (22) $$

where

$$ F_{\lambda_r \lambda_r} = S_{\lambda_r} (\mathbb{T}_r^{j \to i} - W_{rr}^{i \leftrightarrow j} \mathbb{T}_r^{j \to i} F_{rr}^j) Q_{\lambda_r}^j + \mathbf{I} $$

$$ F_{\lambda_r \mathbf{E}_c} = S_{\lambda_r} (W_{rr}^{i \leftrightarrow j} \mathbb{T}_r^{j \to i} F_{rc}^j) Q_{\mathbf{E}_c}^j $$

$$ d_{\lambda_r} = S_{\lambda_r} (W_{rr}^{i \leftrightarrow j} \mathbb{T}_r^{j \to i} d_r^j) $$

$$ F_{rr}^i = \mathbb{D}_r^i K_{rr}^{i-1} \mathbb{D}_r^i T $$

$$ F_{rc}^i = \mathbb{D}_r^i K_{rr}^{i-1} K_{rc}^i $$

$$ d_r^i = \mathbb{D}_r^i K_{rr}^{i-1} f_r^i $$

and the operator $S_{\lambda_r}$ is defined by

$$ S_{\lambda_r} u = \sum_{i=1}^{N_s} Q_{\lambda_r}^{i T} \sum_{j \in \text{neighbor}(i)} \mathbb{T}_r^{i \to j} u^{j \to i} \quad (23) $$

Thanks to the definition of the operator $S_{\mathbf{E}_c}$, the second line of (Eq. 21) provides

$$ F_{\mathbf{E}_c \lambda_c} = S_{\mathbf{E}_c} \lambda_c^i = 0 \quad (24) $$

We can now assemble (Eq. 14), (Eq. 22) and (Eq. 24) to derive the FETI-DPEM2(a) Full Interface Problem

$$ \begin{bmatrix} F_{\lambda_r \lambda_r} & -F_{\lambda_r \mathbf{E}_c} \\ -F_{\mathbf{E}_c \lambda_r} & F_{\mathbf{E}_c \mathbf{E}_c} \end{bmatrix} \begin{bmatrix} \lambda_r \\ \mathbf{E}_c \end{bmatrix} = \begin{bmatrix} -d_{\lambda_r} \\ d_{\mathbf{E}_c} \end{bmatrix} \quad (25) $$

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_{rc}^{i \leftrightarrow j}$ 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

$$ \begin{bmatrix} \lambda_r^{i \to j} \\ \lambda_c^{i \to j} \end{bmatrix} + \begin{bmatrix} \lambda_r^{j \to i} \\ \lambda_c^{j \to i} \end{bmatrix} = - \begin{bmatrix} W_{rr}^{i \leftrightarrow j} & W_{rc}^{i \leftrightarrow j} \\ 0 & 0 \end{bmatrix} \begin{bmatrix} E_r^{j \to i} \\ E_c^{j \to i} \end{bmatrix} \quad (26) $$

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_{\lambda_r \mathbf{E}_c}$, which now corresponds to

$$ F_{\lambda_r \mathbf{E}_c} = S_{\lambda_r} (W_{rr}^{i \leftrightarrow j} \mathbb{T}_r^{j \to i} F_{rc}^j - W_{rc}^{i \leftrightarrow j} \mathbb{T}_c^{j \to i}) Q_{\mathbf{E}_c}^j $$

3) **FETI-FDP2 method:** In the framework of a two-dimensional configuration [21], we not only showed that the term $W_{rc}^{i \leftrightarrow j}$ 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 $\lambda_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

$$ \tilde{\lambda}_c^{i \leftrightarrow j} = \tilde{\lambda}_c^{j \to i} = \frac{1}{2} (\lambda_c^{i \to j} + \lambda_c^{j \to i}) \quad (27) $$

which results in a matrix $F_{\lambda_c \lambda_c} = 2I$.

In the current 3D configuration, we propose a more general framework. Instead of searching for $\lambda_c^{i \to j}$ and $\lambda_c^{j \to i}$, we search for global corner Lagrange multipliers $\lambda_c^{i \leftrightarrow j} = \lambda_c^{i \to j} + \lambda_c^{j \to i}$, where $i \leftrightarrow j$ is the index of the interface $\Gamma^{ij}$ among all the $N_I$ interfaces. Thus

$$ \begin{bmatrix} \lambda_r^{i \to j} + \lambda_c^{j \to i} \\ \lambda_c^{i \to j} \end{bmatrix} = - \begin{bmatrix} W_{rr}^{i \leftrightarrow j} & W_{rc}^{i \leftrightarrow j} \\ W_{cr}^{i \leftrightarrow j} & W_{cc}^{i \leftrightarrow j} \end{bmatrix} \begin{bmatrix} E_r^{j \to i} \\ E_c^{j \to i} \end{bmatrix} \quad (28) $$

A new projection Boolean matrix $Q_{\lambda_c}^{i \leftrightarrow j}$ is introduced to extract the double Lagrange multipliers $\lambda_c^{i \leftrightarrow j}$ associated to
$\Gamma^{ij}$ from the full list of $\lambda_c$, i.e., $\lambda_c^{i \leftrightarrow j} = \mathbb{Q}_{\lambda_c}^{i \leftrightarrow j} \lambda_c$. A new operator $S_{\lambda_c}$ is also defined such that [11]

$$S_{\lambda_c} u = \sum_{i \leftrightarrow j=1}^{N_I} \mathbb{Q}_{\lambda_c}^{i \leftrightarrow j} u^{T} u^{i \leftrightarrow j} = \sum_{i \leftrightarrow j=1}^{N_I} \mathbb{Q}_{\lambda_c}^{i \leftrightarrow j} u^{T} u^{j \leftrightarrow i} \quad (29)$$

By applying this new operator to the second line of (Eq. 28), we obtain

$$-F_{\lambda_c \lambda_r} \lambda_r - F_{\lambda_c E_c} E_c + F_{\lambda_c \lambda_c} \lambda_c = -d_{\lambda_c} \quad (30)$$

where

$$F_{\lambda_c \lambda_r} = S_{\lambda_c} (W_{cr}^{i \leftrightarrow j} T_r^{j \to i} F_{rr}^j) \mathbb{Q}_{\lambda_r}^j$$

$$F_{\lambda_c E_c} = S_{\lambda_c} (W_{cr}^{i \leftrightarrow j} T_r^{j \to i} F_{rc}^j - W_{cc}^{i \leftrightarrow j} T_c^{j \to i}) \mathbb{Q}_{E_c}^j$$

$$F_{\lambda_c \lambda_c} = S_{\lambda_c}$$

$$d_{\lambda_c} = S_{\lambda_c} (W_{cr}^{i \leftrightarrow j} T_r^{j \to i} d_r^j)$$

The matrix $F_{E_c \lambda_c}$ can not be introduced through the operator $S_{E_c}$. It is a Boolean projection matrix which puts in correspondence all the global corner Lagrange multipliers $\lambda_c^{i \leftrightarrow j}$ related to the global corner $E_c^k$ for every $k = 1...N_{E_c}$. The size of this matrix is equal to $N_{E_c} \times N_{\lambda_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.

Combining (Eq. 14), (Eq. 22) and (Eq. 30), we arrive at the FETI-FDP2 Full Interface Problem

$$\begin{bmatrix} F_{\lambda_r \lambda_r} & -F_{\lambda_r E_c} & 0 \\ -F_{E_c \lambda_r} & F_{E_c E_c} & F_{E_c \lambda_c} \\ -F_{\lambda_c \lambda_r} & -F_{\lambda_c E_c} & F_{\lambda_c \lambda_c} \end{bmatrix} \begin{bmatrix} \lambda_r \\ E_c \\ \lambda_c \end{bmatrix} = \begin{bmatrix} -d_{\lambda_r} \\ d_{E_c} \\ -d_{\lambda_c} \end{bmatrix} \quad (31)$$

### 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 $\lambda_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

$$(F_{\lambda_r \lambda_r} + F_{\lambda_r E_c} \hat{F}_{E_c E_c}^{-1} \hat{F}_{E_c \lambda_r}) \lambda_r = \hat{d}_{\lambda_r} \quad (32)$$

where

$$\hat{F}_{E_c \lambda_r} = -F_{E_c \lambda_r} + F_{E_c \lambda_c} F_{\lambda_c \lambda_c}^{-1} F_{\lambda_c \lambda_r}$$

$$\hat{F}_{E_c E_c} = F_{E_c E_c} + F_{E_c \lambda_c} F_{\lambda_c \lambda_c}^{-1} F_{\lambda_c E_c}$$

$$\hat{d}_{E_c} = d_{E_c} + F_{E_c \lambda_c} F_{\lambda_c \lambda_c}^{-1} d_{\lambda_c}$$

$$\hat{d}_{\lambda_r} = -d_{\lambda_r} + F_{\lambda_r E_c} \hat{F}_{E_c E_c}^{-1} \hat{d}_{E_c}$$

We can now easily compute $E_c$ by solving the following system of equations

$$\hat{F}_{E_c E_c} E_c = \hat{d}_{E_c} - \hat{F}_{E_c \lambda_r} \lambda_r \quad (33)$$

Finally, replacing $E_c$ and $\lambda_r$ in (Eq. 12) gives access to $E_r^i$ in each subdomain $\Omega^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 matrix $\hat{F}_{E_c E_c}$ 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.

### A. Configuration description

We are interested in the computation of the field scattered by a target consisting of four dielectric ellipsoids. Each scatterer has similar external dimensions ($D_x, D_y, D_z$) = (5, 8, 5) cm. Their positions are visible in Figure 1. Their relative permit-

<figure><img src="image_2.png" /><figcaption>Fig. 1. Schematic map of the scatterer inside the investigation domain.</figcaption></figure>

tivities $\epsilon_r^1, \epsilon_r^2, \epsilon_r^3$ and $\epsilon_r^4$ 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_0$.

The computational domain $\Omega$ is a parallelepiped whose size depends on the wavelength. For example, at 4GHz, the investigation domain $\Omega$ is a $[0.32 \times 0.34 \times 0.22]$ m$^3$ box ($\approx 61\lambda^3$, where $\lambda$ 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 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 $\hat{c}_{FEM}^{sc}$. The same frontal solver is used to solve
<figure><img src="image_3.png" /><figcaption>Fig. 2. Map of the computational domain, showing the way the partitioning has been performed for $N_s = 42$.</figcaption></figure>

the Reduced Interface Problems associated to each of the three aforementioned methods, with $\alpha = jk_0$. A $L_2$-norm error is introduced to compare the various scattered fields

$$L_2(\vec{\mathcal{E}}^{\text{sc}}) = \frac{\| \vec{\mathcal{E}}^{\text{sc}} - \vec{\mathcal{E}}_{\text{FEM}}^{\text{sc}} \|^2}{\| \vec{\mathcal{E}}_{\text{FEM}}^{\text{sc}} \|^2} \quad (34)$$

and the associated errors are provided in Table I.

<table><caption>TABLE I
$L_2$-NORM ERROR OBTAINED FOR VARIOUS PARTITIONING WHEN THE INTERFACE PROBLEMS ARE SOLVED WITH A FRONTAL METHOD ($f = 4$ GHz, $w_{\text{PML}} = \lambda$, $\alpha = jk_0$)</caption><thead><tr><th>$N_s$</th><th>FETI-DPEM2<sup>(a)</sup></th><th>FETI-DPEM2<sup>(b)</sup></th><th>FETI-FDP2</th></tr></thead><tbody><tr><td>5</td><td>5.11E-002</td><td>9.12E-013</td><td>2.92E-012</td></tr><tr><td>7</td><td>8.98E-002</td><td>1.59E-012</td><td>7.66E-013</td></tr><tr><td>10</td><td>1.73E-001</td><td>9.78E-013</td><td>3.46E-012</td></tr></tbody></table>

It is obvious that the relative error obtained with the FETI-DPEM2<sup>(a)</sup> method is significantly higher than the other ones. This comforts us in the fact that, as proposed in [21], the term $W_{rc}^{i \leftrightarrow j}$ of (Eq. 20) cannot be neglected. From now on, we will thus stop considering the FETI-DPEM2<sup>(a)</sup> method. The second obvious conclusion is that the FETI-DPEM2<sup>(b)</sup> 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 $\lambda$-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 $\Omega^{15}$.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Modulus of the Lagrange multipliers in $\Omega^{15}$ with ($f = 4$ GHz, $N_s = 42$, $w_{\text{PML}} = \lambda$, $\alpha = jk_0$). (Full line) Exact values extracted from $\mathcal{E}_{\text{FEM}}^{\text{sc}}$. (Dotted line) Values extracted from $\mathcal{E}_{\text{FDP2}}^{\text{sc}}$ after 10 GMRES iterations. The vertical grey lines denote the location of the PML area.</figcaption></figure>

The choice of this specific subdomain is due to the fact that, among all of them, the relative $L_2$-error in $\Omega^{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 $\alpha$ is simply set to $jk_0$ as in the Despres approach [35], the approximation
<figure><img src="image_5.png" /><figcaption>Fig. 4. Evolution (in dB) of the GMRES error along the iterations, when the external boundary conditions and the methods change ($f = 4$ GHz, $N_s = 42$, $\alpha = jk_0$).</figcaption></figure>

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 higher-order 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 $\alpha$ 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 $\alpha = jk_0(1+j\chi)$, $\chi$ being a real-valued positive coefficient. The optimal value for $\chi$ depends in particular from the mean curvature on the interface [37]. As the partitioning is performed in an irregular fashion, we simply set $\chi = 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 $\alpha = jk_0$ everywhere,

2) the EMDA approach, which states $\alpha = jk_0(1 + j\chi)$ everywhere,

3) a new Mixed approach, with $\alpha = jk_0(1 + j\chi)$ in the PML areas and $\alpha = 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. The number of eigenvalues which have negative real part

<figure><img src="image_6.png" /><figcaption>Fig. 5. Sketch of the various methods chosen for approximating the internal transmission coefficients. (Upper row) In presence of PML. (Lower row) In presence of RBC.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution (in dB) of the GMRES error along the iterations for the FETI-DPEM2(b) and FETI-FDP2 methods, in presence of RBC, when the internal boundary conditions change in terms of values and locations ($f = 4$ GHz, $N_s = 42$).</figcaption></figure>

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.

<figure><img src="image_8.png" /><figcaption>Fig. 7. Evolution (in dB) of the GMRES error along the iterations for the FETI-FDP2 method, in presence of PML, when the internal boundary conditions change in terms of values and locations ($f = 4$ GHz, $N_s = 42$, $w_{PML} = \lambda$).</figcaption></figure>

A similar case is considered, but this time with smaller PML ($w_{PML} = 0.4\lambda$). 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.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Eigenvalues of the Reduced Interface Problem related to the FETI-FDP2 method, in presence of PML, calculated for (a) the Despres approach, (b) the EMDA approach and (c) the Mixed approach. ($f = 4$ GHz, $N_s = 42$, $w_{PML} = \lambda$)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Evolution (in dB) of the GMRES error along the iterations for the FETI-DPEM2(b) and FETI-FDP2 methods, in presence of PML, when the internal boundary conditions change in terms of values and locations ($f = 4$ GHz, $w_{PML} = 0.4\lambda$).</figcaption></figure>

### 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

<figure><img src="image_11.png" /><figcaption>Fig. 10. Evolution (in dB) of the BiCGstab error along the iterations, in presence of PML, when the Mixed approach is applied for the two FETI methods ($f = 4$ GHz, $N_s = 42$, $w_{PML} = \lambda$).</figcaption></figure>

<table><caption>TABLE II
NUMBER OF ITERATIONS AND RELATED $L_2$-ERROR OBTAINED WITH THE BICGSTAB METHOD, IN PRESENCE OF RBC ($f = 4$ GHz, $N_s = 42$).</caption><thead><tr><th rowspan="2"></th><th colspan="2">FETI-DPEM2<sup>(b)</sup></th><th colspan="2">FETI-FDP2</th></tr><tr><th>N<sub>iter</sub></th><th>L<sub>2</sub>-error</th><th>N<sub>iter</sub></th><th>L<sub>2</sub>-error</th></tr></thead><tbody><tr><td>Despres/Mixed</td><td>111</td><td>8.60E-005</td><td>95</td><td>1.63E-004</td></tr><tr><td>EMDA</td><td>-</td><td>-</td><td>337</td><td>2.22E-004</td></tr></tbody></table>

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_{rr}^i$ is provided as well as the construction time and the factorization time required for the matrix $\hat{F}_{E_c E_c}$. In the FETI-DPEM2(b) method, only the computation of $F_{E_c E_c}$ is required. The addition of the new set of Lagrange multipliers $\lambda_c$ in the FETI-FDP2 method introduces 3 supplemental matrix-matrix multiplications, in order to obtain $\hat{F}_{\lambda_r E_c}$, $\hat{F}_{E_c E_c}$ and $\hat{d}_{E_c}$ 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_{\lambda_c \lambda_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_2$-error is thus not available. It is of interest to note that the $L_2$-error, 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.
<table><caption>TABLE III
RESULTS OBTAINED AT DIFFERENT FREQUENCIES IN PRESENCE OF RBC, WITH THE MIXED APPROACH. THE FETI-DPEM2(b) AND FETI-FDP2 METHODS ARE COMBINED WITH EITHER GMRES (GM) OR BICGSTAB (BCG) ITERATIVE SOLVERS.</caption><thead><tr><th rowspan="2">f (GHz)</th><th rowspan="2">Size of Ω (λ³)</th><th rowspan="2">Number of unknowns</th><th rowspan="2">N<sub>s</sub></th><th rowspan="2">N<sub>λr</sub></th><th rowspan="2">N<sub>Ec</sub></th><th colspan="2">FETI-DPEM2<sup>(b)</sup></th><th colspan="2">FETI-FDP2</th><th rowspan="2">K<sub>tr</sub><sup>i</sup> fact. time (sec)</th><th rowspan="2">F̂E<sub>tr</sub> E<sub>c</sub> constr/fact time (sec)</th><th rowspan="2">one GM iteration time (sec)</th></tr><tr><th>Iterations GM/BCG</th><th>L<sub>2</sub>-error</th><th>Iterations GM/BCG</th><th>L<sub>2</sub>-error</th></tr></thead><tbody><tr><td>4</td><td>61</td><td>403 638</td><td>42</td><td>89 464</td><td>2 563</td><td>9/111</td><td>4.02E-004</td><td>9/95</td><td>4.99E-004</td><td>21</td><td>58/3</td><td>10</td></tr><tr><td>8</td><td>143</td><td>943 997</td><td>99</td><td>219 188</td><td>6 768</td><td>10/104</td><td>9.34E-004</td><td>9/94</td><td>1.15E-003</td><td>43</td><td>145/67</td><td>23</td></tr><tr><td>12</td><td>232</td><td>1 893 660</td><td>199</td><td>461 562</td><td>15 209</td><td>12/224</td><td>n.a.</td><td>12/130</td><td>n.a.</td><td>103</td><td>406/756</td><td>56</td></tr><tr><td>16</td><td>327</td><td>2 580 141</td><td>145</td><td>501 628</td><td>13 011</td><td>14/473</td><td>n.a.</td><td>13/246</td><td>n.a.</td><td>172</td><td>632/464</td><td>89</td></tr><tr><td>20</td><td>547</td><td>3 188 231</td><td>177</td><td>616 942</td><td>15 943</td><td>17/427</td><td>n.a.</td><td>17/215</td><td>n.a.</td><td>205</td><td>778/860</td><td>109</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Map in the plane z = 0 of the amplitude of the scattered field computed at 20 GHz, in presence of RBC, computed with the FETI-FDP2 method and the Mixed approach. The location of the ellipsoidal scatterers are indicated with white lines.</figcaption></figure>

## 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 Robin-type 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.

## VI. ACKNOWLEDGMENTS

The authors would like to thank the editor and the referees for their judicious comments and suggestions as well as the deep analysis they have performed on the paper. They would also like to thank Dr. M. Xue for fruitful discussions.

## REFERENCES

[1] A.F. Peterson, S.L. Ray, and R. Mitra. *Computational Methods for Electromagnetics*. Oxford University Press, 1998.
[2] J. Jin. *The Finite Element Method in Electromagnetics*. John Wiley & Sons, 2002.

[3] J.L. Volakis, A. Chatterjee, and L.C. Kempel. *Finite Element Method for Electromagnetics with Application to Antennas, Microwave Circuits, and Scattering*. IEEE Press, 1998.

[4] F. Brezzi and L.D. Marini. A three-field domain decomposition method. *Contemporary Mathematics*, 157:27–34, 1993.

[5] K. Zhao, V. Rawat, and J.-F. Lee. A domain decomposition method for electromagnetic radiation and scattering analysis of multi-target problems. *IEEE Transactions on Antennas and Propagation*, 56(8):2211–2221, 2008.

[6] C. Farhat and J. Mandel. The two-level FETI method for static and dynamic plate problems - Part I: An optimal iterative solver for biharmonic systems. *Computer Methods in Applied Mechanics and Engineering*, 155(1-2):129–151, 1998.

[7] C. Farhat, P.S. Chen, J. Mandel, and F.X. Roux. The two-level FETI method - Part II: Extension to shell problems, parallel implementation and performance results. *Computer Methods in Applied Mechanics and Engineering*, 155(1-2):153–179, 1998.

[8] C. Farhat, A. Macedo, M. Lesoinne, F.X. Roux, F. Magoulès, and A. de la Bourdonnaie. Two-level domain decomposition methods with Lagrange multipliers for the fast iterative solution of acoustic scattering problems. *Computer Methods in Applied Mechanics and Engineering*, 184(2–4):213–239, 2000.

[9] C. Farhat, P. Avery, R. Tezaur, and J. Li. FETI-DPH: A dual-primal domain decomposition method for acoustic scattering. *Journal of Computational Acoustics*, 13(3):499–524, 2005.

[10] Y. Boubendir, X. Antoine, and C. Geuzaine. A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation. *Journal of Computational Physics*, 231(2):262–280, 2012.

[11] C. Farhat, A. Macedo, and M. Lesoinne. A two-level domain decomposition method for the iterative solution of high frequency exterior Helmholtz problems. *Numerische Mathematik*, 85(2):283–308, 2000.

[12] Y.-J. Li and J.-M. Jin. A vector dual-primal finite element tearing and interconnecting method for solving 3-D large-scale electromagnetic problems. *IEEE Transactions on Antennas and Propagation*, 54(10):3000–3009, 2006.

[13] Y.-J. Li and J.-M. Jin. A new dual-primal domain decomposition approach for finite element simulation of 3-D large-scale electromagnetic problems. *IEEE Transactions on Antennas and Propagation*, 55(10):2803–2810, 2007.

[14] Y.-J. Li and J.-M. Jin. Implementation of the second-order ABC in the FETI-DPEM method for 3D EM problems. *IEEE Transactions on Antennas and Propagation*, 56(8):2765–2769, 2008.

[15] V. Dolean, S. Lanteri, and R. Perrussel. A domain decomposition method for solving the three-dimensional time-harmonic Maxwell equations discretized by discontinuous Galerkin methods. *Journal of Computational Physics*, 227(3):2044–2072, 2008.

[16] V. Dolean, L. Gerardo-Giorda, and M.J. Gander. Optimized Schwarz methods for Maxwell's equations. *Siam Journal On Scientific Computing*, 31(3):2193–2213, 2009.

[17] V. Dolean, M. El Bouajaji, M.J. Gander, and S. Lanteri. Optimized Schwarz methods for Maxwell's equations with non-zero electric conductivity. *Domain Decomposition Methods In Science and Engineering XIX*, 78:269–276, 2011.

[18] C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, and D. Rixen. FETI-DP: a dual-primal unified FETI method - Part I: A faster alternative to the two-level FETI method. *International Journal For Numerical Methods In Engineering*, 50(7):1523–1544, 2001.

[19] M.-F. Xue and J.-M. Jin. Nonconformal FETI-DP methods for large-scale electromagnetic simulation. *IEEE Transactions on Antennas and Propagation*, 60(9):4291–4305, 2012.

[20] V. Dolean, S. Lanteri, and R. Perrussel. Optimized Schwarz algorithms for solving time-harmonic Maxwell's equations discretized by a discontinuous Galerkin method. *Ieee Transactions On Magnetics*, 44(6):954–957, 2008.

[21] I. Voznyuk, H. Tortel, and A. Litman. Scattered field computation with an extended FETI-DPEM2 method. *Progress In Electromagnetics Research*, 139:247–263, 2013.

[22] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves. *Journal of Computational Physics*, 114:185–200, 1994.

[23] J.-M. Jin and D.J. Riley. *Finite Element Analysis of Antennas and Arrays*. IEEE Press. John Wiley & Sons, 2009.

[24] Y. Boubendir, A. Bendali, and M.B. Fares. Coupling of a non-overlapping domain decomposition method for a nodal finite element

method with a boundary element method. *International Journal for Numerical Methods in Engineering*, 73:1624–1650, 2008.

[25] Z. Peng and J.F. Lee. Non-conformal domain decomposition method with second-order transmission conditions for time-harmonic electromagnetics. *Journal of Computational Physics*, 229(16):5615–5629, 2010.

[26] M.-F. Xue and J.-M Jin. A preconditioned dual-primal finite element tearing and interconnecting method for solving three-dimensional time-harmonic maxwell's equations. *Journal of Computational Physics*, 274:920–935, 2014.

[27] V. Rawat. *Finite element domain decomposition with second order transmission condition for time harmonic electromagnetic problem*. PhD thesis, The Ohio State University, 2009.

[28] Y. Saad and M. H. Schultz. GMRES: a generalised minimal residual algorithm for solving nonsymmetric linear systems. *SIAM Journal of Scientific and Statistical Computing*, 7:856–869, 1986.

[29] Z. S. Sacks, D.M. Kingsland, R. Lee, and J.-F. Lee. A perfectly matched anisotropic absorber for use as an absorbing boundary condition. *IEEE Transactions on Antennas and Propagation*, 43:1460–1463, 1995.

[30] Y.-J. Li and J.-M. Jin. Parallel implementation of the FETI-DPEM algorithm for general 3D EM simulations. *Journal of Computational Physics*, 228:3255–3267, 2009.

[31] P. R. Amestoy, I. S. Duff, and J.-Y. L'Excellent. Multifrontal parallel distributed symmetric and unsymmetric solvers. *Comput. Methods in Appl. Mech. Eng.*, 184:501–520, 2000. See internet address : http://graal.ens-lyon.fr/MUMPS.

[32] C. Geuzaine and J.F. Remacle. *GMSH: A three dimensional finite element mesh generator*. See internet address : http://www.geuz.org/gmsh/.

[33] Karipis. *METIS - Family of Multilevel Partitioning Algorithms*. See internet address : http://glaros.dtc.umn.edu/gkhome/views/metis.

[34] Y. Saad. *Iterative Methods for Sparse Linear Systems*. SIAM, 2003.

[35] B. Despres, P. Joly, and J. Roberts. A domain decomposition method for the harmonic Maxwell's equations. In *IMACS, International Symposium on Iterative Methods in Linear Algebra*, pages 475–484, 1992.

[36] M. Gander, F. Magoules, and F. Nataf. Optimized schwarz methods without overlap for the helmholtz equation. *SIAM Journal on Scientific Computing*, 24(1):38–60, 2002.

[37] Y. Boubendir, X. Antoine, and C. Geuzaine. A non-overlapping quasi-optimized Schwarz domain decomposition algorithm for the Helmholtz equation. In R. Bank, M. Holst, O. Widlund, and J. Xu, editors, *Domain Decomposition Methods in Science and Engineering XX*, number 91 in Lecture Notes in Computational Science and Engineering, pages 519–526. Springer, 2013.

[38] C. Eyraud, J.-M. Geffrin, and A. Litman. 3d-aggregate quantitative imaging : experimental results and polarization effects. *IEEE Transactions on Antennas and Propagation*, 59(4):1237–1244, 2011.

[39] C. Eyraud, R. Vaillon, A. Litman, J.-M. Geffrin, and O. Merchiers. Polarization effects in 3d vectorial-induced current reconstructions. *Journal of the Optical Society of America A*, 30(10):1967–1974, 2013.

[40] A. Litman and L. Crocco. Testing inversion algorithms against experimental data: 3d targets. *Inverse Problems*, 25(2):020201, 2009.

[41] C. Eyraud, A. Litman, A. Herique, and W. Kofman. Microwave imaging from experimental data within a bayesian framework with realistic random noise. *Inverse Problems*, 25(2):024005, 2009.

<figure><img src="image_13.png" /></figure>

Ivan Voznyuk was born in Russia in 1988. He received the M.S. degree in Applied Mathematics from the Novosibirsk State Technical University, Russian Federation and the Ph.D. in Physics from the University of Aix-Marseille, France in 2014. During his thesis, which was carried out at Institut Fresnel, he worked on the development of advanced electromagnetic numerical schemes for forward and inverse scattering applications. He is currently a postdoctoral researcher at CEA-LETI, Grenoble, France where his current research interests include the conception and the optimisation of various types of meta-materials.
PLACE
PHOTO
HERE

**Hervé Tortel** was born in France in 1970. He received the M.S. and Ph.D. degrees from the University of Provence, Marseille, France, in 1995. In 1996, he joined the Institut Fresnel as an Assistant Professor. His research interests include direct and inverse scattering problems in electromagnetism.

PLACE
PHOTO
HERE

**Amélie Litman** was born in France in 1972. She received the Ph.D. degree in Applied Mathematics from the University of Paris XI, Paris, France, in 1997. During 1997-1998, she was with the Eindhoven University of Technology, Eindhoven, The Netherlands, working in a postdoctoral position. From 1998 to 2002, she was with Schlumberger, France, where she worked on the development of inversion algorithms for oil prospection. In October 2002, she joined Institut Fresnel, Université Aix Marseille, CNRS, Ecole Centrale Marseille, France, as an Assistant Professor. Her research interests include forward and inverse scattering techniques.