Pierre Gosselet
email: gosselet@lmt.ens-cachan.fr

Christian Rey
email: rey@lmt.ens-cachan.fr

Non-overlapping domain decomposition methods in structural mechanics

The modern design of industrial structures leads to very complex simulations characterized by nonlinearities, high heterogeneities, tortuous geometries... Whatever the modelization may be, such an analysis leads to the solution to a family of large ill-conditioned linear systems. In this paper we study strategies to efficiently solve to linear system based on non-overlapping domain decomposition methods. We present a review of most employed approaches and their strong connections. We outline their mechanical interpretations as well as the practical issues when willing to implement and use them. Numerical properties are illustrated by various assessments from academic to industrial problems. An hybrid approach, mainly designed for multifield problems, is also introduced as it provides a general framework of such approaches.

INTRODUCTION

Hermann Schwarz is often referred to as the father of domain decomposition methods. In a 1869-paper he proposed an alternating method to solve a PDE equation set on a complex domain composed the overlapping union of a disk and a square (fig. 1), giving the mathematical basis of what is nowadays one of the most natural ways to benefit the modern hardware architecture of scientific computers.

In fact the growing importance of domain decomposition methods in scientific computation is deeply linked to the growth of parallel processing capabilities (in terms of number of processors, data exchange bandwidth between processors, parallel library efficiency, and of course performance of each processor). Because of the exponential increase of computational resource requirement for numerical simulation of more and more complex physical phenomena (non-linearities, couplings between physical mechanisms or between physical scales, random variables...) and more and more complex problems (optimization, inverse problems...), parallel processing appears to be an essential tool to handle the resulting numerical models.

Parallel processing is supposed to take care of two key-points of modern computations, the amount of operations and the required memory. Let us consider the simulation of a physical phenomenon, classically modeled by a PDE L(x) = f, x ∈ H(Ω). To take advantage of the parallel architecture of a calculator, a reflexion has to be carried out on how the original problem could be decomposed into collaborating subprocesses. The criteria for this decomposition will be: first, the ability to solve independent problems (on independent processors); second, how often processes have to be synchronized; and last what quantity of data has to be exchanged when synchronizing. When tracing back the idle time of resolution processes and analyzing hardware and software solutions, it is often observed that inter-processor communications are the most penalizing steps.

If we now consider the three great classes of mathematical decomposition of our reference problem which are operator splitting (for instance L = i L i [START_REF] Fortin | Méthodes de lagrangien augmenté -applications à la résolution numérique de problèmes aux limites[END_REF]), function-space decomposition (for instance H(Ω) = span(v i ), an example of which is modal decomposition) and domain decomposition (Ω = Ω i ), though the first two can lead to very elegant formulations, only domain decompositions ensure (once one subdomain or more have been associated to one processor) that independent computations are limited to small quantities and that the data to exchange is limited to the interface (or small overlap) between subdomains which is always one-to-one exchange of small amount of data.

So domain decomposition methods perfectly fit the criteria for building efficient algorithms running on parallel computers. Their use is very natural in engineering (and more precisely design) context : domain decomposition methods offer a framework where different design services can provide the virtual models of their own parts of a structure, each assessed independently, domain decomposition can then evaluate the behavior of the complete structure just setting specific behavior on the interface (perfect joint, unilateral contact, friction). Of course domain decomposition methods also work with one-piece structure (for instance fig. 2), then decomposition can be automated according to criteria which will be discussed later.

From an implementation point of view, programming domain decomposition methods is not an overwhelming task. Most often it can be added to existing solvers as a upper level of current code using the existing code as a black-box. The only requirement to implement domain decomposition is to be able to detect the interface between subdomains and use a protocol to share data on this common part. In this paper we will mostly focus on domain decomposition methods applied to finite element method [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]109], anyhow they can be applied to any discretization method (among others meshless methods [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF][START_REF] Breitkopt | Meshfree and particle based approaches in computational Mechanics[END_REF][START_REF] Liu | An Introduction to Meshfree Methods and Their Programming[END_REF] and discrete element methods [START_REF] D'addetta | A particle center based homogenization strategy for granular assemblies[END_REF][START_REF] Bolander | Irregular lattice model for quasistatic crack propagation[END_REF][START_REF] Delaplace | Fine description of fracture by using discrete particle model[END_REF]).

Though domain decomposition methods were more than one century old, they had not been extensively studied. Recent interest arose as they were understood to be well-suited to modern engineering and modern computational hardware. An important date in reinterest in domain decomposition methods is 1987 as first international congress dedicated to these methods occurred and DDM association was created (see http://www.ddm.org).

Yet the studies were first mathematical analysis oriented and emphasized on Schwarz overlapping family of algorithms. As interest in engineering problems grew, non-overlapping Schwarz and Schur methods, and coupling with discretization methods (mainly finite element) were more and more studied. Indeed, these methods are very natural to interpret mechanically, and moreover mechanical considerations often resulted in improvement to the methods. Basically the notion of interface between neighboring subdomains is a strong physical concept, to which is linked a set of conservation principles and phenomenological laws: for instance the conservation of fluxes (action-reaction principle) imposes the pointwise mechanical equilibrium of the interface and the equality of incoming mass (heat...) from one subdomain to the outgoing mass (heat...) of its neighbors; the "perfect interface" law consists in supposing that displacement field (pressure, temperature) is continuous at the interface, contact laws enable disjunction of subdomains but prohibit interpenetration.

In this context two methods arose in the beginning of the 90's : so-called Finite Element Tearing and Interconnecting (FETI) method [START_REF] Farhat | A method of finite tearing and interconnecting and its parallel solution algorithm[END_REF] and Balanced Domain Decomposition (BDD) [104]. From a mechanical point of view BDD consists in choosing the interface displacement field as main unknown while FETI consists in privileging the interface effort field. BDD is usually referred to as a primal approach while FETI is a dual approach. One of the interests of these methods, beyond their simple mechanical interpretation, is that they can easily be explained from a purely algebraic point of view (ie directly from the matrix form of the problem). In order to fit parallelization criteria, it clearly appeared that the interface problem should be solved using an iterative solver, each iteration requiring local (ie independent on each subdomain) resolution of finite element problem, which could be done with a direct solver. Then these methods combined direct and iterative solver trying to mix robustness of the first and cheapness of the second. Moreover the use of an iterative solver was made more efficient by the existence of relevant preconditioners (based on the resolution of a local dual problem for the primal approach and a primal local problem for the dual approach).

When first released, FETI could not handle floating substructures (ie substructures without enough Dirichlet conditions), thus limiting the choice of decomposition, while the primal approach could handle such substructures but with loss of scalability (convergence decayed as the number of floating substructures increased). A key point then was the introduction of rigid body motions as constraints and the use of generalized inverses. Because of its strong connections with multigrid methods [108], the rigid body motions constraint took the name of "coarse problem", it made the primal and dual methods able to handle most decompositions without loss of scalability [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Mandel | Balancing domain decomposition[END_REF]102]. From a mechanical point of view, the coarse problem enables non-neighboring subdomains to interact without requiring the transmission of data through intermediate subdomains, it then enables to spread global information on the whole structure scale.

Once equipped with their best preconditioners and coarse problems, mathematical results [START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF][START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF][START_REF] Brenner | Lower bounds in domain decomposition[END_REF] provide theoretical scalability of the methods. For instance for 3D elasticity problems, if h is the diameter of finite elements and H the diameter of subdomains, condition number κ of the interface problem reads (C is a real constant):

κ ≃ C 1 + log H h 2 (1)
which proves that the condition number only depends logarithmatically on the number of elements per subdomain. Many numerical assessment campaigns confirmed the good properties of the methods, their robustness compared to iterative solvers applied to the complete structure and their low cost (in terms of memory and CPU requirements) compared to direct solvers. Thus because they are well-suited to modern hardware (like PC clusters) they enable to achieve computations which could not be realized on classical computers because of too high memory requirement or too long computational time: these methods can handle problems with several millions of degrees of freedom.

Primal and dual methods were extended to heterogeneous problems by a cheap intervention on the preconditioners [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF] and on the initialization [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF], and to forth order elasticity (plates and shells) problems by the adjunction of so-called "second level problem" in order to regularize the displacement field around the corners of subdomains [START_REF] Farhat | The two-level FETI method for static and dynamic plate problemspart I: An optimal iterative solver for biharmonic systems[END_REF]103,[START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF]. As it became clear that the regularization of the displacement field was sufficient to suppress rigid body motions, specific algorithms which regularized a priori the subdomain problems were proposed: FETIDP [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF] and its primal counterpart BDDC [START_REF] Cros | A preconditioner for the schur complement domain decomposition method[END_REF], first in the plates and shells context, then in the second order elasticity context [START_REF] Klawonn | Some computational results for dual-primal feti methods for three dimensional elliptic problems[END_REF]. Now FETIDP and BDDC are considered as efficient as original FETI and BDD.

Methods were employed in many other contexts: transient dynamics [START_REF] Fragakis | The mosaic of high-performance domain decomposition methods for structural mechanics -part ii: Formulation enhancements, multiple righthand sides and implicit dynamics[END_REF], multifield problems (multiphysic problems such as porous media [START_REF] Gosselet | Une approche hybride de décomposition de domaine pour les problèmes multiphysiques : application à la poroélasticité[END_REF] and constrained problems such as incompressible flows [START_REF] Li | A dual-primal feti method for solving stokes/navier-stokes equations[END_REF][START_REF] Goldfeld | Balancing Neumann-Neumann for (in)compressible linear elasticity and (generalized) Stokes -parallel implementation[END_REF]), Helmotz equations [START_REF] De La Bourdonnaye | Advances in Computational Mechanics with High Performance Computing, chapter A method of finite element tearing and interconnecting for the Helmholtz problem[END_REF][START_REF] Farhat | FETI-H: a scalable domaine decomposition method for high frequency exterior Helmholtz problems[END_REF][START_REF] Farhat | A two-level domain decomposition method for the iterative solution of high frequency exterior Helmholtz problems[END_REF] and contact [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF][START_REF] Dostal | An overview of scalable feti-dp algorithms for variational inequalities[END_REF]. The use of domain decomposition methods in structural dynamic analysis is a rather old idea though it can now be confronted to well established methods in static analysis; the Craig-Bampton algorithm [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] is somehow the application of the primal strategy to such problems, the dual version of which was proposed in [START_REF] Rixen | A dual craig-bampton method for dynamic substructuring[END_REF], moreover ideas like the adjunction of coarse problems enabled to improve these methods.

Because of the strong connection between primal and dual approaches, some methods try to propose frameworks which generalize the two methods. The hybrid approach [START_REF] Gosselet | Méthodes de décomposition de domaine et méthodes d'accélération pour les problèmes multichamp en mécanique non-linéaire[END_REF][START_REF] Gosselet | A monolithic strategy based on an hybrid domain decomposition method for multiphysic problems. application to poroelasticity[END_REF] enables to select specific treatment (primal or dual) for each interface degree of freedom; if all degrees of freedom have the same treatment, the hybrid approach is exactly a classical approach. For certain multifield problems the hybrid approach enables to define physic-friendly solvers. The hybrid approach can also be obtained from specific optimality considerations [START_REF] Victorita Dolean | New constructions of domain decomposition methods for systems of pdes[END_REF]. Mixed approaches [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Nouy | Une stratégie de calcul multiéchelle avec homogénéisation en temps et en espace pour le calcul de structures fortement hétérogènes[END_REF][START_REF] Series | Une méthode de décomposition de domaine avec deux multiplicateurs de Lagrange[END_REF] consist in searching a linear combination of interface displacement and effort field, depending on the artificial stiffness introduced on the interface one can recover the classical approaches (null stiffness for the dual approach, infinite stiffness for the primal approach). Moreover, the mixed approaches enable to provide the interface mechanical behavior and provide a more natural framework to handle complex interfaces (contact, friction...) than classical approaches.

In this paper we aim at reviewing most of non-overlapping domain decomposition methods. To adopt a rather general point of view we introduce a set of notations strongly linked to mechanical considerations as it gives the interface the main role of the methods. We try to include all methods inside the same pattern so that we can easily highlight the connections and differences between them. We adopt a practical point of view as we describe the mechanical concepts, the algebraic formulations, the algorithms and the practical implementation of the methods. At each step we try to emphasize keypoints and not to avoid theoretical and practical difficulties.

This paper is organized as follows. In section 2 we introduce the mechanical framework of our study, the common notations and the notion of interface assembly operators and mechanical operators which will play a central role in the methods. Section 3 provides a rather extensive review of the nonoverlapping domain decomposition methods in the framework of discretized problems: basic primal and dual approaches (with their variations), three-field method for conforming grids, mixed and hybrid approaches. A keypoint of the previous methods is the adjunction of optional constraints to form a "coarse problem" which transmits global data through the whole structure, the strategies to introduce these optional constraints are studied in section 4, which leads to the definition of "recondensed" strategies FETIDP and BDDC. Section 5 deals with practical issues which are very often common to most of the methods. Assessments are given in section 6 to illustrate the methods and outline their main properties. Section 7 concludes the paper. As Krylov iterative solvers are often coupled to domain decomposition methods, main concepts and algorithms to use them are given in appendix A.

Formulation of an interface problem

To present as smoothly as possible non-overlapping domain decomposition methods we first consider a reference continuous mechanics problem, decompose the domain in two subdomains in order to introduce interface fields, then in order to describe correctly the interface we study a N -subdomain decomposition. Since our aim is not to prove theoretical performance results but make the reader feel some key-points of substructuring, we do not go too far in continuous formulation and quickly introduce discretized systems.

Reference problem

Ω g f u 0 ∂ g Ω ∂ u Ω Figure 3. Reference problem
Let us consider domain Ω in R n (n=1, 2 or 3) submitted to a classical linear elasticity problem (see figure 3) : displacement u 0 is imposed on part ∂ u Ω of the boundary of the domain, effort g is imposed on complementary part ∂ g Ω, volumic effort f is imposed on Ω, elasticity tensor is a [START_REF] Germain | Mécanique[END_REF][START_REF] Salençon | Mécanique des milieux continus[END_REF]. The system is governed by the following equations :

             div(σ) + f = 0 in Ω σ = a : ε(u) in Ω ε(u) = 1 2 grad(u) + grad(u) T in Ω σ.n = g on ∂ g Ω u = u 0 on ∂ u Ω (2)
In order to have the problem well posed, we suppose mes(∂ u Ω) > 0. We also suppose that tensor a defines a symmetric definite positive bilinear form on 2nd-order symmetric tensors. Under these assumptions, problem (2) has a unique solution [START_REF] Duvaut | Mécanique des milieux continus[END_REF].

Two-subdomain decomposition

Ω (1)
Ω (2) Υ Let us consider a partition of domain Ω into 2 substructures Ω (1) and Ω (2) . We define interface Υ between substructures (figure 4) :

Υ = ∂Ω (1)
∂Ω (2) (3) System ( 2) is posed on domain Ω, we write its restrictions to Ω (1) and Ω (2) :

s = 1 or 2,                div(σ (s) ) + f (s) = 0 in Ω (s) σ (s) = a (s) : ε(u (s) ) in Ω (s) ε(u (s) ) = 1 2 grad(u (s) ) + grad(u (s) ) T in Ω (s) σ (s) .n (s) = g (s) on ∂ g Ω ∂Ω (s) u (s) = u (s) 0 on ∂ u Ω ∂Ω (s) (4)
and the interface connection conditions, continuity of displacement u (1) = u (2) on Υ [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF] and equilibrium of efforts (action-reaction principle)

σ (1) n (1) + σ (2) n (2) = 0 on Υ (6)
Of course system [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF][START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF][START_REF] Bernardi | A new non conforming approach to domain decomposition: the Mortar Element Method[END_REF] is strictly equivalent to global problem (2).

N-subdomain decomposition

Let us consider a partition of domain Ω into N subdomains denoted Ω (s) . We can define the interface between two subdomains, the complete interface of one subdomain and the geometric interface at the complete structure scale:

   Υ (i,j) = Υ (j,i) = ∂Ω (i) ∂Ω (j) Υ (s) = j Υ (s,j) Υ = s Υ (s) (7)
When implementing the method, one (possibly virtual) processor is commonly assigned to each subdomain, hence because we can tell "local" computations (realized independently on each processor) from "global" computations (realized by exchanging data between processors), we often refer to values as being global or local. Then Υ (s) is the local interface Ω (1) Ω ( 2)

Ω (3) Υ (a) Primal (geometric) interface Ω (1)
Ω ( 2)

Ω (3) Υ (1,2)
Υ (1,3) Υ (2,3) (b) Dual interface (connectivity)

Figure 5. Definition of the interface for a N -subdomain decomposition and Υ the global interface. Because exchanges are most often one-to-one, Υ (i,j) is the (ij)-communication interface.

Using more than two subdomains (except when using "band"-decomposition) leads to the appearance of "multiple-points" also called "crosspoints" (which are nodes shared by more than two subdomains). These crosspoints lead to the existence of two descriptions of the interface (figure 5): so-called geometric interface Υ and so-called connectivity interface made out of the set of one-to-one interfaces (Υ (i,j) ) 1 i<j N . Each of the two most classical methods exclusively uses one of these descriptions so the geometric interface Υ is often referred to as the primal interface while the connectivity interface Υ is referred to as the dual interface.

How crosspoints are handled is a fundamental key in the differentiation of domain decomposition methods. In the remainder of the paper we will always refer to data attached to the dual interface using underlined notation.

Remark 2.1 Reader may have observed that the above presented connectivity description is redundant for crosspoints: let x be a crosspoint if x belongs to Υ (1,2) and Υ (2,3) , it of course belongs to Υ (1,3) . In the general case of a m-multiple crosspoint, there are m 2 connectivity relationships while only (m-1) would be sufficient and necessary. We will present strategies to remove these redundancies in the algebraic analysis of the methods.

Remark 2.2 Cross-points may also introduce, at the continuous level, punctual-interfaces in 2d or edge-interfaces in 3d, which are interface with zero measure. Most often from a physical point of view these are not considered as interfaces. Anyhow after discretization all relationships are written node-to-node and the problem no longer exists.

Discretization

We suppose that the reference problem has been discretized, leading to the resolution of n × n linear system:

Ku = f (8)
Because of its key role in structural mechanics we will often refer to finite element discretization though any other technique would suit. The key points are the link between matrix K and the domain geometry and the sparse filling of matrix K (related to the fact that only narrow nodes have non-zero interaction). We restrict to the case of element-oriented decompositions (each element belongs to one and only one substructure) which are conforming to the mesh which implies three conditions [START_REF] Rixen | Substructuring and dual methods in structural analysis[END_REF]:

• there is a one-to-one correspondence between degrees of freedom by the interface ;

• approximation spaces are the same by the interface ;

• models (beam, shell, 3d...) are the same by the interface.

Under these assumptions connection conditions simply write as node equalities. For nonconforming meshes, a classical solution is to use mortar-elements for which continuity and equilibrium are verified in a weak sense [START_REF] Achdou | Substructuring preconditioners for finite element methods on nonmatching grids[END_REF][START_REF] Achdou | Méthode itérative de sous-structuration pour les éléments avec joints[END_REF]101].

Boolean operators

In order to write communication relation between subdomains we have to introduce several operators.

The first one is the "local trace" operator t (s) which is the discrete projection from Ω (s) to Υ (s) . It enables to cast data from a complete subdomain to its interface, and once transposed to extend data set on the interface to the whole subdomain (setting internal degrees of freedom to zero). In the remainder of the paper we will use subscript b for interface data and subscript i for internal data.

Then data lying on one subdomain interface has to be exchanged with its neighboring subdomains. It can be either realized on the primal interface or the dual interface, leading to two (global) "assembly" operators: the primal one A (s) , and the dual one A (s) . The primal assembly operator is a strictly boolean operator while the dual assembly operator is a signed boolean operator (see figure 6): if a degree of freedom is set to 1 on one side of the interface, its corresponding degree of freedom on the other side of the interface is set to -1. Non-boolean assembly operators can be used in order to average connection conditions when using non-conforming domain decomposition methods [START_REF] Bernardi | A new non conforming approach to domain decomposition: the Mortar Element Method[END_REF].

Remark 2.3 Our (t (s) , A (s) , A (s) ) set of operators is not the most commonly used in papers related to domain decomposition. The interest of this choice is to be sufficient to explain most of the available strategies with only three operators. Other notations use "composed operators" like (B (s) = A (s) t (s) or L (s) T = A (s) t (s) ) which are not sufficient to describe all methods and which, in a way, omit the fundamental role played by the interface.

Boolean operators have important classical properties. Please note the first one which expresses the orthogonality of the two assembly operators. s A (s) A (s) T = 0 (9a)

A (s) T A (s) = I Υ (s) (9b)
A (s) T A (s) = diag(multiplicity -1) Υ (s) (9c)
A (s) A (s) T = I on Υ (s) 0 elsewhere (9d)

Remark 2.4 An interesting choice of description would have been to use redundant local interface (defining some kind of t (s) ). This choice would stick to most classical implementations where the local interface of one subdomain is defined neighborwise. Dual assembly operator would write easily as a simple signing operator, but handling multiple points would be slightly more difficult for the primal assembly operator (see section 5.6).

Remark 2.5 Redundancies can easily be removed from the dual description of the interface. One has just to modify the connectivity table, so that one multiple point is connected only once to each subdomain. This can be carried out introducing two different assembly operators the "non-redundant" one and the "orthonormal" one (see figure 7) [START_REF] Fragakis | A unified framework for formulating domain decomposition methods in structural mechanics[END_REF]. Only relationship 9c is modified (then A (s) T A (s) = I Υ (s) ). The interest of the use of these assembly operators will be discussed in section 3.3.

(d) Dual interface

t (1) = 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 t (2) = 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 t (3) = 1 0 0 0 0 1 0 0 0 0 1 0 A (1) =    0 0 0 0 1 0 1 0 0 0 0 1    A (2) =    1 0 0 0 0 1 0 0 0 0 1 0    A (3) =    0 0 1 0 0 0 1 0 0 0 1 0    A (1) =       0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 1       A (2) =       1 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0 0       A (3) =       0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 0 0 -1 0       Figure 6.
Local numberings, interface numberings, trace and assembly operators 1

(1) b 2 (1) b 3 (1) b 1 (2) b 2 (2) b 3 (2) b 1 (3) b 2 (3) b 3 (3) b (a) Local interface 1 Υ 2 Υ 3 Υ 4 Υ 5 Υ 6 Υ (b) Redundant connectivity 1 Υ 2 Υ 3 Υ 4 Υ 5 Υ (c) Non-redundant connectivity 1 Υ 2 Υ 3 Υ 4 Υ 5 Υ (d) Orthonormal connectivity A (1) =       0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 1       A (2) =       1 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0 0 0 0       A (3) =       0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 0 0 -1 0       A (1) N =      0 0 0 0 1 0 1 0 0 0 0 0 0 0 1      A (2) N =      1 0 0 0 0 -1 0 0 0 0 1 0 0 -1 0      A (3) N =      0 0 -1 0 0 0 -1 0 0 0 -1 0 0 0 0      A (1) O =        0 0 0 0 1 √ 2 0 1 √ 2 0 0 0 0 0 0 0 2 √ 6        A (2) O =        1 √ 2 0 0 0 0 -1 √ 2 0 0 0 0 1 √ 2 0 0 -1 √ 6 0        A (3) O =        0 0 -1 √ 2 0 0 0 -1 √ 2 0 0 0 -1 √ 2 0 0 -1 √ 6 0        Figure 7
. Suppressing redundancies of dual interface

Basic equations

In order to rewrite equation ( 8) in a domain-decomposed context, we have to introduce the reaction unknown which is the discretization of σ (1) n (1) = -σ (2) n (2) in equation ( 6).

λ (s) is the reaction imposed by neighboring subdomains on subdomain (s). Commonly λ (s) is defined on the whole subdomain (s) while it is non-zero only on its interface, so

λ (s) = t (s) T λ (s) b . ∀s, K (s) u (s) = f (s) + λ (s) (10a) s A (s) t (s) u (s) = 0 (10b) s A (s) t (s) λ (s) = 0 (10c)
Equation (10a) corresponds to the (local) equilibrium of each subdomain submitted to external conditions f (s) and reactions from neighbors λ (s) . Equation (10b) corresponds to the (global) continuity of the displacement field through the interface. Equation (10c) corresponds to the (global) equilibrium of the interface (action-reaction principle).

This three-equation system [START_REF] Brenner | Lower bounds in domain decomposition[END_REF] is the starting point from a rich zoology of methods, most of which possess strong connections we will try to emphasis on. Before going further in the exploration of these methods, we propose to introduce local condensed operators that represent a subdomain on its interface, then a set of synthetic notations.

Local condensed operators

Philosophically, local condensed operators are operators that represent how neighboring subdomains "see" one subdomain: a subdomain can be viewed as a black-box, the only information necessary for neighbors is how it behaves on its interface. Associated to this idea is the classical assumption that local operations are "exactly" performed. From an implementation point of view, when solving problems involving local matrices, a direct solver is employed. As we will see, the use of exact local solvers will be coupled with the use of iterative global solvers leading to a powerful combination of speed and precision of computations.

In this section we will always refer to the local equilibrium of subdomain (s) under interface loading:

K (s) u (s) = λ (s) = t (s) T λ (s) b (11)
Primal Schur complement S (s)

p : If we renumber the local degrees of freedom of subdomain (s) in order to separate internal and boundary degrees of freedom, system [START_REF] Brezzi | A three-field domain decomposition method[END_REF] writes

K (s) ii K (s) ib K (s) bi K (s) bb u (s) i u (s) b = 0 λ (s) b (12)
From the first line we draw

u (s) i = -K (s) ii -1 K (s) ib u (s) b (13)
then the Gauss elimination of u (s)

i leads to

K (s) bb -K (s) bi K (s) ii -1 K (s) ib u (s) b = S (s) p u (s) b = λ (s) b (14)
which is the condensed form of the local equilibrium of subdomains expressed in terms of interface fields. Operator S (s)

p is called local primal Schur complement. Its computation is realized by the inversion of matrix K (s)

ii which corresponds to Dirichlet conditions imposed on the interface of subdomain (s), so the primal Schur complement is always well defined, and commonly called the "local Dirichlet operator". Note that the symmetry, positivity, and definition properties are inherited by matrix S (s) p from matrix K (s) .

An important result is that the kernel of matrices K (s) and S is the identity matrix on the interface):

K (s) R (s) = 0 =⇒ S (s) p t (s) R (s) = S (s) p R (s) b = 0 (15)
S (s) p R (s) b = 0 =⇒ K (s) -K (s) ii -1 K (s) ib I (s) b R (s) b = K (s) R (s) = 0 (16)
Primal Schur complement can also be interpreted as the discretization of the Stecklov-Poincaré operator. From a mechanical point of view, it is the linear operator that provides the reaction associated to given interface displacement field.

If we consider that the subdomain is also loaded on internal degrees of freedom, then the condensation of the equilibrium on the interface reads:

K (s) u (s) = f (s) ⇒ S (s) p u (s) b = b (s) p (17) b (s) p = f (s) b -K (s) bi K (s) ii -1 f (s) i (18) b (s)
p is the condensed effort imposed on the substructure.

Dual Schur complement S (s)

d : the dual Schur complement is a linear operator that computes interface displacement field from given interface effort. From equation [START_REF] Brezzi | A three-field domain decomposition method[END_REF] and [START_REF] Chapman | Deflated and augmented Krylov subspace techniques[END_REF] we have:

t (s) K (s) + t (s) T λ (s) b = S (s) p + λ (s) b = S (s) d λ (s) b = u (s) b (19)
where K (s) + is the generalized inverse of matrix K (s) , and where it is assumed that no rigid body motion is excited. If we denote by R (s) the kernel of matrix K (s) this last condition reads:

R (s) T λ (s) = 0 or equivalently R (s) b T λ (s) b = 0 ( 20
)
Remark 2.6 A generalized inverse (or pseudo-inverse) of matrix M is a matrix, denoted M + , which verifies the following property: ∀y ∈ range(M ), M M + y = y. Note that this definition leads to non-unique generalized inverse, however all results presented are independent of the choice of generalized inverse.

Of course in order to take into account, inside [START_REF] Cros | A preconditioner for the schur complement domain decomposition method[END_REF], the possibility of the substructure to have zero energy modes, an arbitrary rigid displacement can be added leading to the next expression where vector α (s) denotes the magnitude of rigid body motions:

u (s) b = S (s) d λ (s) b + R (s) b α (s) (21)
Hybrid Schur complement: S (s) pd this operator corresponds to an interface where degrees of freedom are partitioned into two subsets. Suppose that the first subset is submitted to given Dirichlet conditions and the second to Neumann condition, S (s) pd is the linear operator that associates resulting reaction on the first subset and resulting displacement on the second subset to those given conditions. We denote by subscript p data defined on the first subset and by subscript d data defined on the second subset (schematically b = p ∪ d and p ∩ d = ∅).

S (s) pd u (s) p λ (s) d = λ (s) p u (s) d (22)
The computation of operator S (s) pd though no more complex than the computation of S (s) p or S (s) d , requires more notations. A synthetic option is to denote by subscript p the sets of internal (subscript i) data and second-interface-subset (subscript d) data, schematically p = i ∪ d. We introduce a modified trace operator:

t (s) d v p = t (s) d v i v d = v d (23)
Then internal equilibrium (11) reads:

K (s) pp K (s) pp K (s) pp K (s) pp u (s) p u (s) p = t (s) d T λ (s) d λ (s) p (24)
Then hybrid Schur complement is:

S (s) pd = K (s) pp -K (s) pp K (s) pp + K (s) pp K (s) pp K (s) pp + t (s) d T -t (s) d K (s) pp + K pp t (s) d K (s) pp + t (s) d T (25)
As can be noticed the diagonal blocks of S (s)

pd look like fully primal and fully dual Schur complements, while extradiagonal blocks are antisymmetric (assuming K (s) is symmetric). Of course if all interface degrees of freedom belong to the same subset, the hybrid Schur complement equals "classical" fully primal or fully dual Schur complement. Moreover it stands out clearly that:

S (s) pd + = S (s) dp =   t (s) p K (s) d d + t (s) p T -t (s) p K (s) d d + K dd K (s) d d K (s) d d + t (s) p T K (s) dd -K (s) d d K (s) d d + K (s) dd   (26)

S (s)

dp is the operator which associates displacement on the first subset and reaction on the second subset to given effort on the first subset and given displacement on the second subset.

As both matrices K (s) pp and K (s) d d may not be invertible, only their pseudo-inverse has been introduced. The invertibility is strongly dependent on the choice of interface subsets.

Block notations

While condensed operators simplify the writing of local operations, the block notations make it easer to understand the global operations of domain decomposition. We propose to denote by superscript ⋄ the row-block repetition of local vectors and the diagonal-block repetition of matrices, block assembly operators are written in one row (column-block) and denoted by special font, for instance:

u ⋄ =   
u (1) . . .

u (N )    f ⋄ =   
f (1) . . .

f (N )    λ ⋄ =   
λ (1) . . .

λ (N )    K ⋄ =       K (1) 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 K (N )       t ⋄ =       t (1) 0 . . . 0 0 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 t (N )       A = A (1) . . . A (N ) A = A (1) . . . A (N )
Remark 2.7 The specific notation for assembly operators aims at emphasizing at their specific role in terms of parallelism for the methods. Moreover, the only operation that requires exchange of data between subdomains is the use of non-transposed assembly operators.

Fundamental system (10) then reads:

K ⋄ u ⋄ = f ⋄ + λ ⋄ (27a) At ⋄ λ ⋄ = 0 (27b) At ⋄ u ⋄ = 0 (27c)
or in condensed form:

S ⋄ p u ⋄ b = b ⋄ p + λ ⋄ b (28a) Aλ ⋄ b = 0 (28b) Au ⋄ b = 0 (28c)
The orthogonal property of assembly operators (9a) simply reads:

AA T = 0 (29)
Relation (9d) reads:

AA T = diag(multiplicity) (30)
Remark 2.8 For improved readability, we will denote by bold font objects defined in a unique way on the interface (ie "assembled" quantities). Schematically, assembly operators enable to go from block notations to bold notations and transposed assembly operators realize the reciprocal operations.

Brief review of classical strategies

We can define general strategies to solve system ( 27) or ( 28):

Primal approaches [104, 105, 74, 102, 75, 106, 103] a unique interface displacement unknown u b satisfying equation (28c) is introduced, then an iterative process enables to satisfy (28b) while always verifying (28a).

Dual approaches [START_REF] Farhat | A method of finite tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | A saddle-point principle domain decomposition method for the solution of solid mechanics problems[END_REF][START_REF] Farhat | The dual Schur complement method with well-posed local Neumann problems[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Farhat | Optimal convergence properties of the FETI domain decomposition method[END_REF][START_REF] Mandel | Convergence of a substructuring method with Lagrange multipliers[END_REF][START_REF] Bhardwaj | Application of the FETI method to ASCI problems: Scalability results on a thousand-processor and discussion of highly heterogeneous problems[END_REF] a unique interface effort unknown λ b satisfying equation (28b) is introduced, then an iterative process enables to satisfy (28c) while always verifying (28a).

Three fields approaches [START_REF] Brezzi | A three-field domain decomposition method[END_REF][START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: algorithm description[END_REF][START_REF] Rixen | Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver[END_REF] a unique interface displacement u b is introduced, then relation (28c) is dualized so that interface efforts λ ⋄ b are introduced as Lagrange multipliers which yet have to verify relation (28b). Then the iterative process looks simultaneously for (λ ⋄ b , u b , u ⋄ b ) verifying exactly equation (28a). As this method is mostly designed for non-matching discretizations it will not be exposed in the remaining of this paper, anyhow a variant of the dual method which is equivalent to the three-field method with conforming grids will be described.

Mixed approaches [START_REF] Glowinski | Augmented lagrangian interpretation of the nonoverlapping Schwartz alternating method[END_REF][START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF]100,[START_REF] Series | Une méthode de décomposition de domaine avec deux multiplicateurs de Lagrange[END_REF] new interface unknown is introduced which is a linear combination of interface displacement and effort,

µ ⋄ b = λ ⋄ b + T ⋄ b u ⋄ b ,
then the interface system is rewritten in terms of unknown µ ⋄ b , this new system is solved iteratively and then λ ⋄ b and u ⋄ b are postprocessed. Of course matrix T ⋄ b is an important parameter of these methods.

Hybrid approaches [START_REF] Klawonn | Dual and dual-primal FETI methods for elliptic problems with discontinuous coefficients[END_REF][START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF][START_REF] Mandel | On the convergence of a dual-primal substructuring method[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF][START_REF] Gosselet | A monolithic strategy based on an hybrid domain decomposition method for multiphysic problems. application to poroelasticity[END_REF] interface is split into parts where primal, dual or mixed approaches are applied, specific recondensation methods may then be applied.

Many different methods can be deduced from these large strategies, the most common will be presented and discussed in section 3. Anyhow since iterative solvers are used to solve interface problems, we recommend the reader to refer to appendix A where most used solvers are presented, including important details about constrained resolutions.

Classical solution strategies to the interface problem

The aim of this section is to give extended review of classical domain decomposition methods, the principle of which has just been exposed. The association with Krylov iterative solvers is an important point of these methods, appendix A provides a summary of important results and algorithms that are used in this section.

Primal domain decomposition method

The principle of primal domain decomposition method is to write the interface problem in terms of one unique unknown interface displacement field u b . The trace of local displacement fields then writes u ⋄ b = A T u b . Because of the orthogonality between assembly operators [START_REF] Erhel | Restarted GMRes preconditioned by deflation[END_REF], equation (28c) is automatically verified. Using equation (28b) to eliminate unknown reaction λ ⋄ b inside (28a), we get the primal formulation of the interface problem:

S p u b = AS ⋄ p A T u b = Ab ⋄ p = b p (31)
Operator S p is the global primal Schur complement of the decomposed structure, it results as the sum of local contributions (with non-block notations S p = s A (s) S (s) p A (s) T ). Using a direct solver to solve system (31) implies the exact computation of local contribution, the sum of these contributions (in a parallel computing context, this step correspond to large data exchange between processors) and the inversion of the global primal Schur complement which size is the global geometric interface (the size of which is far from being neglectable) and which sparsity is very poor (each interface degree of freedom is connected to degrees of freedom belonging to the same subdomains). Using an iterative solver is much less expensive since the only required operations are matrix-vector products which can be realized in parallel because of the assembled structure of global primal Schur complement; moreover excellent and rather cheap preconditioner exists. Note that if global matrix K is symmetric positive definite then so is operator S p and then popular conjugate gradient algorithm can be used to solve the primal interface problem, in other cases solvers like GMRes or orthodir have to be employed.

Preconditioner to the primal interface problem

An efficient preconditioner S-1 p is an interface operator giving a good approximation of the inverse of S p . Various strategies are possible. For instance, a direct preconditioning method is based on the construction of an approximate Schur complement from a simplified structure defined by degrees of freedom "near" the interface. Anyhow such a method does not respect the repartition of the data through processors. A good parallel preconditioner has to minimize data exchange.

Since operator S p is the sum of local contributions, the most classical strategy is then to define S-1 p as a scaled sum of the inverse of local contributions:

S-1 p = ÃS ⋄ p + ÃT = ÃS ⋄ d ÃT (32)
Since S ⋄ p + = S ⋄ d requires the computation of local problems with given effort on the interface, this preconditioner is called the Neumann preconditioner. Scaled assembly operator à can be defined the following way [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF]:

à = AM ⋄ A T -1 AM ⋄ ( 33
)
where M ⋄ is a parameter which enables to take into account the heterogeneity of the subdomains connected by the interface. It should make matrix AM ⋄ A T easily invertible and give a representation of the stiffness of the interface, most commonly:

• M ⋄ = I ⋄ for homogeneous structures, • M ⋄ = diag(K ⋄ bb )
for compressible heterogeneous structures,

• M ⋄ = µ ⋄ for incompressible heterogeneous structures (µ ⋄ is the diagonal matrix the coefficients of which are the shearing modulus of interface degrees of freedom).

The (s) notation makes it easier to understand implementation of scaled assembly operators:

S-1 p = s M (s) A (s) S (s) d A (s) T M (s) (34)
• M (s) = diag( 1 multiplicity ) for homogeneous structures, • M (s) = diag( diag(K (s) bb ) i j diag(K (j) bb ) i
) for compressible heterogeneous structures (assuming i represents the same degree of freedom shared by the j subdomains)

• M (s) = diag( µ (s) i j µ (j) i
) for incompressible heterogeneous structures (assuming i represents the same degree of freedom shared by the j subdomains)

The following partition of unity result clearly holds:

ÃA T = I Υ (35) s M (s) = I Υ (36)

Coarse problem

The use of dual Schur complement is associated to an optimality condition, as said earlier vector being multiplied by the pseudo inverse should lie inside the image of S ⋄ p . Since preconditioning is applied to residual r, the optimality condition reads:

R ⋄ b T ÃT r = 0 ( 37
)
and introducing classical notation G = ÃR ⋄ b , G T r = 0. Such a condition can then be interpreted as an augmented-Krylov algorithm (see section A.7). Once equipped with that augmentation problem, the primal Schur complement method is referred to as the "balanced domain decomposition" (BDD [START_REF] Mandel | Balancing domain decomposition[END_REF]102]). Algorithm 3.1 summarizes the classical BDD approach, and figure 8 provides a schematic representation of the first iteration of the preconditioned primal approach. Algorithm 3.1 Primal Schur complement with conjugate gradient

1: Set P = I -G(G T S p G) -1 G T S p 2: Compute u 0 = G(G T S p G) -1 G T b p 3: Compute r 0 = b p -S p u 0 = P T b p 4: z 0 = S-1
p r 0 set w 0 = z 0 5: for j = 0, . . . , m do 6: p j = S p P w j (notice S p P = P T S p = P T S p P ) 7:

α j = (z j , r j )/(p j , w j ) 8: u j+1 = u j + α j w j 9:
r j+1 = r jα j p j 10:

z j+1 = S-1 p r j+1 11:
For 0 i j,

β i j = -(z j+1 , p i )/(w i , p i ) 12:
w j+1 = z j+1 + j i=1 β i j w i 13: end for

Error estimate

The reference error estimate is the one linked to the convergence over the complete structure:

Ku-f f
. Assuming local inversions are exact, we reach the following result:

Ku -f f = S p u b -b p f (38)
During the iterative process S p u bb p is the norm of residual r as computed line 9 of algorithm 3.1, so the global convergence can be controlled by the convergence of the interface iterative process.

P-FETI method

The P-FETI method is a variation of BDD proposed by [START_REF] Fragakis | The mosaic of high-performance domain decomposition methods for structural mechanics -part i: Formulation, interrelation and numerical efficiency of primal and dual methods[END_REF][START_REF] Fragakis | The mosaic of high-performance domain decomposition methods for structural mechanics -part ii: Formulation enhancements, multiple righthand sides and implicit dynamics[END_REF] inspired by the dual approach (the reader should refer to the dual method before going further inside P-FETI). Its principle is to provide another assembly operator which incorporate rigid body elimination by a dual-like projector.

S- The choice of matrix Q is guided by the same considerations as in the dual method. It is worth noting that when Q is chosen equal to the Dirichlet preconditioner of the dual method (Q = ÃT S ⋄ p Ã) then the P-FETI method is equivalent to the classical balanced domain decomposition.

1 p = HS ⋄ d H T (39)
H T = ÃT -A T QG G T QG -1 G T (40)

Dual domain decomposition method

The principle of dual domain decomposition method is to write the interface problem in terms of one unique unknown interface effort field λ b . The trace of local effort fields then writes λ ⋄ b = A T λ b . Because of the orthogonality between assembly operators (29), equation (28b) is automatically verified. In order to eliminate unknown interface displacement field using (28c), we first obtain it from equation (28a) (or equivalently (27a)): as seen in ( 21) the inversion of local systems may require the use of generalized inverse and the introduction of rigid body motions the magnitude of which is denoted by vector α (s) , the use of generalized inverse is then submitted to compatibility condition.

u ⋄ b = S ⋄ d (b ⋄ p + A T λ b ) + R ⋄ b α ⋄ (41) R ⋄ b T (b ⋄ p + A T λ b ) = 0 (42)
The first line is then premultiplied by A (same expressions could be obtain from non condensed notations).

S d = AS ⋄ d A T b ⋄ d = S ⋄ d b ⋄ p = t ⋄ K ⋄+ f ⋄ G = AR ⋄ b e ⋄ = R ⋄ b T b ⋄ p = R ⋄T f ⋄ we get the dual formulation of the interface problem: S d G G T 0 λ b α ⋄ = -b d -e ⋄ (43)
This is the basic dual Schur complement method, also called Finite Element Tearing and Interconnecting method (FETI [START_REF] Farhat | A method of finite tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF]). For similar reasons to the primal Schur complement method, this system is most often solved using an iterative solver, then we will soon discuss the preconditioning issue and how the G T λ b + e ⋄ = 0 constraint is handled. Let us first remark that global dual Schur complement S d is non-definite as soon as redundancies appear in the connectivity description of the interface, anyhow it is easy to prove [START_REF] Farhat | A method of finite tearing and interconnecting and its parallel solution algorithm[END_REF] that local contributions λ ⋄ b = A T λ b are unique (non-definition only affect the "artificial" splitting of forces on multiple points), and that because the right hand side lies in range(A) iterative process converges; other considerations on the splitting of physical efforts between subdomains will lead to improved initialization (see section 3.2.4 and [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF]).

Preconditioner to the dual interface problem

Like it is done in the primal approach, the most interesting preconditioners are researched as assembly of local contributions, and the global dual Schur complement being a sum of local contributions, optimal preconditioner is a scaled sum of local inverses.

S-1

d = ÃS ⋄ d + ÃT = ÃS ⋄ p ÃT (44)
Because this preconditioner uses local primal Schur complement, which corresponds to the local resolution of imposed displacement problems, it is commonly called the Dirichlet preconditioner. One interesting point is the possibility to give approximation of the local Schur complement operator leading to the following preconditioners:

S ⋄ p ≈ K ⋄ bb lumped preconditioner (45)
S ⋄ p ≈ diag(K ⋄ bb ) superlumped preconditioner (46)
These preconditioners have very low computational cost (they do not require the computation and storage of the inverse of local internal matrices

K ⋄ ii -1
), even if their numerical efficiency is not as strong as the Dirichlet preconditioner, they can lead to very reduced computational time.

Scaled assembly operator à can be defined the following way [START_REF] Klawonn | FETI and Neumann-Neumann iterative substructuring methods: connections and new results[END_REF]:

à = AM ⋄-1 A T + AM ⋄-1 (47)
where M ⋄ is the same parameter as for the primal approach. Such a definition is not that easy to implement, an almost equivalent strategy is then used, easily described using the (s) notation:

S-1 d = s M (s) A (s) S (s) p A (s) T M (s) (48)
• M (s) = diag( 1 multiplicity ) for homogeneous structures,

• M (s) = diag( diag(K (r) bb ) i j diag(K (j) bb ) i
) for compressible heterogeneous structures (assuming i represents the same degree of freedom shared by the j subdomains and (r) is the subdomain connected to (s) on degree of freedom i),

• M (s) = diag( µ (r) i j µ (j) i
) for incompressible heterogeneous structures (assuming i represents the same degree of freedom shared by the j subdomains and (r) is the subdomain connected to (s) on degree of freedom i).

We have the following partition of unity result:

ÃA T = I Υ ( 49
)
and the following complementarity between primal and dual scalings [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF]:

ÃT A + A T Ã = I ⋄ (50)
A (s) T M (s) A (s) + A (s) T M (s) A (s) = I Υ (s) (51)

Coarse problem

Admissibility condition G T λ b + e ⋄ = 0, can be handled with an initialization / projection algorithm (see section A.6): λ b = λ 0 + P λ * with G T λ 0 = -e ⋄ and G T P = 0.

λ 0 = -QG G T QG -1 e ⋄ (52)
P = I -QG G T QG -1 G T (53)
The easiest choice for operator Q is the identity matrix, projector P is then orthogonal, this choice is well suited to homogeneous structures. For heterogeneous structures, matrix Q has to provide information on the stiffness of subdomains, then Q is chosen to be a version of the preconditioner leading to "superlumped projector" (Q = Ã diag(K ⋄ bb ) ÃT ), "lumped projector" (Q = ÃK ⋄ bb ÃT ) and "Dirichlet projector" (Q = ÃS ⋄ p ÃT ). Superlumped projector is often a good compromise between numerical efficiency and computational cost. Algorithm 3.2 presents a classical implementation of FETI method, and figure 9 provides a schematic representation of the first iteration of the preconditioned dual approach.

Error estimate

The convergence of the dual domain decomposition method is strongly linked to physical considerations. After projection, the residual can be interpreted as the jump of displacement between substructures:

r = P T (-b d -S d λ) = Au ⋄ = ∆(u) (54)
∆(u) |Υ (i,j) = u (i) |Υ (i,j) -u (j) |Υ (i,j) (55)
Anyhow, such an interpretation cannot be connected to the global convergence of the system. In order to evaluate the global convergence, a unique interface displacement field has to be defined (most often using a scaled average of local displacement fields) and used to evaluate the global residual. When using the Dirichlet preconditioner, it is possible Algorithm 3.2 Dual Schur complement with conjugate gradient

1: Set P = I -QG(G T QG) -1 G T 2: Compute λ 0 = -QG(G T QG) -1 e 3: Compute r 0 = P T b d -S d λ 0 )
4: z 0 = P S-1 d r 0 set w 0 = z 0 5: for j = 0, . . . , m do 6:

p j = P T S d w j 7:
α j = (z j , r j )/(p j , w j ) 8:

λ j+1 = λ j + α j w j 9:
r j+1 = r jα j p j 10:

z j+1 = P S-1 d r j+1 11:
For 0 i j, β i j = -(z j+1 , p i )/(w i , p i )

12:

w j+1 = z j+1 + j i=1 β i j w i 13: end for 14: α ⋄ = (G T QG) -1 G t r m 15: u ⋄ = K ⋄+ λ m + R ⋄ α ⋄
to cheaply evaluate that convergence criterion. Average interface displacement u b can be defined as follow:

u b = A ÃT ∆ (56)
then, from equation [START_REF] Farhat | The two-level FETI method for static and dynamic plate problemspart I: An optimal iterative solver for biharmonic systems[END_REF], convergence criterion reads: Kuf = AS ⋄ p ÃT r . So when using the Dirichlet preconditioner, the evaluation of the global residual only requires the use of a geometric assembly after the local Dirichlet resolution.

Interpretation and improvement on the initialization

Let us come back to the original dual system [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF].

K ⋄ u ⋄ = f ⋄ + t ⋄T A T λ b At ⋄ u ⋄ = 0 (57)
And suppose this system is being initialized with non zero effort λ b0 :

K ⋄ u ⋄ = f ⋄ + t ⋄T A T λ b (58) let λ b = λ b + λ b0 K ⋄ u ⋄ = f ⋄ + t ⋄T A T λ b + t ⋄T A T λ b0 = f ⋄ + t ⋄T A T λ b (59)
with

f ⋄ = f ⋄ + t ⋄T A T λ b0 (60)
So initialization λ b0 can be interpreted as modification t ⋄T A T λ b0 of the intereffort between substructures: local problems are defined except for an equilibrated interface effort field; the only field that makes mechanical sense (and that is uniquely defined) is the assembly of interface efforts.

At ⋄ f ⋄ = At ⋄ f ⋄ = f b global interface effort (61) because At ⋄ t ⋄T A T λ b0 = AA T λ b0 = 0 (62)
Non-zero initialization then can be interpreted as a repartition of global interface effort f b . Two strategies can be defined in order to realize that splitting. Classical effort splitting Though splitting is hardly ever interpreted as a specific initialization, it is commonly realized, based on the difference of stiffness between neighboring substructures (that idea is strongly connected to the definition of scaled assembly operators): the aim is to guide the stress flow inside the stiffer substructure, sticking to what mechanically occurs.

Global interface effort f b is then split according to the stiffness scaling (M ⋄ = diag(K ⋄ bb )), which leads to modified local effort f ⋄ b .

f ⋄ b = M ⋄ A AM ⋄ A T -1 f b (63)
Complete effort f ⋄ is constituted by f ⋄ inside the substructure ((It ⋄T t ⋄ )f ⋄ ) and split effort on its interface (t ⋄T f ⋄ b ).

f ⋄ = (I -t ⋄T t ⋄ )f ⋄ + t ⋄T f ⋄ b ( 64
)
Because of the complementarity between scaled assembly operators (50), final effort reads

f ⋄ = f ⋄ -t ⋄T A T AM ⋄-1 A T + AM ⋄-1 t ⋄ f ⋄ (65)
Interface effort splitting [START_REF] Gosselet | On the initial estimate of interface forces in FETI methods[END_REF] If we start from condensed dual system (28)

S ⋄ p u ⋄ b = b ⋄ p + A T λ b Au ⋄ b = 0 (66)
condensed efforts can be split along the interface as long as global condensed effort remains unique. Assembled condensed interface effort reads b p = Ab ⋄ p , if it is split according to the stiffness of the substructures:

b ⋄ p = M ⋄ A AM ⋄ A T -1 b p (67)
We have, using the complementarity between scalings:

b ⋄ p = b ⋄ p -A T AM ⋄-1 A T + AM ⋄-1 b ⋄ p ( 68
)
Or in a non-condensed form:

f ⋄ = f ⋄ -t ⋄T A T AM ⋄-1 A T + AM ⋄-1 b ⋄ p ( 69
)
As will be shown in assessments, the classical splitting leads to almost no improvement of the method while the condensed splitting can be very efficient for heterogeneous structures. In fact the initialization associated to this splitting can be proved to be optimal in a mechanical sense; it can also be obtained from the assumptions used for the primal approach. Primal approach initialization is realized supposing that interface displacement field is zero on the condensed problem; from (66) we get:

A T λ b0 + b ⋄ p ≃ 0 (70)
which could only be the solution if null interface displacement was the solution. Then local interface efforts are split into an equilibrated part and its remaining ̺ ⋄ :

b ⋄ p = A T γ + ̺ ⋄ γ = AD ⋄ A T + AD ⋄ b ⋄ p ( 71
)
D ⋄ is a symmetric definite matrix, remaining ̺ ⋄ is orthogonal to range(D ⋄ A). If the system is initialized by:

λ b00 = -γ (72)
then initial residual A T λ b00 + b ⋄ p = -̺ ⋄ is minimal in the sense of the norm associated to

D ⋄ . If D ⋄ = diag(K bb ⋄ ) -1
then initialization is equivalent to the splitting of condensed efforts according to the stiffness of substructures ; diag(K bb ⋄ ) being an approximation of S p ⋄ that norm can be interpreted as an energy.

The initialization by the splitting of condensed efforts has to be made compatible with solid body motions by the computation of:

λ b0 = P λ b00 -QG G T QG -1 e ⋄ (73)
Remark 3.1 If D ⋄ = S ⋄ d was not computationally too expensive then improved initialization with Dirichlet preconditioner would lead to immediate convergence.

Remark 3.2

The recommended choice D ⋄ = diag(K bb ⋄ ) -1 is computationally very cheap, the heaviest operation is the computation of condensed efforts (one application of Dirichlet operator). Then if the Dirichlet preconditioner is used, new initialization is just as expensive as one preconditioning step but it can lead to significant reduction of iterations, so it should be employed. Of course if light preconditioner is preferred classical splitting should be used.

Three fields method / A-FETI method

The A-FETI method [START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: performance evaluation[END_REF][START_REF] Park | An algebraically partitioned FETI method for parallel structural analysis: algorithm description[END_REF] can be explained as the application of the three-field strategy [START_REF] Brezzi | A three-field domain decomposition method[END_REF] to conforming grids, this method is widely studied in [START_REF] Rixen | Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver[END_REF]. Back to (28), we have

S ⋄ p u ⋄ b = b ⋄ p + λ ⋄ b ( 74
)
Aλ ⋄ b = 0 (75) Au ⋄ b = 0 (76)
A-FETI method is based, like in the primal approach, on the introduction of unknown interface displacement field u b , the continuity of displacement then reads:

u ⋄ b = A T u b (77)
but local displacements are not eliminated like in the primal approach, complete system reads:

  S ⋄ p -I 0 -I 0 A T 0 A 0     u ⋄ b λ ⋄ b u b   =   b ⋄ p 0 0   (78)
In order to eliminate interface displacement u b a specific symmetric projector is introduced:

B = I -A T AA T -1 A (79)
B realizes the orthogonal projection on ker(A) (AB = 0). Since Aλ ⋄ b = 0 then λ ⋄ b can be written as

λ ⋄ b = Bµ ⋄ b (80)
µ ⋄ b is a new interface effort, corresponding (recall AA T = diag(multiplicity)) to an average of original λ ⋄ b . Introducing last result and using B T A T = 0 to eliminate interface displacement we have

S ⋄ p -B -B T 0 u ⋄ b µ ⋄ b = b ⋄ p 0 ( 81
)
Then using classical elimination of local displacement by the inversion of the first line of the previous system, we get

u ⋄ b = S ⋄ p + b ⋄ p + Bµ ⋄ b + R ⋄ b α ⋄ (82) R ⋄ b T b ⋄ p + Bµ ⋄ b = 0 (83)
which leads to

B T S ⋄ d B B T R ⋄ b R ⋄ b T B 0 µ ⋄ b α ⋄ = -B T S ⋄+ p b ⋄ p -R ⋄ b T b ⋄ p ( 84
)
This system is very similar to the classical dual approach system, and in consequence is solved in the same way (using projected algorithm). Anyhow the main difference is that Lagrange multiplier µ ⋄ b is defined locally on each subdomain and not globally on the interface.

It was proved in [START_REF] Rixen | Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver[END_REF] that A-FETI is mathematically equivalent to classical FETI with special choice of the Q matrix parameter of the rigid body motion projector. In fact if Q = diag( 1 multiplicity ) then FETI leads to the same iterates as A-FETI. Moreover operator B is an orthonormal projector which realizes the interface equilibrium of local reactions µ ⋄ b , it can be analyzed as an orthonormal assembly operator as described in figure 7.

To sum up, A-FETI can be viewed as the conforming grid version of the three-field approach, a specific case of classical FETI, and a dual approach with non-redundant description of the connectivity interface with orthonormal assembly operator.

Mixed domain decomposition method

Mixed approaches offer a rich framework for domain decomposition methods. It enables to give a strong mechanical sense to the method, mostly by providing a behavior to the interface. The mixed approach is one of the bases of the LaTIn method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Nouy | Une stratégie de calcul multiéchelle avec homogénéisation en temps et en espace pour le calcul de structures fortement hétérogènes[END_REF][START_REF] Ladevèze | On a mixed and multiscale domain decomposition method[END_REF], a very interesting strategy designed for nonlinear analysis; as we have restrained our paper to linearized problems, we do not go further inside this method which would deserve extended survey. Several studies were realized on mixed approaches, these methods possess strong similarities, we here mostly refer to works on so-called "FETI-2-fields" method [100,[START_REF] Series | Une méthode de décomposition de domaine avec deux multiplicateurs de Lagrange[END_REF].

The principle of the method is to rewrite the interface conditions:

Aλ ⋄ b = 0 Au ⋄ b = 0 ( 85
)
in terms of a new local interface unknown, which is a linear combination of interface effort and displacement.

µ ⋄ b = λ ⋄ b + T ⋄ b u ⋄ b ( 86
)
µ ⋄ b is homogeneous to an effort and T ⋄ b can be interpreted as an interface stiffness. Mixed methods thus enable to give a mechanical behavior to the interface, in our case (perfect interfaces) it can be mechanically interpreted as the insertion of springs to connect substructures. New interface condition reads:

A T Aλ ⋄ b + T ⋄ b A T Au ⋄ b = 0 ⋄ (87)
or

A T Aµ ⋄ b -A T AT ⋄ b -T ⋄ b A T A u ⋄ b = 0 ⋄ (88)
It is of course necessary to study the condition for this system being equivalent to system [START_REF] Rey | Reuse of krylov spaces in the solution of large-scale non linear elasticity problems[END_REF]. It is important to note that two conditions lying on the global interfaces (geometric and connectivity) were traced back to the local interfaces, so up to a zero-measure set (multiple points) the conditions have the same dimension. The new condition is equivalent to the former if facing local interfaces do not hold the same information which is the case

if matrix A T AT ⋄ b -T ⋄ b A T A is invertible.
An easy method to construct such matrices will be soon discussed If unknown µ ⋄ b is introduced inside local equilibrium equation, the local system reads:

S ⋄ p + T ⋄ b u ⋄ b = b ⋄ p + µ ⋄ b (89)
If we assume that T ⋄ b is chosen so that S ⋄ p + T ⋄ b is invertible then we have:

u ⋄ b = S ⋄ p + T ⋄ b -1 µ ⋄ b + b ⋄ p ( 90
)
Then substituting this expression inside interface condition (88), interface system reads:

A T A -A T AT ⋄ b -T ⋄ b A T A S ⋄ p + T ⋄ b -1 µ ⋄ b = A T AT ⋄ b -T ⋄ b A T A S ⋄ p + T ⋄ b -1 b ⋄ p (91)
so mixed approaches have the originality to rewrite global interface conditions on the local interfaces and to look for purely local unknown (which means that the size of the unknown is about twice the size of the unknown in classical primal or dual methods). This general scheme for mixed methods has, as far as we know, never been employed. A first reason is that it leads to certain programming complexity, second the manipulation of zero-measure interfaces is not easy for methods aiming at introducing strong mechanical sense and hard to justify from a mathematical point of view. So most often a simplified method is preferred, which takes only into account non-zero-measure interfaces. Such an approach simplifies the connectivity description of the interface, every relationship on the interface only deals with couples of subdomains. In order to have the clearer expression possible, we present the algorithm in the two subdomains case. Interface equilibrium reads:

u (1) b -u (2) b = 0 λ (1) b + λ (2) b = 0 (92) which is equivalent to    λ (1) b + λ (2) b + T (1) u (1) b -u (2) b = 0 λ (1) b + λ (2) b + T (2) u (2) b -u (1) b = 0 ( 93
)
under the condition of invertibility of T (1) + T (2) . Introducing unknown µ

(s) b = λ (s)
b + T (s) u (s) interface system reads:

µ (1) b + µ (2) b -T (1) + T (2) u (2) = 0 µ (1) b + µ (2) b -T (1) + T (2) u (1) = 0 (94) Local equilibrium reads:    S (1)
p + T (1) u

(1) b = µ (1) b + b (1) p S (2)
p + T (2) u (2) b = µ (2) b + b (2) p (95)
Assuming T (s) is chosen so that matrix S

p + T (2) -1

T (1) + T (2) S (1)
p + T (1) -1 I   µ (1) µ (2) =   T (1) + T (2) S (2)
p + T (2) -1 b (2) p T (1) + T (2) S (1)
p + T (1) -1 b (1) p   (96)
This expression enables to give better interpretation of the stiffness parameters T (s) . Suppose

T (1) = S (2)
p and

T (2) = S (1)
p then matrix ( 96) is equal to identity and solution is directly achieved. So the aim of matrix T (s) is to provide one substructure with the interface stiffness information of the other substructures.

If we generalize to N -subdomain system (91), we can deduce that the optimal choice for T (s) is the Schur complement of the remaining substructures on the interface of domain (s) (some kind of S (s) p where s denotes all the substructures but s). Of course such a choice is not computationally feasible (mostly because it does not respect the localization of data), and approximations have to be considered. In decreasing numerical efficiency and computational cost order, we have:

• Approximate the Schur complement of the remaining of the substructure by the Schur complement of the neighbor;

• approximate the Schur complement of the neighbor by the Schur complement of the nearer nodes of the neighbor ("strip"-preconditioners which idea is developed in [START_REF] Paz | An interface strip preconditioner for domain decomposition methods: application to hydrology[END_REF] in another context);

• approximate the Schur complement of the neighbor by the stiffness matrix of the interface of the neighbor (strategy of dual approach lumped preconditioner).

The second strategy is quite a good compromise: it respects data localization, it is not computationally too expensive and yet it enables the propagation of the information beyond the interface. Of course an important parameter is the definition of elements "near the interface", which can be realized giving an integer n representing the number of layers of elements over the interface.

Coarse problem

Because the interface stiffness parameter T ⋄ regularizes local operators S ⋄ p , local operator S ⋄ p + T ⋄ is always invertible. Such a property can be viewed as an advantage because it simplifies the implementation of the method introducing no kernel and generalized inverse; but it also can be considered as a disadvantage because no more coarse problem enables global transmission of data among the structure. Then the communications inducted by this method are always neighbor-to-neighbor which means that the transmission of a localized perturbation to far substructures is always a slow process. It is then necessary to add an optional coarse problem (see section A.7). Most often the optional coarse problem is constituted of would-be rigid body motions (if subdomains had not been regularized). Another possibility, which is proposed inside the LaTIn method is to use rigid body motions and extension modes of each interface as coarse problems, this leads to much larger coarse space. The coarse matrix corresponds to the virtual works of first order of deformation of substructures; so mechanically it realizes and propagates a numerical first order homogenization of the substructures.

Hybrid approach

The hybrid approach (see [START_REF] Gosselet | A monolithic strategy based on an hybrid domain decomposition method for multiphysic problems. application to poroelasticity[END_REF] for a specific application) is a proposition to provide a unifying scheme for primal and dual approaches though it could easily be extended to other strategies. It relies on the choice for each interface degree of freedom of its own treatment (for now primal or dual). So let us define two subsets of interface degrees of freedom: the first is submitted to primal conditions (subscript p) and the second to Neumann conditions (subscript d). Local equilibrium then reads (p

= i ∪ d, b = d ∪ p, p ∩ d = ∅): K ⋄ pp K ⋄ pp K ⋄ pp K ⋄ pp u ⋄ p u ⋄ p = f ⋄ p f ⋄ p + t ⋄ d T λ ⋄ d λ ⋄ p ( 97
)
Preferred interface unknowns are unique displacement on the first subset u p and unique effort on the second subset λ d . Local contributions then reads:

u ⋄ p = A T p u p ( 98
)
λ ⋄ d = A T d λ d (99)
which ensure the continuity of displacement on the p degrees of freedom and the actionreaction principle on the d degrees of freedom, of course operators A p and A d have been restricted respectfully to the p and d subsets. Remaining interface conditions read:

A d u ⋄ d = A d t ⋄ d u ⋄ p = 0 (100) A p λ ⋄ p = 0 (101)
To obtain the global interface system, first local unknown u ⋄ p has to be eliminated:

u ⋄ p = K ⋄ pp + f ⋄ p + t ⋄ d T λ ⋄ d -K ⋄ pp u ⋄ p + R ⋄ pα ⋄ (102)
Applying continuity condition to previous result and equilibrium condition to the second row of (97), interface system reads:

  S pd G p G d -G T p G T d 0     u p λ d α ⋄   =   b p -b d -e ⋄   (103)
with the following notations:

S pd = A p 0 0 A d S ⋄ pd A p 0 0 A d T G p = A p K pp ⋄ R ⋄ p, G d = A d t ⋄ d R ⋄ p e ⋄ = R ⋄ pT f ⋄ p b p = A p f ⋄ p -K ⋄ pp K ⋄ pp + f ⋄ p , b d = A d t ⋄ d K ⋄ pp + f ⋄ p
This interface problem corresponds to the constrained resolution of one linear system. The constraint is linked to the possible non-invertibility of matrix K ⋄ pp and thus to the choice of primal subset. Notice that G p represents the reaction of primal degrees of freedom to zero energy modes of K ⋄ pp and then should be zero in most cases (it may be non zero in buckling cases). Moreover this system may represent classical primal approach (if all interface degrees of freedom are in subset p) or classical dual approach (if all interface degrees of freedom are in subset d). Operator S pd is a primal/dual Schur complement, it is the sum of local contributions S ⋄ pd (25). The above system is nonsymmetric semi-definite (because of redundancies on the dual subset) positive, it has to be solved by GMRes-like algorithm.

Hybrid preconditioner

Inspired by primal and dual preconditioners, we propose to approximate the inverse of the sum of local contributions by a scaled sum of local inverses.

S-1

pd = Ãp 0 0 Ãd S ⋄ dp Ãp 0 0 Ãd T ( 104
)
Scaled assembly operator are defined in the same way as in primal and dual approaches.

Coarse problems

As said earlier, depending on the choice of subset p, local operator K ⋄ pp (involved in the computation of S ⋄ pd ) may not be invertible and, like in dual approach, a first coarse correction has been incorporated inside the hybrid formulation. Anyhow local operator K ⋄ d d involved in preconditioning step may also not be invertible and, like in primal approach, a second coarse problem has to be added to make the preconditioner optimal. Then, the optimal version of the hybrid system incorporates two coarse problems handled by specific initialization/projection algorithm. The admissibility condition reads:

G T x = e (105) with G = G p G d , x = u p λ d and b = b p -b d . If r = b -S pd
x stands for the residual before preconditioning, the optimality condition reads:

H T r = H p H d T r p r d = 0 ( 106
)
with

H p = Ãp t ⋄ p R ⋄ d and H d = Ãd K ⋄ d dR ⋄
d (as said before most often H d = 0). To sum up constraints:

G T x = -e ⋄ H T S pd x = H T b (107)
Handling such constraints is described in section A.8.

Figure 10 provides a schematic representation of the first iteration of the preconditioned hybrid approach, in the specific case of a nodal partition of the interface. Assessments will deal with partition at the degree of freedom level.

Error estimate

Because GMRes-like solver is used, the euclidian norm of the residual is directly available, such a norm is the sum of displacement gap on the d part of the interface and an effort gap on the p part of the interface. For now, no other estimator with better physical sense is available.

Adding optional constraints

The aim of optional constraints is to ensure that the research space for the iterative solver possesses a certain regularity. The choice of constraints will be discussed in a later section. Going back to the interface system:

  S ⋄ p -I ⋄ A 0 0 A   u ⋄ b λ ⋄ b =   b ⋄ p 0 0   (108)
any constraint of the form C T Au ⋄ b = 0 or C T Aλ ⋄ b = 0 is trivially verified by the solution fields, it is just a restriction of continuity/equilibrium conditions. From an iterative process point of view, these conditions will be reached once converged ; the principle of optional constraints is to have every iteration verify that condition.

There are two classical solutions to ensure these optional constraints: either to realize a splitting of research space and ensure, using a projector, that the resolution is limited to convenient subspace, or to realize a condensation of constraints and make iterations on a smaller space. In other words suppose there are n c independent constraints in a ndimension space the first strategy researches n-sized solution in a (nn c )-ranked space, while the second solution researches (n-n c )-sized solution in a (n-n c )-dimension space then deduce the n-sized solution. From a numerical point of view both solutions are equivalent, they are just two ways of handling the same constraints, anyhow from implementation and computational points of view they have strong differences.

We will essentially focus on application to primal and dual domain decomposition methods.

Augmentation strategy

For this strategy the constraint is reinterpreted in terms of constraint on the residual. Typically, the primal approach can be augmented by constraints on the effort field while the dual approach can have constraints on the displacement field:

C T Aλ ⋄ b = -C T Ab ⋄ p -AS ⋄ p A T u b (109) C T Au ⋄ b = C T AS ⋄ d A T λ b + b ⋄ d ( 110
)
The constraint is then handled as an augmentation inside the iterative solver, its is classically realized using a projector (see sections A.7 and A.8).

Recondensation strategy

This strategy was recently introduced in the framework of the dual approach, leading to the FETIDP algorithm [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part i: a faster alternative to the two-level FETI method[END_REF][START_REF] Farhat | A scalable dual-primal domain decomposition method[END_REF][START_REF] Lesoinne | FETI-DP: An efficient, scalable, and unified Dual-Primal FETI method[END_REF][START_REF] Klawonn | Some computational results for dual-primal feti methods for three dimensional elliptic problems[END_REF]. Because for now only constraints on the u ⋄ b field have been considered we will restrain to this kind of constraints, the application of constraints on λ ⋄ b is straightforward.

Basic method

We first consider constraints which impose continuity of specific degrees of freedom, in other terms we suppose matrix C is identity on certain degrees of freedom and zero elsewhere; we will show how any constraint can be rewritten in such a form. Because these degrees of freedom will be submitted to a primal treatment we will denote them with subscript p while the remaining of the interface will be denoted with subscript d. Constraint reads:

0 = C T Au ⋄ b = 0 I p T A d u ⋄ d A p u ⋄ p = A p u ⋄ p (111)
Like in the hybrid approach this constraint is verified using a unique displacement field on the primal part of the interface: u ⋄ p = A T p u p . Interface system then reads like in the hybrid approach, with the additional assumption that the constraints are so that the local problem possesses enough Dirichlet conditions to make it invertible.

S pd u p λ d = b p -b d (112)
Introducing following notations for blocks composing S ⋄ pd :

S ⋄ pd = s ⋄ pp s ⋄ pd -s ⋄ dp s ⋄ dd (113)
then

S pd = s pp s pd -s dp s dd = A p 0 0 A d s ⋄ pp s ⋄ pd -s ⋄ dp s ⋄ dd A p 0 0 A d T (114)
Unknown u p is condensed on the remaining of the interface:

u p = s -1 pp b p -s pd λ d ( 115
)
s d λ d = s dd + s dp s -1 pp s pd λ d = -b d + s dp s -1 pp b p (116)
The latest equation is solved using an iterative solver, operator s d has the same properties as the restriction of the dual operator to the d-part of the interface (semi-definition, symmetry and positivity). Operator s d is the assembly of local contributions

s d = A d s ⋄ d A T d = A d s ⋄ dd + s ⋄ dp A T p s -1 pp A p s ⋄ pd A T d ( 117
)
Using operator s d requires the computation of the inverse of matrix s pp = A p s ⋄ pp A T p , which is an assembled matrix. Then this formulation includes a global coarse problem set on primal variables.

The recommended preconditioner for such an approach is directly inspired by the dual approach: it consists in solving local Dirichlet problems with scaled imposed displacement on the d-part of the interface and null displacement on the p part of the interface and extracting the average reaction of the d-part of the interface. Then the preconditioner reads:

s-1

d = 0 p à S ⋄ p 0 p à T (118)
Figure 11 provides schematic representation of the first iteration of preconditioned FETID method.

More complex constraints

We now consider constraints which are not limited to one degree of freedom, for instance one can consider that we want to ensure that the average jump of displacement on one edge is equal to zero, which involves all the degrees of freedom of the edge.

The classical solution [START_REF] Klawonn | Some computational results for dual-primal feti methods for three dimensional elliptic problems[END_REF] is to realize a change of basis of degrees of freedom (denoted by matrix T ⋄ ) so that each constraint is represented by one "modified" degree of freedom. The change of basis is the same local operation realized on every subdomain, then we can define a global change T so that

AT ⋄ = T A C T Au ⋄ b = C T AT ⋄ û⋄ b = C T T Aû ⋄ b = ĈT Aû ⋄ b ( 119
) with Ĉ = 0 I p (120)
After the change of basis is realized the same algorithm can be employed. Because constraints most often respect a certain locality of data (for instance independent constraints on each edge), change of basis is not a too expensive operation, and does not make too poor the sparsity of local matrices.

Adding "constraints" to the preconditioner

This subsection deals with the dualization of the recondensation strategies [START_REF] Cros | A preconditioner for the schur complement domain decomposition method[END_REF][START_REF] Dohrmann | A preconditioner for substructuring based on constrained energy minimization[END_REF]. For instance the balanced domain decomposition with constraints (BDDC) is a primal version of FETIDP: during the preconditioning step, the continuity of displacement is ensured at specific degrees of freedom (which can be the result of a local change of basis), so that the local Neumann operator remains fully invertible. So solving classical primal approach problem

S p u b = AS ⋄ p A T u b = b p (121)
the preconditioner reads

S-1 p = I p 0 0 Ãd s -1 pp -s -1 pp A p s ⋄ pd -s ⋄ dp A T p s -1 pp s ⋄ d I p 0 0 Ãd T ( 122
)
Figure 11 provides a schematic representation of the first iteration of preconditioned BDDC method.

Classical issues

In this section we try to provide answers to questions that commonly arise when using domain decomposition methods.

Rigid body motion detection

Handling floating substructures is definitely a very accurate issue. This difficulty is one of the reason of the success of methods which regularize the stiffness of subdomains like FETIDP or mixed approaches, leading to fully invertible matrices. Anyhow basic primal and dual methods remain competitive (mostly because zero-energy modes provide a very natural coarse problem), hence providing an efficient algorithm for the computation of rigid body motion is essential. Many strategies can be used [START_REF] Farhat | On the computation of the null space and generalized invers of large matrix, and the zero energy modes of a structure[END_REF][START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF], and this review does not claim to be exhaustive.

First we have to discuss the exact composition of the possible kernel of substructures. What can be found is:

• rigid body motions of floating substructures,

• internal mechanisms of substructures,

• weird things due to nonlinearities (buckling, exotic behaviors, ...)

• numerical zero-energy modes.

An internal mechanism can for instance occur when a substructure is made out of two parts connected by one linear edge (pivot) or one singular point (kneecap). Methods exist to avoid such substructures either inside the decomposition algorithm or as external programs regularizing a given decomposition [START_REF] Fragakis | A unified framework for formulating domain decomposition methods in structural mechanics[END_REF], and then should be employed.

The last two kinds of kernel are non-standard and can only be detected using fully algebraic methods (like Gauss pivoting). The problem with algebraic methods is their high sensitivity to the condition number of the stiffness matrix of substructures. The condition number can be influenced by the aspect ratio and the material composition of the substructures, even after adimensionalization the quality of the methods is hard to warranty.

Finally we only develop here two strategies to handle zero-energy modes. The first one belongs to the purely algebraic methods, it is very simple to implement and can lead to satisfying results for not-too-complex problems, it can handle more than solid body motions but it is strongly dependent on a priori selected degrees of freedom. The second method is purely geometric, it is very robust but only suited to detect solid body motions.

In order to simplify notations, we consider the research of the zero-energy modes of matrix K (which should be a local stiffness matrix).

Simple algebraic approach

This method is based on fundamental relationship between the kernel of a matrix and the kernel of Schur complement [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF]. The principle is to preselect a small set of degrees of freedom which we will denote by subscript N (the other degrees of freedom are denoted with subscript O). Then compute explicitly primal Schur complement associated to these degrees of freedom:

S = K N N -K N O K -1 OO K ON ).
If N -degrees of freedom are selected so that K 00 is invertible (if only solid body motions have to be detected then it is sufficient to take the degrees of freedom associated to three non-aligned nodes) then S is well defined and its kernel is linked to the kernel of matrix K.

Since N is a "small" set (12 degrees of freedom is often sufficient), computing the kernel of matrix S using "exact" algorithm like singular value decomposition is rather cheap and then kernel of matrix K can be computed using equation [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF].

Geometric approach

The basic idea of this method is that kinematically admissible solid body motions can be deduced from boundary conditions imposed on one subdomain. Let R c be a basis of candidate rigid body motions (would be solid body motions if no Dirichlet conditions were applied on the subdomain, R c is a 6-column matrix in 3D and 3-column matrix in 2D), and let E be the matrix of Dirichlet boundary conditions: each column of E represents one (combination of) blocked degree of freedom. Kinematically admissible rigid body motions R are linear combinations of candidate rigid body motions (hence R = R c Q) which do not make Dirichlet boundary conditions work (ie E T R = 0). In order to find such linear combination, we compute singular value decomposition of matrix E T R c = U DV T and set Q = V 0 where V 0 is the submatrix of V associated to negligible singular values. Because E T R c is a matrix only made out of geometric considerations, it is well conditioned and the criterion to detect "zero" singular values is well defined (there is a large gap between zero and non-zero singular values).

Generalized inverse

The two methods presented above led to r-ranked basis R of the kernel of matrix K. To compute generalized inverse K + , the most classical way is to select r degrees of freedom, so that if they were added Dirichlet conditions, rigidity matrix would be well defined. For instance after pivoting and renumber one could get R = R * I r and then the r last degrees of freedom would suit. Then we have:

K = K rr K rr K rr K rr and K + = K -1 rr 0 0 0 ( 123
)
This is just one instance of generalized inverse, other can be built using penalization or other modification to matrix K. Though not theoretically prohibited, choosing "blocked" degrees of freedom on the interface is often a bad idea from a practical point of view.

Choice of optional constraints

As seen in section 4, constraints can be imposed either using augmentation (using one or two projectors, sections A.7 and A.8) or using recondensation algorithms. In the case of recondensation algorithms, constraints have to be sufficient in order to suppress rigid body motions and then regularize the local stiffness matrix. In the case of augmentation algorithm, constraints have to be independent from rigid body motions which are already handled by the formulation. Because constraints are expensive to handle, they have to be chosen with care; anyhow, except in a few cases, there are no general results on how to choose them. Mechanical comprehension of studied phenomena and anticipation of convergence difficulties may lead to efficient strategies. In the case of solving several linear systems (even with different left hand sides) interesting strategies exist [START_REF] Saad | On the Lanczos method for solving symmetric linear systems with several right hand sides[END_REF][START_REF] Erhel | An augmented conjugate gradient method for solving consecutive symmetric positive definite linear systems[END_REF][START_REF] Gosselet | On a selective reuse of krylov subspaces in newton-krylov approaches for nonlinear elasticity[END_REF][START_REF] Rey | Solution to large nonlinear systems: acceleration strategies based on domain decomposition and reuse of krylov subspaces[END_REF][START_REF] Risler | On the reuse of Ritz vectors for the solution to nonlinear elasticity problems by domain decomposition methods[END_REF].

The next two subsections deal with very common strategies, while the last subsection describes another framework for constraints inspired by the LaTIn method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Nouy | Une stratégie de calcul multiéchelle avec homogénéisation en temps et en espace pour le calcul de structures fortement hétérogènes[END_REF].

Forth order elasticity

As plate and shell models are often used in structural mechanics, forth order problems have been carefully studied [START_REF] Farhat | The two-level FETI method for static and dynamic plate problemspart I: An optimal iterative solver for biharmonic systems[END_REF][START_REF] Farhat | The two-level FETI method -part II: Extension to shell problems. parallel implementation and performance results[END_REF][START_REF] Roux | Parallel implementation of direct solution strategies for the coarse grig solvers in 2-level FETI method[END_REF]103,[START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF]. Such problems are characterized by the appearance of singularities at the corner of substructures (so-called "corner modes") which are destroying the scalability of usual methods. The classical solution consists in enforcing the continuity of the (most often only normal) displacement field at the corners in order to regularize the problem. From a practical point of view, corners are most often defined as multiple points (nodes shared by more than two substructures), that set can be enriched by extremities of edges.

The implementation, just like the singularity, is strongly linked to the chosen formulation.

Dual approach: since the projected residual corresponds to the displacement jump between substructures [START_REF] Glowinski | Augmented lagrangian interpretation of the nonoverlapping Schwartz alternating method[END_REF], one just has to use augmentation algorithm with one constraint for each pair of neighbor at each corner. Matrix C is then made out of columns with one coefficient 1 on one corner degree of freedom and 0 elsewhere. Because the dual description of a m-multiple point leads to (m -1) relationships, such a coarse space is rather large.

Primal approach: In order to regularize the displacement field, the constraints have to be imposed on the preconditioned residual (assuming Neumann-Neumann preconditioner is employed). The aim is then to have the local contributions of preconditioned residual equal to zero on corner points. So if C ⋄ denotes the local interface matrix made out of columns with one coefficient 1 on one corner degree of freedom and 0 elsewhere, and C⋄ the same matrix scaled according to the scaling used inside the preconditioner, the constraints read

C = ÃS ⋄ d C⋄ .
Then a m-multiple degree of freedom leads to m constraints.

Recondensed approaches: FETIDP or BDDC were first designed in this context, from the consideration that (extended) corners constraints were sufficient to suppress rigid body motions, then the first level constraints could be avoided. So the methods directly apply since they consist in constraining the displacement field. Here whatever the multiplicity of a corner may be, it always leads to just 1 constraint.

Second order elasticity

Because classical methods are already scalable in the frame of second order elasticity, optional constraints are not often used in such a context. Furthermore, it is hard to predict what constraint should be imposed. In some cases, efficient strategies have been proposed, such as in [START_REF] Gosselet | On a selective reuse of krylov subspaces in newton-krylov approaches for nonlinear elasticity[END_REF] for nonlinear problems using Newton-Raphson solver where approximations of eigen vectors are used. The question of optional constraints arose when willing to extend recondensed algorithms (FETIDP and BDDC) to such problems, mostly because the previous definition of corners lead to significant problems in 3D (too many constraints, poor convergence...). The first solution was proposed in [START_REF] Lesoinne | FETI-DP: An efficient, scalable, and unified Dual-Primal FETI method[END_REF], the idea was to select 3 non-aligned nodes on each face (interface between 2 subdomains) which maximized the surface of the triangle they defined. The current solution, the scalability of which is mathematically and numerically proved, is to enforce average convergence on edges [START_REF] Klawonn | Some computational results for dual-primal feti methods for three dimensional elliptic problems[END_REF], which is realized by a change of basis described in (119). In order to take into account heterogeneity on the interface, the average may be scaled by a coefficient representing the stiffness of the subdomains. For more difficult problems, first order moments of edges can also be added.

Link with homogenization theory

This paragraph intends to present a mechanical vision of optional constraint which, though hard to implement in the framework of the presented method, may lead to better understanding of what optional constraints and associated coarse problems can provide to the methods. This analysis is inspired by the multiscale version of the LaTIn method [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics -New Approaches and Non-Incremental Methods of Calculation[END_REF].

The underlying question when choosing optional constraints (except from specific numerical questions like in the forth order elasticity) is "what global information should be transmitted to the whole structure ?" or more precisely "what should far substructures know from one substructure". A meaningful answer is provided by Saint-Venant principle and homogenization theory: at a first order development, a substructure can be represented by its rigid body motions and its constant strain states (simple traction and shearing states). Such an idea adds six (3 in 2D) more constraints per subdomains; as these constraints are somehow complex to build, they can be approximated by interface modes (but of course the number of constraints then grows quickly).

Linear Multiple Points Constrains

Multiple points constraints are relationships defined between some degrees of freedom, they are often used in order to connect nonconforming meshes, to represent boundary conditions (for instance periodicity), to model contact or apply control laws. In the case of linear(ized) constraints we can write, on the whole structure scale:

Ku = f (124) Cu = a (125)
What seems most suited to the domain decomposition context is to dualize the constraint and introduce Lagrange multiplier µ in order to enforce the condition. System then reads:

K C T C 0 u µ = f a (126)
After decomposition we have

   K ⋄ -I ⋄ C T At ⋄ 0 0 0 At ⋄ 0 C 0 0      u ⋄ λ ⋄ µ   =    f ⋄ 0 0 a    (127)
with C so that Cu ⋄ = Cu = a which implies (since

u ⋄ b = A T u b ): Cu = C i C b u ⋄ i u b = Cu ⋄ = C i C b u ⋄ i u ⋄ b (128) then C b A T = C b (129)
Or in other words, if matrix C deals with interface degrees of freedom, the associated constraints have to be correctly distributed between sharing subdomains. The constraint can be interpreted as specific (non-boolean) assembly operator which explains the chosen notation. Using MPCs with domain decomposition methods was studied in [START_REF] Daniel Rixen | Extended preconditioners for the feti method applied to constrained problems[END_REF] in the dual method context.

In order to provide general methodology to apply MPCs, we then incorporate inside hybrid domain decomposition method (equations ( 97) to (100)):

  K ⋄ pp K ⋄ pp A T p C T p A p K ⋄ pp A p K ⋄ pp A T p C T p C p C p 0     u ⋄ p u p µ   =   f ⋄ p + t ⋄ d T A T d λ d A p f ⋄ p a   (130)
A d t ⋄ d u ⋄ p = 0 (131)
The elimination of u ⋄ p leads to, with classical hybrid notations:

    S pd C T p C T d G p G d -C p C d C µ G α -G T p G T d G T α 0        u p λ d µ α ⋄    =    b p -b d h -e ⋄    (132)
with

C p = C p -C pK ⋄ pp + K ⋄ pp A T p (133)
C d = C pK ⋄ pp + t ⋄T d A T d (134)
C µ = C pK ⋄ pp + C T p (135)
G α = C pR p (136) h = a -C pK ⋄ pp + f ⋄ p (137)
Various strategies can be used in order to solve system (132), which combine elimination of constraints (rigid body motions and/or MPCs) by projection methods (FETI-like approaches) and/or by recondensation methods (FETIDP like approaches). All these methods correspond to solving rather complex coarse problems but they are suited to traditional preconditioners. We propose, after [START_REF] Daniel Rixen | Extended preconditioners for the feti method applied to constrained problems[END_REF], to use classical projection method to handle rigid body motions then use iterative solver to find simultaneously (u p , λ d , µ) and provide efficient preconditioner to this problem.

System reads with trivial notations (for simplicity reasons we suppose that the rigid body motion constraints have been symmetrized, which is always possible and which is naturally the case if G p = 0 like in many applications):

S pdµ G G T 0 x α ⋄ = b -e ⋄ (138)
with

S pdµ =   A p 0 0 A d t ⋄ d 0 C p   S ⋄ pd   A p 0 0 A d t ⋄ d 0 C p   T +   0 0 C p 0 0 0 -C T p 0 0   (139)
As can be seen, MPCs have very different actions wether they are set on primal interface degrees of freedom or not: matrix C p modifies the structure of the hybrid system while dual and internal constraints lead to classical hybrid approach with modified dual trace and assembly operator

A = A d t ⋄ d C p
. The definition of efficient preconditioner inspired by classical methods is then much simplified if no constraints are set on primal degrees of freedom (C p = 0), which is what we suppose now:

S pdµ = A p 0 0 A S ⋄ pd A p 0 0 A T (140)
The S ⋄ pd notation is due to the association of the trace operator with the assembly operator, which in fact is equivalent to defining "extended interface dual degrees of freedom" made out of dual degrees of freedom and degrees of freedom involved in MPCs. Since, in this hypothesis, system reads like classical hybrid approach, that is (modified) assembly of local contributions, the proposed preconditioner is a scaled assembly of local inverses.

S-1

pdµ = Ãp 0 0 Ã S ⋄ dp Ãp 0 0 Ã T (141)
The primal scaled assembly operator can be directly imported from the primal approach. Concerning the dual approach, according to previous definitions, we have:

à = AM ⋄ p -1 A T + A ⋄ M ⋄ p -1 (142)
Where M ⋄ p is a diagonal matrix chosen like in the classical methods. The matrix to inverse reads:

AM ⋄ p -1 A T = A d t ⋄ d M ⋄ p -1 t ⋄ d T A T d A d t ⋄ d M ⋄ p -1 C T p C pM ⋄ p -1 t ⋄ d T A T d C pM ⋄ p -1 C T p ( 143
)
The idea is then to make this system easy to inverse, having the off-diagonal blocks equal to zero. We have

C p = C i C d with C d A T d = C d ; if we choose C d = C d A d t ⋄ d M ⋄ p t ⋄ d T A T d -1 A d t ⋄ d M ⋄ p t ⋄ d T (144)
then

A d t ⋄ d M ⋄ p -1 C
T p = 0 and the scaled assembly operator is non expensive to compute. In others words, one simply has to split constraints on interface degrees of freedom between subdomains according to the scaling used inside the preconditioner.

Choice of decomposition

Decomposing a given structure in order to use the algorithms presented in this paper is a complex problem. Algorithms and softwares were proposed [START_REF] Farhat | Top/domdec -a software tool for mesh partitioning and parallel processing[END_REF][START_REF] Cross | A parallelisable algorithm for optimising unstructured mesh partitions[END_REF][START_REF] Karypis | Multilevel algorithms for multi-constraint graph partitioning[END_REF], which mostly refer to graph theory. Such approaches enable to take into account load balance between processors (supposing each processor is assigned to one subdomain, this is equivalent to making each local problem as easy to solve as others) and to minimize the dimension of the interface (so that the global condensed problem is as small as possible). They also enable to avoid internal rigid body motions (so called mechanisms).

Anyhow another important point is to have local problems as well conditioned as possible, so having subdomains with good aspect ratio (ratio between largest and smallest characteristical dimensions of the subdomain) is considered as an important point. Anyhow one has to realize that good aspect ratio is often linked to large local bandwidth and then to somehow expensive local problems to solve.

An even more difficult point to take into account high heterogeneities (see figure 13): using stiffness scaling enables to correctly handle heterogeneity when interface between subdomains matches interface between materials, anyhow when interface between subdomains "crosses" interface between materials then numerical difficulties may occur. For now scaled-average optional constraints seem to be the best solution to handle these difficulties but it leads to large coarse problems.

Finally finding the best decomposition is still a rather open problem and mechanical sense is often a necessary complement to efficient automatic decomposing algorithm. Nonsymmetry occurs in many physical modeling: plasticity, nonlocal models for fracture [START_REF] Germain | Méthodes de calcul non local : Application aux structures composites[END_REF], frictional contact [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multicontact problems[END_REF]. The use of the domain decomposition methods presented in this paper just requires more care in the implementation because some simplifications are not available (for instance coarse problem matrices are nonsymmetric), and of course the use of well suited iterative solver like GMRes, OrthoDir or BiCG because Schur complements are no longer symmetric. Globally methods show good numerical performance results. Anyhow a real problem is the absence of theoretical results to ensure good convergence properties (this is mainly due to the fact that proofs for classical methods rely on the construction of an interface inner-product related to Schur complements which is no longer possible).

Nonlinear problems

For now we have considered the solution to linear systems. To adapt the method to nonlinear problems, a classical solution is to use Newton-Raphson linearization scheme: linearized stiffness matrices are computed independently on each subdomain then the linearized system is solved using domain decomposition [START_REF] Rey | Reuse of krylov spaces in the solution of large-scale non linear elasticity problems[END_REF][START_REF] Gosselet | On a selective reuse of krylov subspaces in newton-krylov approaches for nonlinear elasticity[END_REF][START_REF] Risler | On the reuse of Ritz vectors for the solution to nonlinear elasticity problems by domain decomposition methods[END_REF][START_REF] Farhat | The second generation FETI methods and their application to the parallel solution of large-scale linear and geometrically non-linear structural analysis problems[END_REF][START_REF] Gosselet | A domain decomposition method for quasi incompressible formulations with discontinuous pressure field[END_REF]]. For such approaches, domain decomposition methods can be seen as efficient black-boxed linear solvers.

One critical point in such a method is that, depending on the formulation (for instance fully Lagragian formulation of a large deformation elasticity problem) rigid body motions may vary from one system to the other. In the proposed context, what has been proved is that translations always belong to zero energy modes, what has been observed is that rotations only appeared as zero-energy-modes in the first system (which corresponds to linear elasticity problem). This might be penalizing because the size of the information transmitted inside the coarse problem is decreased after the first system; moreover, rotations are often converted to "negative-energy modes" which, if they are in small number, can be handled by fully-reorthogonalized conjugate gradient (though convergence will be slower). One classical solution is to reinject disappeared rotations as optional constraints (via augmentation algorithms).

In the case where nonlinearity is localized in few substructures, an interesting strategy can be to carry out subiterations of the nonlinear solver independently in those substructures [START_REF] Cai | Nonlinearly preconditioned inexact newton algorithms[END_REF][START_REF] Cresta | Comparison of multiscale and parallel nonlinear strategies based on domain decomposition for post buckling analysis[END_REF].

Implementation issues

Implementing domain decomposition methods from existing code is not a too complex task. We give a few details on our software architecture though practical solutions are far from being unique. Our code is a plug-in to ZeBuLoN object-oriented finite element computational software [START_REF]Z-set user manual[END_REF][START_REF]Z-set developper manual[END_REF], it takes advantage of Frederic Feyel's previous work [START_REF] Feyel | Application du calcul parallèle aux modèles à grand nombre de variables internes[END_REF][START_REF] Feyel | Quelques multi-problèmes en mécanique des matériaux et structures[END_REF]. Our implementation aims at being as generic as possible, so for now hybrid domain decomposition method has been developed (and mixed approaches are under construction), and separation between formulation and solver (so that any iterative solver can be used to solve the interface problem). All classical projectors and preconditioners are implemented.

The most basic pieces of the code are:

• topological description of the interface, ie ability to realize trace operations;

• exchange library (PVM, MPI), ie ability to realize assembly operations;

• classical FE code, ie ability on any given subdomain to get any local field given sufficient boundary conditions.

Organization of the topological information

In order to implement hybrid approach, we propose to define a specific class for "interface degrees of freedom" which wraps classical degrees of freedom and provide information on:

• the number of subdomains which share this degree of freedom

• the kind of treatment (primal, dual...) which is applied Then degrees of freedom are set together neighborwise, defining "interface" class made out of a list of pointer to "interface degrees of freedom", and the global number of the neighbor (that number enables to identify the subdomain and to realize exchanges).

A "subdomain" is defined as a collection of "interfaces" and a classical domain in the sense of usual FE software (mainly mesh), it possesses its own global identification number.

Note that with such a description of the decomposition, the local interfaces are redundant, multiple degrees of freedom appear in several interfaces. It is then necessary to take certain care to define some operations (transposed trace for dual degrees of freedom). A specific class can then be used to ease the management of multiple degrees of freedom. In order to easily connect the domain decomposition formulation to an iterative solver, we propose to define "interface vectors" (displacement, intereffort...), "interface matrices" (trace of rigid body motions..., can be seen as a collection of interface vectors) and "interface operators" (square interface matrices), with all classical operations (basically sum, difference, product and transposed product).

The particularity of these objects is to be defined on the interface and then data is shared between subdomains, so all the previous operations sometimes require to assemble data (an interesting idea is to have a boolean member indicating the assembled state of data). Because of the choice of description of the interface, the assembly operation requires certain care for primal multiple degrees of freedom (these degrees of freedom shall have the same value at their different occurrences). Note that an object like "interface matrix" can highly be optimized (mainly in terms of memory storage).

Interface operators mainly know how to multiple vectors and matrices, they are used to define Schur complements, scaling operators, projectors. Composite design pattern can be used to simplify the succession of operations.

In order to let the user choose the various configurations of the domain decomposition method, we use inheritance and object factory design pattern. 5.6.3 Articulation between formulation and solver (fig. 16) What we propose is to have client/server relationship between solver and formulation: basically an iterative solver needs to know how to initialize, how to multiply, how to precondition, how to do inner products, how to evaluate convergence. All these operations are implemented inside the "interface formulation" object which is linked to one subdomain (topology and stiffness) and creates "interface algebraic objects" in order to define required operations.

Assessments

The assessments we propose here are based on the code described in the previous section. Basically, we have implemented the hybrid approach which lets us assess the classical primal and dual approaches with most classical preconditioners and projectors.

We first present a sequence of academical tests in order to recover classical numerical performance results (scalability and relative efficiency of the different approaches): two dimensional plane stress problem, plate problem, three dimensional problem with heterogeneity and unstructured decomposition. In these problems, H denotes the characteristic size of the subdomains and h the characteristic size of the elements. We also present results on non-academic problem (bitraction test specimen).

For all those tests, in order to compare all the approaches (including the hybrid approach), GMRes solver is used and convergence is monitored by the norm of the residual as given by the solver (with ε set to 10 -6 ). In other cases (when hybrid approach is not assessed), the convergence is monitored through the evaluation of global primal residual |Kuf |/|f | < ε with ε set to 10 -6 . Depending on the method, a different coarse problem may be introduced, we denote by CS:a+b the size of the coarse problems (a for the admissibility coarse problem and b for the optimal preconditioning coarse problem) or number of iterations total number of constraints . Note that the hybrid approach deals with two independent coarse problems, so their solutions is much cheaper than the solution to a unique large coarse problem.

We test the primal approach with or without optimality coarse problem, the dual approach with lumped or Dirichlet preconditioner and identity or superlumped or Dirichlet projector (denoted by P(I), P(W) and P(D) respectively). As for such tests no physical consideration can guide the choice of hybrid treatment to the interface, in order to show the potential of the hybrid approach, we present results where all degrees of freedom of one direction (U 1 , U 2 or U 3 ) are treated in the same way. For instance "D-P" stands for a dual treatment for degrees of freedom associated to direction U 1 , and a primal treatment for degrees of freedom associated to direction U 2 .

Two dimensional plane stress problem

We consider a simple second order two-dimensional problem, the structure is an homogeneous square decomposed in square substructures meshed with linear square finite element (Q1 Lagrange). The behavior is linear elastic (Young modulus E = 200000 MPa and Poisson coefficient ν = 0.3), the loading consists in clamping on the left side and punctual effort on the top right corner (figure 17).

Table 1 shows the number of iterations of available strategies for different number of elements per subdomain (for a 16-subdomain decomposition), and table 2 for different number of subdomains (for given ratio H/h = 16). Globally all approaches (primal, dual and hybrid) equipped with their best preconditioner and projector behave similarly and are scalable. Note that even in its optimal configuration the hybrid approach requires a Finally for such simple problems, the simplified versions of the dual approach gives excellent results.

Bending plate

We consider a forth order plate problem, the structure is an homogeneous square decomposed in squared substructures meshed with square Mindlin plate element. The behavior is linear elastic (Young modulus E = 200000 MPa and Poisson coefficient ν = 0.3), the loading consists in clamping on the left side and punctual normal effort on the top right corner.

Table 3 presents the number of iterations for the dual and primal approaches, with or without optional corner constraints (the subscript indicates the total size of coarse problems, ie rigid body motions and corner modes). As predicted, corner constraints are essential in order to make the algorithms scalable. Anyhow the dimension of the coarse space associated to corners quickly explodes which makes the methods less interesting from a CPU time point of view, which justifies the FETIDP philosophy which leads to much smaller coarse spaces. scaling is employed. What appears clearly is the good behavior of the approaches face to heterogeneity except the identity projector of the dual approach (which is definitely not suited to heterogeneous problems), and the efficiency of the condensed initialization. For such a problem the superlumped projector leads to very good results, anyhow for more complex cases Dirichlet projector is necessary and shall be improved at no extra computational cost by the condensed initialization. Table 5 enables to highlight the fundamental role played by the decomposition: scalabily result only holds for structured decomposition; moreover there might by a huge performance gap between two decompositions with the same number of subdomains (factor 3 for 27 subdomains, factor 4 for 64 subdomains).

Bitraction test specimen

In order to assess "real life" problems, we consider the bitraction test specimen presented in figure 2 (this structure, courtesy of ONERA -Pascale Kanouté -, was optimized with ZeBuLoN software in order to have stress field as homogeneous as possible in its center As shown in table 6 all methods give excellent performance results on non-academical problem. Note the good behavior of the hybrid approach which gives equivalent performance with much smaller coarse problem, even if no physical consideration could guide the choice of the treatment of interface degrees of freedom.

Conclusions

In this paper, we have reviewed most used non-overlapping domain decomposition methods. These methods are perfectly suited to modern computational hardware, they are based on very close concepts which we tried to outline. We introduced the hybrid framework to include as many methods as possible: the principle is to assign to each interface degree of freedom its own treatment, for now primal and dual treatments have been implemented, and mixed and recondensed should follow. Hybriding methods also enable to define physicfriendly approaches for multifield problems.

Because of the conceptual proximity of all methods, assessments showed very close numerical performance results. Once equipped with convenient preconditioner and coarse problem, all the methods proved their ability to handle second and forth order elasticity in presence of strong heterogeneities. Though from a computational point of view some combination may be more interesting: dual approach with lumped preconditioner or simplified projector (if these are sufficient to ensure fine rate of convergence), hybrid approach (which generates smaller coarse space). We also outlined the importance of the decomposition even for very simple problems. The methods have also proved their efficiency on industrial cases, some of them were implemented in computational softwares.

In this paper we limited to the solution to linearized systems, which anyhow enables to solve nonlinear problems. Another strategy is to commute the nonlinear solver and the domain decomposition method so that nonlinear problems can be solved independently on each subdomain. Another evolution of domain decomposition philosophy is the decomposition of the time interval [START_REF] Jl | Résolution d'edp par un schéma en temps pararéel[END_REF] for evolution problems.

First author wish to thank ECCOMAS and ARCME for their kind proposal to write this paper.

A Krylov iterative solvers

Krylov iterative solvers for the resolution of linear systems have been widely studied. The aim of this section is just to briefly present important results and algorithms, reader interested in wider documentation can refer to [START_REF] Saad | Iterative methods for sparse linear systems[END_REF], and to [START_REF] Barrett | Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods[END_REF] for shorter explanation.

Krylov methods belong to the projection class of iterative algorithms, which consist in approximating solution S -1 b of system Sx = b by vector p(S)b where p is a smartly built polynomial.

In this section we consider the iterative solution to system Sx = b. S is a n × n matrix and b a vector in range(S). The i th iteration leads to approximation x i of the solution, associated residual is r i = b -Sx i = S(xx i ). Initialization is x 0 (most often x 0 = 0). Canonical (orthonormal) basis of R n reads (e 1 , . . . , e n ).

A.1 Principle of Krylov solvers

Krylov solvers are based on the iterative construction of so-called "Krylov subspace" K m (S, r 0 ) defined by: K m (S, r 0 ) = span r 0 , . . . , S m-1 r 0

The solution to linear system consists in searching x m under the following constraints:

x m ∈ x 0 + K m (S, r 0 ) r m ⊥ ? K m (S, r 0 ) (2)
where the choice of the orthogonality relationship enables to define various approaches.

A.2 Most used solvers

We herein present two of the principal Krylov solvers. First GMRes [START_REF] Saad | GMRes: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] which is suited to any type of matrix, then conjugate gradient which is adapted to symmetric definite positive matrices.

Of course, the iterative solution to a linear system assumes that a convergence criterion is employed, and that a limit of precision is set so that the system is supposed to have converged once the criterion is below this precision. We note ε this limit value of the criterion.

A.3 GMRes

Algorithm GMRes (alg. A.1) consists in an oblique projection based on the construction of Krylov subspace K m (S, v 0 ) with v 0 = r 0 / r 0 2 . The research principle is:

x m ∈ x 0 + K m (S, r 0 ) r m ⊥ SK m (S, r 0 ) (3)
which is equivalent to finding x m ∈ x 0 + K m (S, r 0 ) minimizing r m 2 .

A particulary striking property of GMRes is not to compute the approximation at each iteration, a smart implementation of GMRes enables to directly access the norm of the residual r j 2 . Only the final approximation x m is computed (by the inversion of a m × m upper triangular matrix). From the computation complexity point of view, each iteration consists in a full orthonormalization of vector w j with respect to K j .

Algorithm GMRes(m) (or restarted GMRes) consists in stopping computation before convergence at a priori fixed step m and restarting computation using previous x m as initialization. This strategy aims at minimizing orthogonalization computations by limiting the size of Krylov subspaces [START_REF] Erhel | Restarted GMRes preconditioned by deflation[END_REF]. This method may lead to stagnation for non positive definite matrices. Compute w j = Sv j 4:

for i = 0, . . . , j do 5:

h ij = (v i , w j ) 6:
w j = w jh ij v i 7:

end for 8:

h (j+1)j = w j 2

9:

if r j 2 ε then v j+1 = w j /h (j+1)j 13:

end if 14: end for 15: Compute y m minimizing r 0 2 e 1 -Hm y 2 and set x m = x 0 + V m y m

A.4 Conjugate gradient

Let S be a symmetric positive definite matrix, a conjugate gradient algorithm consists in an orthogonal projection. The research principle is:

x m ∈ x 0 + K m (S, r 0 ) r m ⊥ K m (S, r 0 ) (4)
which is equivalent to finding x m ∈ x 0 + K m (S, r 0 ) minimizing x mx S . Because of the properties of S, conjugation (orthogonality) properties appear, leading to algorithm A. α j = (r j , r j )/(Sw j , w j ) 4:

x j+1 = x j + α j w j 5:

r j+1 = r jα j Sw j 6:

β j = (r j+1 , r j+1 )/(r j , r j ) 7:

w j+1 = r j+1 + β j w j 8: end for (r m ) (residual basis) is orthogonal, (w m ) (research direction basis) is S-orthogonal. Step 6-7 of algorithm A.2 is the S-orthogonalization of w j+1 with respect to w j which theoretically implies the orthogonality of w j+1 with respect to all previous research directions. However this orthogonality property is numerically lost as the number of iterations increases, it is then more suited to use a full orthogonalization of research directions leading to algorithm A.3.

Full reorthogonalization is often compulsory for complex simulations. Various implementations are available (among others Gram-Schmidt, modified Gram-Schmidt, iterative Gram-Schmidt [START_REF] Lingen | Efficient Gram-Scmidt orthonormalisation on parallel computers[END_REF][START_REF] Frayssé | On the influence of the orthogonalization scheme on the parallel performance of GMRes[END_REF]) depending on the chosen ratio between precision and computational cost. Our experience leads us to prefer modified Gram-Schmidt algorithm (the one used in algorithm A.1) to classical Gram-Schmidt (the one used in algorithm A.3). Note that once α j = (r j , r j )/(Sw j , w j ) 4:

x j+1 = x j + α j w j 5:

r j+1 = r jα j Sw j 6:

For 0 i j, β i j = -(r j+1 , Sw i )/(w i , Sw i )

7:
w j+1 = r j+1 + j i=1 β i j w i 8: end for fully reorthogonalized, conjugate gradient is almost as expensive as GMRes. Anyhow conjugate gradient provides the approximation at each iteration, which can be very useful (see for instance section 3.2.3, where such an information enables the computation of relevant convergence criterion).

A.5 Study of the convergence, preconditioning

Because of their error-minimization property, conjugate gradient and GMRes have convergence theorems with known minimal convergence rate, for instance for conjugate gradient:

x -x m S 2 √ κ -1 √ κ + 1 m x -x 0 S ( 5
)
where κ is the condition number of matrix S. Condition number is the ratio between the biggest and the smallest eigenvalues.

Moreover performance results of Krylov iterative solvers are strongly linked to the spectrum of matrix S. More precisely only the active spectrum (set of eigenvalues which the right hand side is not orthogonal to the associated eigenvectors) influences the convergence; condition number κ can be replaced with active condition number κ act inside relation [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF] leading to better convergence range. More precise study would lead to the introduction of Ritz spectrum and effective condition number [107].

These simple considerations are sufficient to explain the interest of preconditioning the system: the idea is to solve equivalent system S -1 Sx = S -1 b where S -1 is a well-chosen matrix providing the system with better spectral properties (if S -1 ≈ S -1 then condition number is optimal, which justifies the notation).

For conjugate gradient, the use of preconditioner may seem problematic since the symmetry is a priori lost. However if preconditioner S -1 is symmetric definite positive, applying conjugate gradient to nonsymmetric system is equivalent to a symmetric resolution ((L -T SL -1 )(Lx) = L -T b with Cholesky factorization S = L T L) and the method is still relevant.

So preconditioning the above two algorithms is simply realized replacing S by S -1 S and b by S -1 b in lines 1 and 3 of algorithm A.1, and in lines 1 and 5 of algorithm A.3 (anyhow the research directions are still S-orthogonal). Of course the main problem remains the definition of an efficient preconditioner.

A.6 Constrained Krylov methods, projector implementation

We may deal (for instance in the dual approach) with constrained systems such as:

S G G T 0 x α = b e (7)
Because constraint G T x 0 = e is compulsory, it is often referred to as "admissibility constraint". A classical solution is to find an initialization x 0 which satisfies constraint and then ensure that the remainder of the solution is researched inside a supplemental space: G T (x ix 0 ) = 0. A projected algorithm naturally arises:

x = x 0 + P x * (8)

G T x 0 = e (9)
G T P = 0 [START_REF] Brenner | Lower bounds in domain decomposition[END_REF] which leads to:

x 0 = QG G T QG -1 e (11)
P = I -QG G T QG -1 G T ( 12
)
where Q is a matrix so that matrix G T QG is invertible. Iterative system then reads:

P T SP x * = P T (b -Sx 0 ) (13)
α can be post-computed α = G T QG -1 G T Q (b -Sx).

A.7 Augmented-Krylov methods, projector implementation

Augmented-Krylov methods [START_REF] Chapman | Deflated and augmented Krylov subspace techniques[END_REF][START_REF] Saad | Analysis of augmented Krylov subspace methods[END_REF] are employed to add optional constraints to the resolution of a system. The principle is to set subspace C in R n of dimension n c represented by n × n c rectangular matrix C (range(C) = C, for more simplicity we suppose that C is full-ranked-column), then to define augmented-Krylov subspace K m (S, r 0 , C) = K m (S, r 0 ) + range(C), and to use the following research principle:

x m ∈ x 0 + K m (S, r 0 , C) r m ⊥ ? K m (S, r 0 , C)

Augmented-Krylov methods can be implemented either by reorthogonalization schemes or by projection methods which are the one we propose to present here. The research space is separated into two subspaces: range(C) and a supplemental subspace. The part of the solution in range(C) is detected during initialization, while the remainder is iteratively looked for, projector P ensures the research is realized inside the correct subspace.

x = x 0 + P x * (15)

C T r 0 = C T (b -Sx 0 ) = 0 (16)
C T SP = 0 [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] Which leads to: Though it can be proved that projected system is better conditioned than original problem [START_REF]Conjugate gradient method with preconditioning by projector[END_REF], the efficiency of the method essentially depends on the choice of matrix C, which is most often an opened problem. Within the framework of domain decomposition methods, this choice can be guided by several considerations. In the framework of multiresolution, reuse of previous numerical information can lead to very interesting performance results [START_REF] Rey | Reuse of krylov spaces in the solution of large-scale non linear elasticity problems[END_REF][START_REF] Gosselet | On a selective reuse of krylov subspaces in newton-krylov approaches for nonlinear elasticity[END_REF][START_REF] Gosselet | A domain decomposition method for quasi incompressible formulations with discontinuous pressure field[END_REF][START_REF] Rey | Solution to large nonlinear systems: acceleration strategies based on domain decomposition and reuse of krylov subspaces[END_REF][START_REF] Rey | A Rayleigh-Ritz preconditioner for the iterative solution to large scale nonlinear problems[END_REF][START_REF] Risler | Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale nonlinear problems[END_REF] A.8 Constrained augmented Krylov methods

x 0 = C C T SC -1 C T b (18)
P = I -C C T SC -1 C T S (19)
We here consider solving constrained system [START_REF] Bhardwaj | Application of the FETI method to ASCI problems: Scalability results on a thousand-processor and discussion of highly heterogeneous problems[END_REF] with C-augmented algorithm. The admissibility constraint is often referred to as first level constraint and augmentation as second level constraint. Two strategies are possible, the first consists in mixing levels together while the second respects the hierarchy between constraints.

One projector strategy: set J = G S T H and ẽ = -e H T b , then system reads:

S J J T 0 x α = b ẽ (22)
the following initialization/projection are employed (Q is a parameter to tune)

x 0 = J J T QJ -1 ẽ (23)

P = I -QJ J T QJ -1 J T ( 24
)
Because Q is not easy to interpret and choose, this method is hardly ever used.

Two-projector strategy: the two conditions are imbricated. First ensure admissibility constraint

x = x 0 + P x * (25)
x 0 = G G T QG -1 e (26)

P = I -QG G T QG -1 G T (27)
then set

x * = x * 0 + P * x * * (28)

x * 0 = C C T P T SP C -1 C T P T (b -Sx 0 ) (29)
P * = I -P C C T P T SP C -1 C T P T S (30) so that optimality constraint is verified. As can be seen such an approach is equivalent to classical augmentation with making second level constraints consistent with the first (setting C * = P C).

Figure 1 .

 1 Figure 1. Schwarz' original problem

Figure 2 .

 2 Figure 2. 16 subdomains bitraction test specimen (courtesy of ONERA -Pascale Kanouté)

Figure 4 .

 4 Figure 4. Two-subdomain decomposition

Figure 8 .

 8 Figure 8. Representation of first iteration of preconditioned primal approach

Figure 9 .

 9 Figure 9. Representation of first iteration of preconditioned dual approach

p

  + T (s) is invertible, we can express displacement u (s) b from local equilibrium equation, and suppress it from global interface conditions, which leads to:   I I -T (1) + T (2) S

Figure 10 .

 10 Figure 10. Representation of first iteration of preconditioned hybrid approach

Figure 11 .

 11 Figure 11. Representation of the first iteration of preconditioned FETIDP

Figure 12 .

 12 Figure 12. Representation of the first iteration of preconditioned BDDC

ΥFigure 13 .

 13 Figure 13. Different kinds of heterogeneity in domain decomposition context

Figure 14 .

 14 Figure 14. Topological interface information

Figure 15 .

 15 Figure 15. Algebraic interface objects

Figure 16 .

 16 Figure 16. Articulation between formulation and solver

Figure 17 .

 17 Figure 17. 16 subdomains decomposed square

Figure 18 . 3D assessments 6 . 4

 1864 Figure 18. 3D assessments

Algorithm A. 1 1 :

 11 GMRes Compute r 0 = b -Sx 0 , v 0 = r 0 / r 0 2 2: for j = 0, . . . , m -1 do 3:

2 . 1 :

 21 The algorithm is based on the construction of various basis of K m (S, r 0 ) : Algorithm A.2 Conjugate gradient Compute r 0 = b -Sx 0 , set w 0 = r 0 2: for j = 0, . . . , m do 3:

Algorithm A. 3 1 :

 31 Reorthogonalized conjugate gradient Compute r 0 = b -Sx 0 , set w 0 = r 0 2: for j = 0, . . . , m do 3:

κ = λ n λ 1 with |λ 1 |

 1 |λ 2 | . . . |λ n | eigenvalues of S

  system then reads:SP x * = b -Sx 0 (20)orP T SP x * = P T (b -Sx 0 ) = P T b(21)

Table 1 .

 1 Scalability results in 2D / 16 subdomains smaller coarse space (for instance, table 1, two 12 × 12 coarse problems against one 36 × 36 coarse problem for primal or dual approaches) for equivalent efficiency. As expected if the optimality coarse problem is suppressed, performance results decay and scalability is lost.

Table 5 .

 5 Homogeneous cube / influence of the decomposition

	❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ ❤ Method Decomposition	Structured	Unstructured
		27 sd. 64 sd. 27 sd. 64 sd.
	Primal Neumann-Neumann	11 108	14 288	42	67
	Dual Dirichlet P(I)	12 108	16 288	44	69
	Dual Dirichlet P(D)	12 108	16 288	43	70
	Hybrid D-D-P P(I)	14 72	17 192	-	-
	Hybrid P-P-D P(I)	15 72	19 192	-	-
	Hybrid D-P-D P(I)	13 72	17 192	-	-

Table 6 .

 6 ). It was decomposed with Metis software into 4 or 16 subdomains. Bitraction test specimen

	Method		4 sd.	16 sd.
	Primal Neumann-Neumann	23 0+12 30 0+69
		Lumped P(I)	32 12	41 69
	Dual	Dirichlet P(I)	25 12	32 69
		Dirichlet P(Q)	24 12	32 69
	Hybrid P(Q)	P-P-D D-D-P	31 8 25 8	44 46 37 46

ACKNOWLEDGEMENTS

We wish to thank Frederic Feyel and Daniel Rixen for coauthoring some of our works. More specifically (but not exclusively) Frederic for his contribution when implementing our methods inside ZeBuLoN software. Moreover, working with them is always a great pleasure.