T Hirschler
email: thibaut.hirschler@insa-lyon.fr

R Bouclier

D Dureisseix

A Duval

T Elguedj

J Morlier

A Dual Domain Decomposition Algorithm for the Analysis of Non-conforming Isogeometric Kirchhoff-Love Shells

Keywords: Isogeometric Analysis, Shell, Mortar coupling, Multipatch, Domain Decomposition, Parallel computing

Originally, Isogeometric Analysis aimed at using geometric models for the structural analysis. The actual realization of this objective to complex real-world structures requires a special treatment of the non-conformities between the patches generated during the geometric modeling. Different advanced numerical tools now enable to analyze elaborated multipatch models, especially regarding the imposition of the interface coupling conditions. However, in order to push forward the isogeometric concept, a closer look at the algorithm of resolution for multipatch geometries seems crucial. Hence, we present a dual Domain Decomposition algorithm for accurately analyzing non-conforming multipatch Kirchhoff-Love shells. The starting point is the use of a Mortar method for imposing the coupling conditions between the shells. The additional degrees of freedom coming from the Lagrange multiplier field enable to formulate an interface problem, known as the one-level FETI problem. The interface problem is solved using an iterative solver where, at each iteration, only local quantities defined at the patch level (i.e. per sub-domain) are involved which makes the overall algorithm naturally parallelizable. We study the preconditioning step in order to get an algorithm which is numerically scalable. Several examples ranging from simple benchmark cases to semi-industrial problems highlight the great potential of the method.

Introduction

A recent challenging issue in computational mechanics aims at reunifying geometric modeling and structural analysis. The idea emerged with the introduction of IsoGeometric Analysis by Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. The initial idea was to use NURBS-based models coming from Computer Aided Geometric Design to directly perform the analysis. This has been shown possible by using the NURBS basis functions to build the approximation spaces for the Finite Element Analysis. It has shown a great performance, not only because the geometry is exactly preserved, but also because NURBS have interesting mathematical properties as for instance higher inter-element continuity than classical Lagrange polynomials. Thus, IGA has been successfully applied to solve demanding problems as shell analysis [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF][START_REF] Bouclier | An isogeometric locking-free NURBS-based solid-shell element for geometrically nonlinear analysis[END_REF], shape optimization problems [START_REF] Wall | Isogeometric structural shape optimization[END_REF][START_REF] Nagy | Isogeometric design of anisotropic shells: Optimal form and material distribution[END_REF][START_REF] Hirschler | Isogeometric sizing and shape optimization of thin structures with a solid-shell approach[END_REF], material with non-linear behaviors [START_REF] Elguedj | B and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements[END_REF][START_REF] Bouclier | Development of a mixed displacement-stress formulation for the analysis of elastoplastic structures under small strains: Application to a locking-free, NURBS-based solid-shell element[END_REF][START_REF] Ambati | Isogeometric Kirchhoff-Love shell formulation for elasto-plasticity[END_REF], contact problems [START_REF] Seitz | Isogeometric dual mortar methods for computational contact mechanics[END_REF][START_REF] Occelli | LR B-Splines implementation in the Altair RadiossTM solver for explicit dynamics IsoGeometric Analysis[END_REF], fluidstructure interaction [START_REF] Kamensky | Immersogeometric cardiovascular fluidstructure interaction analysis with divergence-conforming B-splines[END_REF][START_REF] Apostolatos | Systematic evaluation of the interface description for fluidstructure interaction simulations using the isogeometric mortar-based mapping[END_REF], etc. However, IGA faces some difficulties related to the generation of analysis-suitable models. In fact, it is known that taking a raw geometric model, as generated by a CAD modeler, is usually not possible due to the presence of trimmed regions and non-conformities between the patches. There are mostly two ways to deal with this issue. Firstly, one can strive to remove all trimmed regions and non-conformity by invoking re-parametrization strategies which can be cumbersome for volume models [START_REF] Xu | Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications[END_REF][START_REF] Al Akhras | Towards an automatic isogeometric analysis suitable trivariate models generationApplication to geometric parametric analysis[END_REF][START_REF] Massarwi | Volumetric untrimming: Precise decomposition of trimmed trivariates into tensor products[END_REF]. The second option consists in using advanced numerical tools for analyzing models with the mentioned defects coming from the geometric modeling [START_REF] Marussig | A Review of Trimming in Isogeometric Analysis: Challenges, Data Exchange and Simulation Aspects[END_REF][START_REF] Teschemacher | Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis[END_REF]. For that, robust procedures and algorithms have to be employed and this work is heading in this direction: it intends to formulate an efficient algorithm for analyzing multipatch shell models with non-conforming interfaces.

Regarding the issue of non-conforming interfaces, the first point to tackle concerns the imposition of the coupling conditions. The coupling can be addressed by different techniques, and in the specific case of arbitrary non-conforming interfaces, one may resort to weak coupling approaches. There exist three principal classes of methods, namely penalty coupling [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF][START_REF] Breitenberger | Analysis in computer aided design: Nonlinear isogeometric B-Rep analysis of shell structures[END_REF][START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF], Mortar coupling [START_REF] Brivadis | Isogeometric mortar methods[END_REF][START_REF] Bouclier | Development of a new, more regular, mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches[END_REF][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: Geometric details, local refinement, inclusion, fracture[END_REF][START_REF] Dornisch | Dual and approximate dual basis functions for B-splines and NURBS Comparison and application for an efficient coupling of patches with the isogeometric mortar method[END_REF][START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF][START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF][START_REF] Wunderlich | Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity[END_REF][START_REF] Schuß | Multi-patch isogeometric analysis for Kirchhoff-Love shell elements[END_REF], and Nitsche coupling [START_REF] Nguyen | Nitsche's method for two and three dimensional nurbs patch coupling[END_REF][START_REF] Schillinger | The non-symmetric Nitsche method for the parameterfree imposition of weak boundary and coupling conditions in immersed finite elements[END_REF][START_REF] Guo | A parameter-free variational coupling approach for trimmed isogeometric thin shells[END_REF][START_REF] Nguyen-Thanh | Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling[END_REF][START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF]. We focus in this work on shell analysis and especially on the case of isogeometric Kirchhoff-Love shell formulations. In this particular context, all three methods have been successfully applied for large and complex industrial structures, even for non-linear analysis [START_REF] Herrema | Penalty coupling of non-matching isogeometric KirchhoffLove shell patches with application to composite wind turbine blades[END_REF][START_REF] Schuß | Multi-patch isogeometric analysis for Kirchhoff-Love shell elements[END_REF][START_REF] Guo | Variationally consistent isogeometric analysis of trimmed thin shells at finite deformations, based on the STEP exchange format[END_REF]. The differences between these approaches mainly result from the formulation of the interface constraints and how they are imposed during the resolution. In this sense, penalty methods and Mortar methods share similarities. These two methods usually impose kinematic constraints between the displacement fields of neighboring patches. The kinematic constraints vary depending on the mechanical behavior of the patches (3D solids, shells, etc.) but also the geometric nature of the junctions (C 0 , G 1 , or C 1 junctions for example). The transfer of the stress resultants is implicitly taken into account which simplifies the mathematical expression of these two coupling methods in contrast to the Nitsche approach where the stress transfer has to be explicitly formulated. Thus, among all methods, various arguments are put forward regarding criteria as for example the simplicity, the accuracy, the robustness, or even the range of applications. Nevertheless, a great potential of the Mortar method seems not to have been yet fully tapped: there exists a close link between this class of coupling method and the family of Domain Decomposition algorithms.

In fact, the additional Degrees Of Freedom associated to the Lagrange multiplier field of the Mortar approach enables to formulate an interface problem which is naturally suitable for parallel computing [START_REF] Stefanica | A Numerical Study of FETI Algorithms for Mortar Finite Element Methods[END_REF][START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF]. Efficient parallel solvers for large-scale systems have been studied over past decades which have led for example to an entire class of methods, namely the non-overlapping Domain Decomposition Methods [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]. The idea is to subdivide the computational domain into non-overlapping sub-domains: instead of solving one large system, multiple independent local systems are solved per sub-domain and continuity conditions are introduced at the sub-domain interfaces using Lagrange multipliers iteratively. More specifically, the Finite Element Tearing and Interconnecting [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] and the Balancing Domain Decomposition [START_REF] Tallec | Domain decomposition methods for large linearly elliptic three-dimensional problems[END_REF] are often seen as originators of the two main branches of those nonoverlapping DDM. BDD is usually referred to as a primal approach since it chooses the interface displacement field as main unknown, whereas FETI is usually referred to as a dual approach because it privileges the interface loads. Kleiss et al. [START_REF] Kleiss | IETI -Isogeometric Tearing and Interconnecting[END_REF] notice that a close link exists between DDM and IGA since multipatch NURBS discretizations can be seen as a natural Domain Decomposition (each sub-domain being a patch). Thus, the Dual-Primal FETI solver from Farhat et al. [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF] has been employed in the IGA framework under the name Isogeometric Tearing and Interconnecting (IETI) [START_REF] Kleiss | IETI -Isogeometric Tearing and Interconnecting[END_REF][START_REF] Hofer | Dual-primal isogeometric tearing and interconnecting solvers for multipatch dG-IgA equations[END_REF][START_REF] Stavroulakis | Non-overlapping domain decomposition solution schemes for structural mechanics isogeometric analysis[END_REF]. In this work, we are looking for a DDM to analyze non-conforming Kirchhoff-Love shells. In the classical Finite Element Method framework, the FETI-DP (and BDDC) has been especially developed for fourth-order problems [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 | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF][START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF]. However, in case of non-conforming interfaces, this approach seems to be difficult to apply: it is less trivial to identify the equivalent of the so-called corner nodes (subset of global DOF where exact continuity is enforced) since non-conforming parametrizations do not share common DOF by definition. We thus resort to the initial one-level FETI problem to build our dual Domain Decomposition algorithm. We also employ the Mortar method presented in our previous work [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF] where the kinematic constraints allow to couple Kirchhoff-Love shells that intersect in non-conforming way, i.e. where the interface can cross some elements and the angle between both shells can be arbitrarily defined. The motivation of this work is to build a solver that is efficient from the computational point of view (suitable for parallel computing) but also accurate (quality of the resolution, especially at the non-conforming junctions).

The developed dual Domain Decomposition solver is presented as follows. Firstly, in section 2 we remind major aspects concerning multipatch modeling. Specifically, we clarify the notion of non-conforming parametrizations in the context of isogeometric shells. Then, we present the main theoretical points regarding Kirchhoff-Love shells. Two versions are used in this work: the standard isogeometric shell from Kiendl et al. [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF] and the embedded formulation from Hirschler et al. [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF] which simplifies complex geometric modeling tasks. Next, we review the Mortar approach of [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF] that enables to properly couple multiple non-conforming patches. Section 3 deals with the formulation of the interface problem. We detail the different steps necessary to get and solve this system, and we give information about implementation aspects. Finally, section 4 highlights the potential of the method on different examples with an increasing level of complexity. This brings us to section 5 where concluding remarks and perspectives are discussed.

Non-conforming multipatch analysis

Multipatch modeling is inevitable when dealing with structures with non-trivial topology using IGA. In this section, we first remind main aspects about multipatch discretizations. Then, we give the basics regarding the Kirchhoff-Love shell formulation. Especially, we point out the case of an embedded adaptation of the classical Kirchhoff-Love shell [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF]. This new Kirchhoff-Love formulation involves non-isoparametric elements which require to perform particular choices and extensions for the developed FETI solver. Finally, we review the Mortar approach introduced in former contribution [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF] that enables to properly couple non-conforming Kirchhoff-Love shells. It constitutes the starting point for the derivation of the dual DD solver.

Multipatch NURBS geometries

The framework of this study is based on NURBS modeling, widely used in the context of IGA. A NURBS patch provides a parametric description of a curve, a surface, or a volume. A key point regarding the construction of multivariate NURBS is the tensor-product step introduced in the definition of the basis functions. For instance, a NURBS surface is a bivariate vector-valued function of the form

S(u, v) = n i=1 m j=1 R i j (u, v)P i j . (1)
The piecewise rational basis functions R i j of degree p in the u direction and degree q in the v direction are obtained by taking weights w i j associated to the bidirectional net of control points P i j , and the products of univariate B-spline basis functions

R i j (u, v) = N p i (u)N q j (v) w i j n k=1 m l=1 N p k (u)N q l (v) w kl . (2)
Finally, the univariate B-spline basis functions N p i and N q j are defined by two non-uniform knots span of size r = n + p + 1 and s = m + q + 1 respectively. More information can be found, for example, in Piegl and Tiller [START_REF] Piegl | The NURBS Book[END_REF] and Rogers [START_REF] Rogers | An Introduction to NURBS: With Historical Perspective[END_REF].

NURBS are recognized as a powerful tool for geometric design. For more than one decade, they have also shown great potential for numerical simulation with the concept of isogeometric analysis introduced by Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]. However, some inherent properties of NURBS complicate the automatic transfer from a NURBS-based CAD model to an analysis-suitable model. More specifically, a single NURBS patch cannot represent complex industrial geometries. Instead, multiple NURBS patches are required where each of them describes a specific part of the overall structure. Furthermore, if the constraint of conforming or even matching interfaces is required, the geometric model may need to be modified in order to perform the analysis. From here on, we use the following terminology to characterize the coupling [START_REF] Bouclier | Development of a new, more regular, mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches[END_REF]:

• it is non-conforming when some elements are overlapped by the interface,

• it is conforming but non-matching when the interface is aligned with the edges of the patches but the control points on both sides do not match,

• it is matching when the discretizations are exactly the same on both sides of the interface.

Figure 1 describes the process for non-matching and non-conforming interfaces. In the simplest case of a nonmatching interface (see figure 1, left), knots insertion (or removal) and degree elevation (or reduction) allow to recover a final matching mesh. In the more general case of a non-conforming interface (see figure 1, right), a complete reparametrization of some parts is inevitable. The local changes introduced to reach the matching condition at the interface have global influence on the whole discretization. For instance, C 0 lines propagate along the patch as depicted in figure 1. This issue is inherent to the construction of NURBS. The tensor-product step during the definition of the basis function (2) complicates the imposition of local modifications. In order to prevent the use of cumbersome reparametrization procedures, we seek to perform the analysis directly on non-conforming multipatch configurations. Finally, let us notice that trimming techniques are extensively used in CAGD [START_REF] Marussig | A Review of Trimming in Isogeometric Analysis: Challenges, Data Exchange and Simulation Aspects[END_REF][START_REF] Teschemacher | Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis[END_REF]. We do not deal with these procedures in this work. However, the following developments are also applicable in trimmed IGA and the constructed Domain Decomposition algorithm with the Mortar coupling considered in this work could be of great interest.

Based on these observations, we can formulate two statements:

• NURBS modeling generates a natural decomposition of the domain (subparts of a global structure).

• NURBS modeling of complex structures naturally results in non-conforming discretizations.

These points motivate the development of efficient solvers that deal with non-conforming multipatch NURBS-based models.

Shell formulations 2.2.1. Continuum mechanics

For the shell analysis, we employ the isogeometric Kirchhoff-Love shell [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF]. It is based on the Kirchhoff kinematic assumptions which consists in neglecting the effect of transverse shear strains. Consequently, within this theory, it is postulated that the director vector remains orthogonal to the mid-surface in the deformed configuration. These assumptions enable to describe the deformation of the shell body only by the displacement u of the mid-surface S. The overall displacement field U of the shell body takes the form:

U(θ 1 , θ 2 , ζ) = u(θ 1 , θ 2 ) + ζ [Φ × A 3 ] (θ 1 , θ 2 ), (3)
where ζ ∈ [-t 2 ; t 2 ], t being the thickness of the shell. The curvilinear coordinates θ 1 and θ 2 are basically the parameters of the bivariate function (e.g. NURBS) representing the mid-surface. The linearized rotation vector Φ is defined as follows [START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF]:

Φ = ϕ 1 A 1 + ϕ 2 A 2 , (4)
V(ξ 1 ,ξ 2 ,ξ 3 ) S(θ 1 ,θ 2 ) + u(θ 1 ,θ 2 )
Figure 2: Overview of the embedded isogeometric Kirchhoff-Love shell elements as introduced in [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF]. The displacement field is approximated by the same parametrization as the embedded surface. However, the embedded surface does not describe the mid-surface of the shell. An additional mapping step further modifies the geometry of the shell.

where the rotation angles ϕ 1 and ϕ 2 are given by :

ϕ 1 = 1 J u ,θ 2 • A 3 and ϕ 2 = - 1 J u ,θ 1 • A 3 . (5)
The subscript (•) ,θ i = ∂(•)/∂θ i indicates the partial derivative with respect to variable θ i . Finally, the normal unit vector A 3 is obtained using the covariant basis vectors A 1 and A 2 :

A 1 = S ,θ 1 , A 2 = S ,θ 2 , A 3 = 1 J A 1 × A 2 with J = | A 1 × A 2 |. (6)
We willfully omit to give all the intermediary steps that lead to the final expressions of the virtual works for the considered Kirchhoff-Love shells. These theoretical developments can be found, for example, in [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | Isogeometric Analysis and Shape Optimal Design of Shell Structures[END_REF][START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF][START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF][START_REF] Bischoff | Models and Finite Elements for Thin-Walled Structures[END_REF]. The virtual works take the following forms:

δW KL int = - Ω 0 n : δε + m : δκ dΩ 0 , δW KL ext = Ω 0 p • δu dΩ 0 + Γ 0 q • δu dΓ 0 , (7)
where ε and n denote the membrane strains and forces, and κ and m denote the bending strains and moments. The external loads p describe distributed loads per unit of area applied on the mid-surface Ω 0 , and q are axial forces per unit of length applied on the edges of the patch Γ 0 .

Embedded Kirchhoff-Love shell

We consider two possibilities for the definition of the mid-surface. Firstly, the classical approach where a NURBS surface is used which will be referred to as the standard Kirchhoff-Love shell. The element formulation is obtained by introducing in equation ( 7) the same NURBS-based discretization for the mid-surface and the displacement field:

S(θ 1 , θ 2 ) = R i j (θ 1 , θ 2 )P i j , u(θ 1 , θ 2 ) = R i j (θ 1 , θ 2 )u i j . (8)
Secondly, we will also employ the embedded approach as introduced in Hirschler et al. [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF]. The main idea is to describe the mid-surface by the composition of a NURBS volume V and an immersed NURBS surface S:

S(θ 1 , θ 2 ) = V S(θ 1 , θ 2 ) . (9)
Mathematically, the surface S and the volume V are defined using a bi-variate NURBS and a tri-variate NURBS, respectively. We have:

S(θ 1 , θ 2 ) = Rab (θ 1 , θ 2 ) Pab , and V(ξ 1 , ξ 2 , ξ 3 ) = R i jk (ξ 1 , ξ 2 , ξ 3 )P i jk . (10)
The approximation space of the displacement field is build using the discretization of the embedded surface:

u(θ 1 , θ 2 ) = Rab (θ 1 , θ 2 )u ab . (11)
Hence, this embedded Kirchhoff-Love shell has the particularity of not using the isoparametric principle in contrast to the standard isogeometric Kirchhoff-Love shell. We will show that it has further influence, especially on the description of the null space and the construction of the pseudo-inverse brought into play later in section 3. Figure 2 summarizes the embedded approach, and further information can found in Hirschler et al. [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF] if needed. Finally, the only difference with the standard Kirchhoff-Love shell is that an additional step (referred as the Free-Form Deformation step) transforms the mid-surface. This embedded approach enables us to build a model of a wing with its internal sub-structure (see section 4.3). Lastly, let us note that the use of two types of elements highlight the versatility of the method tackled in this paper. The dual Domain Decomposition solver is independent of the element formulation, and thus can be applied in general cases. At the end, the equilibrium configurations of the shell followed from the principle of minimum potential energy lead to the typical linear system:

K u = f , (12)
where the vector u concatenates the degrees of freedom (DOF) of the displacement field. Again, the interested readers can refer to [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Kiendl | Isogeometric Analysis and Shape Optimal Design of Shell Structures[END_REF][START_REF] Cirak | Subdivision surfaces: a new paradigm for thin-shell finite-element analysis[END_REF][START_REF] Kiendl | Isogeometric Kirchhoff-Love shell formulations for general hyperelastic materials[END_REF][START_REF] Bischoff | Models and Finite Elements for Thin-Walled Structures[END_REF] for the expressions of the stiffness matrix K and the load vector f for Kirchhoff-Love-type models.

Mortar coupling 2.3.1. Formulation of the coupled problem

Let us now consider the case of two shells occupying, respectively, domains Ω (1) and Ω (2) , and joining at interface Γ. We assume the case of non-overlapping subdomains, i.e. the junction is described by a curve. In case of fastened shells, the continuity between both displacements has to be ensured at the interface. Furthermore, in case of a shell junction with rigid hinge, an additional constraint is required. The continuity between the rotations in the tangential direction associated with the interface curve also needs to be enforced. These two kinematic constraints are formulated as follows:

u (1) = u (2) on Γ, [3 displacements] Φ (1) • t = Φ (2) • t on Γ, [1 rotation] ( 13
)
where t is a unit tangent vector associated to the interface curve. The constraints (13) cannot be imposed in a strong manner for arbitrary non-conforming parametrizations. A static condensation approach presented in Coox et al. [START_REF] Coox | A flexible approach for coupling NURBS patches in rotationless isogeometric analysis of Kirchhoff-Love shells[END_REF] can be applied in case of patches joining at their edges only. However, if a junction crosses the middle of a patch (i.e. for a non-conforming interface) then this approach may not appear appropriate. To overcome this issue, the continuity constraints are ensured in a weak sense. More precisely, a Mortar method is used in this work because of its direct link with FETI-like solvers [START_REF] Stefanica | A Numerical Study of FETI Algorithms for Mortar Finite Element Methods[END_REF][START_REF] Lacour | Analyse et résolution numérique de méthodes de sous-domaines non conformes pour des problèmes de plaques[END_REF]. Two Lagrange multipliers λ ∈ L d and λ ∈ L r are introduced (L d and L r being ad-hoc spaces for the displacements and the rotations respectively). The coupled problem is obtained by formulating the following Lagrangian:

L u (1) , u (2) , λ, λ = 1 2 a 1 u (1) , u (1) -l 1 u (1) + 1 2 a 2 u (2) , u (2) -l 2 u (2) + λ, u (1) -u (2) + λ, Φ (1) • t -Φ (2) • t , (14)
where bilinear forms a s and linear forms l s constitute the standard variational forms of the elasticity problem on each subdomain. Operator •, • is a bilinear form defined such that:

v, w = Γ v • w dl . ( 15
)

Approximation spaces

The variational principle written in the discrete form gives the coupled linear system to be solved:

                         K (1) 0 C (1)T Z (1)T 0 K (2) C (2)T Z (2)T C (1) C (2) 0 0 Z (1) Z (2) 0 0                                                   u (1) u (2) Λ Λ                          =                          f (1) f (2) 0 0                          , (16)
where vectors u (1) , u (2) , Λ, and Λ collect the DOF corresponding to the discretizations of u (1) , u (2) , λ, and λ respectively:

u (s) h (θ 1 , θ 2 ) = i R (s) i (θ 1 , θ 2 ) u (s) i s = 1, 2, λh (ξ) = i Ri (ξ) λi , λh (ξ) = i Ri (ξ) λi , (17)
The displacement coupling matrices C (s) and the rotation coupling matrices Z (s) consist in sparse signed rectangular operators.

The linear system ( 16) is a saddle point problem, implying that the approximation spaces of the Lagrange multipliers should be carefully chosen in order to satisfy the inf-sup condition. Recent works present optimal approximation spaces for Mortar coupling in isogeometric analysis [START_REF] Brivadis | Isogeometric mortar methods[END_REF][START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF][START_REF] Wunderlich | Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity[END_REF][START_REF] Schuß | Multi-patch isogeometric analysis for Kirchhoff-Love shell elements[END_REF]. In particular, Schuß et al. [START_REF] Schuß | Multi-patch isogeometric analysis for Kirchhoff-Love shell elements[END_REF] deal with the G 1 coupling of non-matching Kirchhoff-Love shells. Based on these works, one could build optimal approximation spaces for the Mortar coupling of shells that intersect in an arbitrary, non-conforming manner (i.e. where the junction crosses a patch, or where the patches intersect with an angle, etc.). This issue is not tackled in this work. Instead, the strategy adopted here is only based on numerical experiments. Let p denote the smallest degree of the discrete displacement field of both subdomains Ω (1) and Ω (2) . We define:

• for the displacement constraints, λh as a vector-valued B-Spline function with degree p -1 since it is mainly related to traction forces,

• for the rotation constraint, λh as a scalar-valued B-Spline function with degree p -2 because it transfers a bending moment,

• an identical mesh refinement for both Lagrange multipliers. We discretize these fields using as many elements as the coarsest of the trace of domains Ω (1) and Ω (2) along the interface.

We assert that with such a choice, we never encountered instabilities in our computations. The numerical examples presented in section 4 highlight the good behavior of this Mortar coupling. One can also notice that the choice of the approximation spaces for the Lagrange multipliers is completely decoupled from the overall approach applied in this work to build the Domain Decomposition solver. With other choices in hand, the remaining part of the method performs identically.

Implementation aspects

Regarding the construction of the coupling operators C (s) and Z (s) , they are obtained by assembling local matrices, as for the stiffness matrix K (s) . Since the construction is done in the same way for each subdomain, we skip the identification through superscript (s). Let l denote the global DOF number of the j th -component associated to the displacement vector u b , such that u l = u b • e j (where e j refers to the basis of the global Cartesian coordinate system). Similarly, let k denote the global DOF number of the Lagrange multipliers. With these notations in hand, one can express the following variations:

∂u ∂u l = R b e j , ∂ λ ∂ λk = Ra e i , ∂ λ ∂ λk = Ra . ( 18
)
x y z

Ω (1) Ω (2) θ 2 θ 1 θ 2 θ 1 θ 2 θ 1 C (1) (ξ) P 1 (1)
P 2

(1)

P 3 (1) θ 2 θ 1 P 1 (2) C (2) (ξ) P 2 (2)
Physical Space : Parameter Spaces

C(ξ)

Non-conforming Interface Figure 3: Non-conforming parametrization: the interface between both surfaces is represented in their parameter spaces by using NURBS curves C(s) . These curves are generated in order to ease the building of the coupling operators for the Mortar approach.

This way, we can express the coupling matrices as:

C kl = ± Γ ∂ λ ∂ λk • ∂u ∂u l dl = ± Γ Ra R b δ i j dl (19)
Z kl = ± Γ ∂ λ ∂ λk ∂Φ ∂u l • t dl = ± Γ Ra A 3 J R b,θ 2 ( A 1 • t) -R b,θ 1 (A 2 • t) • e j dl (20)
The symbol ± indicates whether the subdomain is the master (+) or the slave (-).

Looking at expressions [START_REF] Teschemacher | Realization of CAD-integrated shell simulation based on isogeometric B-Rep analysis[END_REF] and [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF], we still need to clarify how to integrate numerically over the interface. For general non-conforming interface, the methodology adopted in this paper is depicted in figure 3. We manage to represent the interface in the parameter spaces through NURBS curves. As mentioned by Bauer et al. [START_REF] Bauer | Embedded structural entities in NURBS-based isogeometric analysis[END_REF], it results from a surface-to-surface intersection problem which is a common task in CAD. Under the consideration that no approximation is involved during the generation of the curves, we have two definitions of the interface curve (one per subdomain). In other words, the following equation applies:

S (1) C(1) (ξ) = S (2) C(2) (ξ) = C(ξ) ∀ξ ∈ [0, 1]. (21)
The embedded curve C allows to perform a change of variables which results in the following treatment to compute the coupling matrices. For instance, for matrix C, we get:

C kl = ± [0,1] Ra (ξ) R b C(ξ) δ i j |T(ξ)| dξ, (22)
which is also true for matrix Z. The vector T is given by:

T(ξ) = d dξ S C(ξ) = C,ξ (ξ) • e 1 A 1 + C,ξ (ξ) • e 2 A 2 . ( 23
)
Vector T is tangential to the interface curve. Therefore, once normalized, one can obtain the unit tangent vector t which comes into play in the rotation coupling matrix [START_REF] Apostolatos | Domain decomposition methods and kirchhoff-love shell multipatch coupling in isogeometric analysis[END_REF]. Finally, numerical integration is performed using a Gaussian quadrature rule. A refinement by knot insertion of the embedded NURBS curves discretizes the interface in a finite number of segments. The number of segments and the number of Gauss points per-segment should be chosen such that the numerical error remains negligible. At the end, the assembly of the coupling matrices is performed with a loop on the Gauss points of each segment:

1 0 f (ξ) dξ = I ξ I+1 ξ I f (ξ) dξ ≈ I gp f (ξ I gp )J I gp ω gp , (24)
where ξ I gp denotes the position of the Gauss point on the parameter space of the embedded curve, J I gp denotes the Jacobian of the last mapping such that the integral over [ξ I , ξ I+1 ] is changed into an integral over [-1, 1], and ω gp denotes the quadrature weight.

Generalization

Let us now generalize the coupling approach for the case of n s sub-domains (i.e. n s patches) Ω (i) , i = 1 . . . n s . We define as K = diag(K (1) , K (2) , . . . , K (n s ) ) the diagonal block matrix that stores all the stiffness matrices of each sub-domain. We concatenate the displacement and the rotation coupling matrices into single coupling matrices such that B (i) = [C (i)T , Z (i)T ] T . Thus, we introduce the global coupling operator B = [B (1) , B (2) , . . . , B (n s ) ]. We also define the global load vector f = ( f (1)T , f (2)T , . . . , f (n s )T ) T . Finally, the total displacement and Lagrange multiplier DOFs are stored in vectors u and λ respectively. With these notations, the generalization of the coupled system (16) reads:

          K B T B 0                     u λ           =           f 0           . ( 25
)

The proposed dual Domain Decomposition solver

The proposed solver belongs to the family of dual Domain Decomposition algorithms, originally introduced by Farhat and Roux [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF] with the one-level FETI algorithm. Since then, a large number of papers and variants have been developed to successfully apply in many contexts: e.g. plate and shells [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 | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF][START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF], heterogeneous problems [START_REF] Gosselet | Simultaneous FETI and block FETI: Robust domain decomposition with multiple search directions[END_REF][START_REF] Bovet | Adaptive multipreconditioned FETI: Scalability results and robustness assessment[END_REF], contact problems [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF], etc. The following developments are based on these researches and provide an efficient algorithm for the parallel computation of non-conforming multipatch Kirchhoff-Love shells.

Formulation of the interface problem

The first step of the approach consists in splitting problem (25) into the following coupled set of equations:

K (i) u (i) = f (i) -B (i)T λ, for i = 1, . . . , n s (26)
i

B (i) u (i) = 0. ( 27
)
The goal is then to formulate a dual interface problem with the Lagrange multiplier DOF as the only unknowns. This interface problem is obtained by introducing the local equilibriums (26) into the coupling condition [START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF]. Generally, the local stiffness matrices K (i) are not invertible. Indeed, if no Dirichlet boundary conditions are imposed on the floating sub-domain Ω i then the matrix K (i) is singular. Therefore, K (i) being symmetric, equation ( 26) has a solution if and only if the load vector f (i) -B (i)T λ belongs to the image of operator K (i) . It leads to an additional equation, called the admissibility condition, which is:

R (i)T f (i) -B (i)T λ = 0, for i = 1, . . . , n s (28)
where R (i) is a rectangular matrix. Its columns describe a basis of the null space ker(K (i) ) of matrix K (i) . From a mechanical point of view, these vectors are the local rigid body modes. Then, for each sub-domain, the equilibrium is given by:

u (i) = K †(i) f (i) -B (i)T λ + R (i) α (i) , (29)
where α (i) is a vector that collects the different amplitudes of the rigid body modes (or null-energy modes) for subdomain i and K †(i) is a pseudo-inverse of the stiffness matrix, i.e. that satisfies

K (i) K †(i) K (i) = K (i) .
Now, the local solutions ( 29) can be introduced into the compatibility condition [START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF]. It yields the following equation i

B (i) K †(i) B (i)T λ - i B (i) R (i) α (i) = i B (i) K †(i) f (i) . (30)
Finally, combining this last equation with the admissibility condition given by ( 28), we end up with the interface problem. The resulting system is generally called the FETI system and is given by:

          F G G T 0                     λ α           =           d e           , ( 31
)
where α concatenates the rigid body unknowns of each sub-domain. The dual Schur complement F, the constraint matrix G, and the right-hand sides d and e are expressed as:

F = i B (i) K †(i) B (i)T , G = -B (1) R (1) B (2) R (2) • • • B (n s ) R (n s ) , d = i B (i) K †(i) f (i) , e = -f (1)T R (1) f (2)T R (2) • • • f (n s )T R (n s ) T .

Resolution

The resolution of system [START_REF] Nguyen | Nitsche's method for two and three dimensional nurbs patch coupling[END_REF] relies on an iterative solver where only matrix-vector products are performed. This way, there is no need for assembling the dual Schur complement which would be computationally very demanding since pseudo-inverses K †(i) are involved. Moreover, it avoids to invert the dual Schur complement which has a dense structure. Even more important, the product between a vector y and the dual Schur complement F results in a sum of local contributions which is ideally suited to parallel computing:

Fy = i B (i) v (i) , where v (i) = K †(i) B (i)T y. ( 32
)
More specifically, F being symmetric positive definite, a Preconditioned Conjugate Projected Gradient (PCPG) algorithm is used. To satisfy (28), a projector P onto the null space of G T , i.e. onto ker(G T ), is introduced. Thus, by definition, we have the relations:

G T P = 0, and P T G = 0.

It allows to eliminate the unknowns α from the first equation of system [START_REF] Nguyen | Nitsche's method for two and three dimensional nurbs patch coupling[END_REF]. To this end, we multiply this equation to the left by P T :

P T Fλ + Gα = d ⇒ P T Fλ = P T d. (33)
The second equation G T λ = e of the system (31) is explicitly fulfilled by introducing the splitting:

λ = λ 0 + P λ, ( 34
)
where the initialization λ 0 satisfies the constraint G T λ 0 = e. This splitting is operated during the resolution by firstly choosing λ 0 as the starting guess, and secondly by projecting the gradient during the update. Finally, these two projection steps transform the interface problem (31) into the projected interface problem:

P T FP λ = P T d -Fλ 0 , (35)
Initialize

λ 0 = QG(G T QG) -1 e r 0 = d -Fλ 0 Iterate k = 1, 2, . . . until convergence w k-1 = P T r k-1 Project z k-1 = Mw k-1 Precondition y k-1 = Pz k-1
Re-project which is solved using the PCPG algorithm (see table 1). We still need to express the projector and the initialization. It brings into play a symmetric matrix denoted Q for which the product G T QG is invertible. Matrix Q can be taken as being the preconditioner, identity or a scaling matrix [START_REF] Rixen | Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver[END_REF]. Finally, we have:

p k = y k-1 - k-1 i=1 β i p i ( p 1 = y 0 ) Conjugate with β i = y k-1 • Fp i p i • Fp i α k = y k-1 • w k-1 p k • F p k Minimize λ k = λ k-1 + α k p k Update r k = r k-1 -α k Fp k Residual
λ 0 = QG G T QG -1 e, (36)
P = I -QG G T QG -1 G T . ( 37
)
For most problems, the simplest choice Q = I is the most computationally efficient, and many studies reported in the literature have been performed with this choice [START_REF] Stefanica | Parallel FETI algorithms for mortars[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]. We apply the same treatment herein. Finally, full or partial re-orthogonalization of the gradient at every iteration is commonly adopted to overcome the poor conditioning of the dual Schur complement [START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF].

Pseudo-inverse and Null space

As of now, we did not mention how to formulate the null space R (i) and the pseudo-inverse K †(i) that come into play in the one-level FETI solver.

Explicit construction of the kernel

The null space of a stiffness matrix is commonly associated with the space of rigid-body modes. For 3D elasticity problem and under the condition that no Dirichlet boundary conditions are imposed, we know that there exist 6 rigidbody motions consisting of 3 translations (modes numbered 1, 2, and 3) and 3 rotations (modes numbered 4, 5, and 6) around 3 non-collinear axes. Considering a shell with mid-surface S(θ 1 , θ 2 ), an arbitrary kernel element r is :

r(θ 1 , θ 2 ) = a + b × S(θ 1 , θ 2 ), ( 38
)
where a ∈ R 3 defines a translation and b ∈ R 3 defines a rotation. In order to build a basis of the null space, one can successively describe the translations and the rotations with respect to the global unit vectors e i . We can re-write equation [START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF] in the discrete form in order to identify the DOF corresponding to these unitary rigid-body modes:

k R k (θ 1 , θ 2 )r i k = e i , and k R k (θ 1 , θ 2 )r i+3 k = e i × S(θ 1 , θ 2 ), (39)
where R k are the NURBS basis functions of the discretized displacement field. It is obvious that taking all DOF r i k equal to e i enables to describe the translation rigid-body modes. For the rotations, it depends on the discretization of the shell mid-surface:

• in the case of the standard Kirchhoff-Love shell, the displacement field is approximated using the same discretization as for the mid-surface which gives:

k R k (θ 1 , θ 2 )r i+3 k = k R k (θ 1 , θ 2 ) e i × P k ⇒ r i+3 k = e i × P k . (40)
• in the case of the embedded Kirchhoff-Love shell, the displacement field and the mid-surface are discretized differently (R k for the displacement and M l for the mid-surface), leading to the linear system of equations:

k R k (θ 1 , θ 2 )r i+3 k = l M l (θ 1 , θ 2 ) e i × P l , ∀ (θ 1 , θ 2 ). (41)
System ( 41) cannot be solved in general. One can choose a set of parameters {(θ 1 , θ 2 ) 1 , (θ 1 , θ 2 ) 2 , . . . } and build an approximation of the null space by interpolation for instance. However, it does not ensure that equation ( 41) is exactly satisfied (i.e. true for every couple of parameters (θ 1 , θ 2 )) which means that this approximation of the null space may not be energy-free. Thus, we cannot explicitly build the kernel of the stiffness matrix for the embedded shell formulation. It means that the approach presented in Dostál et al. [START_REF] Dostál | Total FETI-an easier implementable variant of the FETI method for numerical solution of elliptic PDE[END_REF], called the Total-FETI, which partly aims at simplifying the implementation of the original FETI method, is not applicable in this context. Without knowing the null space modes a priori, we cannot opt for a geometrical strategy, based for example on fixing nodes and on a regularization of the stiffness matrix as in Kozubek et al. [START_REF] Kozubek | Total FETI domain decomposition method and its massively parallel implementation[END_REF]. We resort to a purely algebraic method presented in the following section to build the null spaces and the pseudo-inverses involved in the algorithm.

Algebraic method

For clarity, we voluntarily skip here the sub-domain identification through superscript (i), but the reader should view the following operators as local quantities. Because we were not able to explicitly build the null space R in the case of the embedded formulation, we resorted to the purely algebraic method presented in Farhat and Géradin [START_REF] Farhat | On the general solution by a direct method of a large-scale singular system of linear equations: application to the analysis of floating structures[END_REF] and already used in Farhat and Roux [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. The pseudo-inverse and the null space are constructed during the factorization of the stiffness matrix. It is based on the partition of matrix K ∈ R n×n as:

K = K 11 K 12 K 12 T K 22 , (42)
where K 11 ∈ R r×r is a submatrix with full rank r equal to the rank of the initial matrix (i.e. r = rk K = rk K 11 ). With such a partition in hand, the following relation holds:

K 22 -K 12 T K -1 11 K 12 = 0. ( 43
)
One can define the pseudo-inverse and the null space as:

K † = K -1 11 0 0 0 , R = -K -1 11 K 12 I , ( 44
)
where I is the identity matrix of size m = nr.

In practice, the partition ( 42) is not known in advance but it is obtained in an unordered form during a Gaussianbased factorization (as for example a LU decomposition). At each step of the factorization, one can check if a null pivot is encountered. If true, the corresponding equation is redundant and is removed from the system. The reduced column can be recycled to later recover the null space. Furthermore, a small tolerance is set to detect the null pivots (i.e. the right term in [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method-part I: A faster alternative to the two-level FETI method[END_REF] is, in practice, not exactly equal to zero). It enables for instance to track the rotation rigidbody modes in the case of the embedded Kirchhoff-Love shell. More information can be found in Farhat and Géradin [START_REF] Farhat | On the general solution by a direct method of a large-scale singular system of linear equations: application to the analysis of floating structures[END_REF]. Alternative methods, based for example on the Singular Value Decomposition (SVD) are also available [START_REF] Golub | Calculating the Singular Values and Pseudo-Inverse of a Matrix[END_REF].

Preconditioning

Last ingredient involved in the FETI algorithm is the preconditioner already introduced in table 1 through operator M. A good preconditioner is crucial for the computational performance of the iterative resolution. Several preconditioners have been studied in the literature for solving different projected interface problems. In the context of FEM matching grids, the Dirichlet preconditioner is said to be optimal regarding the asymptotic bound and thus offers excellent numerical scalability [START_REF] Rixen | Extended preconditioners for the FETI method applied to constrained problems[END_REF]. In case of Mortar interface, we can re-write this Dirichlet preconditioner as:

F-1 I = n s i=1 B (i) B (i)T -1 B (i) 0 0 0 S (i) bb B (i)T B (i) B (i)T -1 , (45)
where matrix S (i) bb denotes the primal Schur complement and is given by:

S (i) bb = K (i) bb -K (i) bi K (i) ii -1 K (i) ib , (46)
where indices b correspond to the boundary DOF on the interface and indices i correspond to the internal DOF. The Dirichlet preconditioner ( 45) is obtained by exactly inverting the local dual Schur complements and summing each local contribution. More specifically, some works deal with the particular case of Mortar coupling and improvements to the initial version of the Dirichlet preconditioner have been studied in this context [START_REF] Stefanica | A Numerical Study of FETI Algorithms for Mortar Finite Element Methods[END_REF][START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF][START_REF] Lacour | Analyse et résolution numérique de méthodes de sous-domaines non conformes pour des problèmes de plaques[END_REF][START_REF] Rixen | Extended preconditioners for the FETI method applied to constrained problems[END_REF]. For instance, scaling factors can be added in order to take into account the relative stiffnesses across the interface [START_REF] Rixen | Theoretical comparison of the FETI and algebraically partitioned FETI methods, and performance comparisons with a direct sparse solver[END_REF]. These heterogeneities are particularly numerous in the case of stiffened structures with the presence of T -junctions, different shell thicknesses and material behaviors [START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF]. Thus, it is primordial to efficiently deal with those heterogeneities to analyze complex shell structures. Therefore, in the present paper, it is proposed to follow the extended Dirichlet preconditioner suggested by Rixen [START_REF] Rixen | Extended preconditioners for the FETI method applied to constrained problems[END_REF] and mentioned by Stefanica [START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF] which can be written as:

F-1 D = n s i=1 BDB T -1 B (i) D (i) 0 0 0 S (i) bb D (i) B (i)T BDB T -1 , (47)
where matrix D (i) is chosen here as the superlumped scaling:

D (i) = diag K (i) -1 .
The global matrix D is the diagonal block assembly of all these local scaling matrices, which leads to a diagonal matrix in the present case.

The extended Dirichlet preconditioner has shown better performance in comparison with the original Dirichlet preconditioner in the case of Mortar interfaces [START_REF] Stefanica | A Numerical Study of FETI Algorithms for Mortar Finite Element Methods[END_REF][START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF][START_REF] Lacour | Analyse et résolution numérique de méthodes de sous-domaines non conformes pour des problèmes de plaques[END_REF]. Rixen [START_REF] Rixen | Extended preconditioners for the FETI method applied to constrained problems[END_REF] shows that this preconditioner is mechanically consistent and mentions its mechanical meaning: once applied to the residual of the interface problem [START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF], it computes the correction for the Lagrange multipliers. Stefanica [START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF] discusses the parallelization properties of this generalized preconditioner. The critical point concerns the two global multiplications by (BDB T ) -1 which cannot be done in parallel. One may build local matrices B(i) = (BDB T ) -1 B (i) once and for all at the beginning of the resolution. In order to preserve the sparsity, biorthogonal Mortars could be helpful in this case [START_REF] Wunderlich | Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity[END_REF][START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF].

Numerical investigation

In this section, we run several examples to highlight the good behavior of the algorithm. Firstly, 2D test cases are presented, then we study several shell test problems, and lastly we analyze two cases of stiffened structures. For every examples, we depict the number of iterations versus the level of refinement of the mesh given a fixed decomposition, with or without the preconditioner [START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF]. The convergence criterion is defined as a relative tolerance on the projected residual (namely w k in table 1) set as ||w k || < 10 -8 ||w 0 ||. The results are given in tables where the headers N tot and N λ denote the numbers of total and interface DOF respectively (i.e., by referring to [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: Geometric details, local refinement, inclusion, fracture[END_REF], N λ is the size of λ, and N tot is the total number of unknowns). Columns with headers I and F-1 D give the iteration number without and with the generalized preconditioner, respectively. The goal of these tables is to assess the numerical scalability of the DD algorithm with the mesh refinement. To show the accuracy of the coupling, we study the convergence of the solution for two examples: the cantilever beam and the Scordelis-Lo roof. In order to do so, we compute the relative energy error as: where the reference energy E mono is computed on a very fine single patch discretization. 1A two non-matching patches with a similar geometry, 1B four non-matching patches with a similar geometry, 1C four non-matching patches with different geometries. For every configuration, the patches are discretized with cubic NURBS.

E h coupling = | E mono -E coupling | E mono , ( 48
)
1C ( ref 0 ) 1B ( ref 0 ) Displacement Stress (Von-Mises) 1A ( ref 0 )
1A 1B 1C ref N tot N λ I F-1 D E h coupling N tot N λ I F-1 D E h coupling N tot N λ I F-1 D E h
As can be seen in figure 4, the displacement and stress fields are smooth and clearly express the bending of the beam. Then, table 2 summarizes the convergence results for the FETI algorithm for this first test case. One can see that the algorithm requires very few iterations to reach convergence when the preconditioner (47) is used. For each decomposition, it converges in a relatively constant number of iterations. It accounts for the numerical scalability of this preconditioner which is in good agreement with the observations of Stefanica [START_REF] Stefanica | A Numerical Study of FETI Algorithms for Mortar Finite Element Methods[END_REF][START_REF] Stefanica | Parallel FETI algorithms for mortars[END_REF] and Rixen [START_REF] Rixen | Extended preconditioners for the FETI method applied to constrained problems[END_REF]. For finer refinement levels, the gain is significant in comparison with the un-preconditioned version. It is especially true for the more complex decomposition 1C.

2A 2B ref N tot N λ I F-1 D N tot N λ I F-
Additionally, the results of the convergence study are given in table 2. The relative energy error is decreasing with the mesh refinement. It indicates that the solutions obtained with the presented solver for non-conforming meshes converge to the single patch solution (reference solution). We observe a same convergence rate for each configuration. Furthermore, it can be noticed that this convergence rate is identical to the one of the single patch discretization.

Heterogeneous domain. The second test case consists in a heterogeneous beam. The description of the problem is given in figure 5. The beam is made from a sequence of hard and soft regions. One side of the beam is fixed and the other side is subjected to a traction load. Each hard region is defined by a 5 m by 4 m rectangle and has a Young's Modulus equal to E = 60 × 10 9 Pa. Each soft region is defined by a 3 m by 4 m rectangle and has a Young's Modulus equal to E = 6 × 10 9 Pa. Same Poisson's ratio equal to ν = 0.30 is chosen for both regions. We consider two cases:

2A three non-conforming patches; two made of hard material separated by one made of soft material, 2B five non-conforming patches; three made of hard material separated by two made of soft material.

Table 3 gives the convergence results for the present case. As for the previous example of the homogeneous cantilever beam, the algorithm requires very few iterations to reach convergence when the preconditioner is used, and this number of iterations does not grow with the refinement level. This does not hold true without the preconditioner. Therefore, the heterogeneity is correctly tackled by the preconditioned. For instance, for configuration 2B with refinement level 3, the algorithm converges in 38 iterations without a preconditioner, while it requires 10 iterations only when the generalized Dirichlet preconditioner is equipped with the superlumped scaling (see table 3). Let us also emphasize here that the good behavior of the preconditioner is mostly due to the superlumped scaling. For the same configuration and without the scaling (i.e. taking matrix D as the identity in expression (47)), the algorithm would have required 30 iterations to converge which is almost as much as without any preconditioner.

3A 3B ref N tot N λ I F-1 D N tot N λ I F-

Shell test cases

Plate bending. The third example consists in the bending of a plate as described in figure 6. The plate is fully fixed at one side while a uniform pressure is applied over the whole structure. We decompose the plate into two patches and we consider two cases:

3A two patches with the same thickness (homogeneous plate), 3B two patches with different thicknesses such that the ratio is t (1) /t (2) = 0.1 (heterogeneous plate), Cubic NURBS are used and the discretization of the patch with the fixed edge is twice finer than the second one (see figure 6). For this particular test case, we observe better results (quality of the stress field) when the Lagrange Multipliers are discretized using the trace space of the finer patch. We opt for this choice here even if the strategy presented in section 2.3.2 leads to acceptable results too. Finally, the parameters of the study are as follows : length L = 20 m, width l = 10 m, Young's Modulus E = 60 × 10 9 Pa, Poisson's ratio ν = 0.20, and uniform pressure P = 1.0 Pa. In the homogeneous case, the thickness of both patches is t = 0.1 m. In the heterogeneous case, the thickness of the fixed patch is t (1) = 0.1 m and the thickness of the free patch is t (2) = 1.0 m. Table 4 presents the performance study of the algorithm for this problem. The results are similar for both the homogeneous and the heterogeneous configurations. More specifically, without the preconditioner the algorithm converges in almost the same number of iterations in both cases. It can be seen that, with preconditioner the algorithm converges in slightly less iterations for the heterogeneous plate. Furthermore, the number of iterations is kept low even for the finer meshes when the preconditioner is used. Ideally, this number of iterations should be constant versus the mesh refinement as observed for the plane stress problems 4.1. Here, we do not have a perfect numerical scalability but it is not surprising. In fact, the one-level FETI equipped with the Dirichlet preconditioner is known to be numerically scalable for second-order problems. However, it does not hold true for fourth-order problems as shown for example by Farhat et al. [START_REF] Farhat | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF]. That is why extensions of the classical FETI method, as for instance the two-level FETI and latter on the FETI-DP, have been introduced to optimally deal with plates and shells [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 | The two-level FETI method Part II: Extension to shell problems, parallel implementation and performance results[END_REF]. Based on these remarks, we understand the increase of the iteration counts versus the mesh refinement and the same behavior is expected for the other shell problems tackled in this work. Nevertheless, for this test case of a plate, the presented dual domain decomposition algorithm performs well and leads to entirely acceptable results.

4A 4B 4C ref N tot N λ I F-1 D E h coupling N tot N λ I F-1 D E h coupling N tot N λ I F-
Scordelis-Lo roof. The Scordelis-Lo roof is part of the shell obstacle course which is widely used to study the performance of shell formulations [START_REF] Belytschko | Stress projection for membrane and shear locking in shell finite elements[END_REF]. It consists of a portion of cylinder subjected to a vertical gravity load and fixed at its two end sections using rigid diaphragms. The study parameters and some results for this test case can be found in the abundant shell literature: [START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF][START_REF] Belytschko | Stress projection for membrane and shear locking in shell finite elements[END_REF][START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and B-Method[END_REF]. We decompose the Scordelis-Lo roof into three non-matching configurations as depicted in figure 7:

4A three patches with a curved interface obtained by using the embedded Kirchhoff-Love formulation, 4B four quasi-identical patches which intersect at a cross point, 4C nine patches with different levels of refinement. The coarsest are quartic patches while the finest are cubic patches.

Figure 7 shows the vertical displacement field and the distribution of a bending moment for each configuration, which appears smooth as expected. The results investigating the numerical scalability are given in table 5. It can be seen that the benefit of the preconditioning step is increasing from the first configuration 4A to the third configuration 4C. For the first configuration 4A, the iteration count is only reduced by a factor of 2 or 3 (except for the coarsest mesh). For the second configuration 4B, the iteration count is reduced by a factor of 4 for the refinement levels above 1. Finally, for the more complex decomposition 4C, the iteration count is drastically reduced when using the generalized Dirichlet preconditioner. Looking at the refinement levels 3 and 4 in table 5 (and columns referred to the configuration 4C), we observe that the number of iterations is reduced by a factor of 10, and 15 respectively. Thus, the benefit provided by the preconditioner is all the more remarkable in case of complex Domain Decomposition. Even if numerical scalability is not fully reached, the growth of the iteration count with the refinement level is very slow for each configuration and it leads to satisfactory results. Additionally, the results of the convergence study are given in table 5. It enables to show that the strategy presented in section 2.3.2 is valid. We do not perform the convergence study for configuration 4C since it contains patches of different degrees (cubic and quartic). For each configuration 4A and 4B, the relative energy error is low and decreases with the mesh refinement. Particularly, the convergence rates for configuration 4A and 4B are similar. Interestingly, this convergence rate is similar to the one we observe for the monopatch discretization (with cubic elements).

Stiffened structures

Stiffened panel. The problem of the stiffened panel is presented in figure 8. It consists of the assembly of a square plate with sub-parts called the stiffeners. We design three stiffeners: two with a parabolic shape and the third one with straight shape. The straight patch intersects the two curvilinear stiffeners. As a result, five non-conforming interfaces are defined: each interface cuts several isogeometric elements. Every patch is discretized using cubic Kirchhoff-Love shells. The material and geometric parameters of the study are given in figure 8. The edges of the panel are fixed (no displacement). The uniform pressure is equal to P = 1000 Pa. We denote by 5A this problem of the stiffened panel.

Once again, we run the analysis for several levels of refinement with and without the preconditioning step. The results of the study are given in table 6. We observe the same behavior than for the previous test cases. Indeed, the iteration count is drastically reduced with the use of the preconditioner and it is all the more true when increasing the refinement level. The numerical scalability is not completely achieved but the growth of the number of iterations remains slow while introducing the preconditioning step. Especially, satisfactory results are obtained even for very fine meshes. For instance, with the refinement level of 4, the dual decomposition domain algorithm converges in 125 iterations even though the total number of DOF is over 100 000. This is a crucial point. Because the algorithm solves an interface problem where the unknowns are the Lagrange Multiplier DOF (see equation ( 35)), we are able to measure, during the resolution, the residual of the coupling condition [START_REF] Zou | Isogeometric Bézier dual mortaring: Refineable higher-order spline dual bases and weakly continuous geometry[END_REF]. The convergence criterion is formulated on this residual which, in other words, ensure the good fulfillment of the coupling condition. Thus, we believe that the presented algorithm is viable and therefore attractive to tackle more complex models as encountered in industrial application.

Aircraft wing. We finally present an example of a wing in order to highlight the potential of our algorithm to handle more sophisticated models. We consider the case of a subsonic aircraft wing adapted from Vassberg et al. [START_REF] Vassberg | Development of a Common Research Model for Applied CFD Validation Studies[END_REF] with the following dimensions: the span is 15 m, the sweep angle is 35 • , the root and tip chords are equal to 7 m and 1 m respectively. We build the wing with its internal substructure made of ribs and spars. To this purpose, we adopt the strategy depicted in figure 9. The starting point is to describe the outer geometry which defines the skin. In our case, we use four NURBS surfaces, two for the top and two others for the bottom. Once it has been done, we define two NURBS volumes that fill the domain delimited by the skin. For each volume, we now have the possibility to embed NURBS surfaces in its parameter space. The NURBS composition of the embedded surfaces with the volume mapping leads to surfaces that are perfectly lying inside the domain delimited by the skin and no geometric approximation at the interfaces is introduced. More information can be found in our previous work [START_REF] Hirschler | The embedded isogeometric Kirchhoff-Love shell: From design to shape optimization of non-conforming stiffened multipatch structures[END_REF]. Thanks to the embedded Kirchhoff-Love shell formulation, the geometric model with NURBS compositions can directly be used for the analysis. Finally, it leads to a model with 16 patches, taken as sub-domains for our solver, as described in figure 10 (4 quadratic standard Kirchhoff-Love shells for the skin, 8 and 4 cubic embedded Kirchhoff-Love shells for the ribs and the spars, respectively). In addition to being tedious to generate, an analysis-suitable model of the wing (i.e. with only matching and conforming interfaces as discussed in figure 1) would lead to a much larger number of patches. It could not be achieved without calling on specific geometric modelers. Here, the modeling task is thus simplified and offers a great flexibility.

The results of the analysis are given in figures 10 and 11, and the scalability study of the algorithm is given in table 6 under the mention 6A. The root side of the wing is fixed and a uniform pressure is applied on the upper skins. Under this loading, the magnitude of the displacement at the tip of the wing is about 130 mm. In order to highlight the correct imposition of the coupling conditions, we depict detail views of the different non-conforming interfaces in figure 11. One can see the continuity of the displacement over these interfaces. Concerning the convergence of the algorithm, we were only able to run the resolution without the preconditioning step for the coarsest discretization. In this case, it already takes almost 2000 iterations to converge. For the highest level of refinement, it does not converge before the stopping criteria of maximal iteration number set to 2000. However, the generalized Dirichlet preconditioner drastically accelerates the convergence as shown in table 6. It enables to analyze refined versions of the wing model in an acceptable number of iterations. As for the other examples tackled in this work, we get a slow growth of the iteration count versus the refinement level. Finally, this example demonstrates the crucial need for a good preconditioner to make the algorithm powerful. The generalized Dirichlet preconditioner [START_REF] Amini | A domain decomposition method for problems with structural heterogeneities on the interface: Application to a passenger ship[END_REF] performs well for all examples and has proved to be robust even for more complex cases as the aircraft wing.

Concluding remarks

A dual Domain Decomposition algorithm is presented for the analysis of non-conforming multipatch isogeometric shell models. The starting point is the definition of a Mortar coupling which weakly imposes the coupling conditions through Lagrange Multipliers. In the literature, one argument against Mortar methods is commonly claimed; it concerns the additional unknowns introduced by the Lagrange multiplier field that increase the overall size of the system to be solved. We believe that this argument is not relevant and the present work fully takes advantage of those additional DOF. In fact, it allows us to formulate an interface problem, namely the one-level FETI problem, where the unknowns of the system are these interface DOF. We highlight for different examples with increasing levels of complexity that it enables to accurately analyze non-conforming multipatch structures, even for very fine level of refinement, since we track the interface residual during the resolution. Moreover, the algorithm is naturally parallelizable. Instead of solving one large linear system, only local problems defined at the patch level have to be solved. In fact, the presented algorithm belongs to the FETI family which has been intensively developed over the past decades. Those Domain Decomposition algorithms have shown great performances to analyze very large problems with millions of DOF. Alternatives to the initial one-level FETI, regarding for example the formulation of the interface problem or the preconditioning step, have been suggested and studied in the literature. All those works can be applied to develop efficient analysis solvers for complex industrial isogeometric models. The present work aims to motivate researches in that direction. The direct use of geometric models with non-conforming interfaces and trimming procedures for the structural analysis could be pushed forward with those new Domain Decomposition solvers.

Figure 1 :

 1 Figure 1: NURBS modeling generates multipatch geometries where each patch represents a specific part of the complete structure. Naturally, these patches present non-matching or even non-conforming parametrization at their interfaces. Additional procedures are required to get a final model with conforming and matching interfaces.

Figure 4 :

 4 Figure 4: Description and results of the cantilever beam problem. Three decompositions 1A, 1B, and 1C are studied. The level of refinement (ref i) is given accordingly to the table 2

4. 1 .

 1 Simple plane stress problems Homogeneous domain. The first test case consists in the bending of a 2D cantilever beam with homogeneous material as described in figure4. Note that for this case, a standard Lagrange multiplier formulation is set up to couple the 2D displacements only. The geometric and mechanical parameters are chosen as: length L = 20 m, height h = 5 m, Young's Modulus E = 60 × 10 9 Pa, Poisson's ratio ν = 0.30, and tangential uniform load P = 1000 Pa. Three different configurations are studied and we denote these decompositions as:

Figure 5 :

 5 Figure 5: Description and results of the heterogeneous beam problem. The beam is made of a sequence of hard and soft regions. Configuration 2A has two hard regions separated by one soft region. Configuration 2B has three hard regions separated by two soft regions.

E

  = 60 GPa, ν = 0.20 L = 20 m, l = 20 m t 1 = 0.1 m, t 2 = {0.1 m for 3A, 1.0 m for 3B} Parameters :

Figure 6 :

 6 Figure 6: Description and results of the bending plate problem. The plate is decomposed into two non-conforming patches. Two cases are studied: 3A the homogeneous case where the patches have the same properties, and 3B the heterogeneous case where the shells have different thicknesses.

Figure 7 :

 7 Figure 7: Description and results for the Scordelis-Lo roof test case. Three different non-matching decomposition are studied: 4A is made of three embedded Kirchhoff-Love shells with curved interfaces, 4B is made of four quasi-identical patches (two patches have one more element per direction), and 4C is made of nine patches with different levels of refinement and degrees. (Note: The displacement has unit of length and the bending moment has unit of moment per unit length, respectively).

E(Figure 8 :

 8 Figure 8: Description and results for the stiffened panel problem. It consists in the assembly of a rectangular plate with three stiffeners; two with a parabolic shape crossed by a straight one. It leads to four non-conforming patches and five coupling interfaces (three stiffener/panel and two stiffener/stiffener interfaces). The main stress resultants for the panel are the bending moments whereas the stiffeners are mainly subjected to membrane forces.

Figure 9 :

 9 Figure 9: Construction of the model of the wing: starting from the surface description of the outer geometry, a first step generates volumes that fill the inside of the wing. Then, the internal substructure made of spars and ribs are defined by embedding NURBS surfaces into the parameter space of these volumes. Thus, the surfaces representing the internal substructure are obtained by NURBS compositions.

Figure 10 :

 10 Figure 10: Description and global results of the wing analysis. The model contains a total of 16 patches taken as sub-domains for our solver: 4 for the skin, 4 for the spars, and 8 for the ribs.

  U magn : [ 110 mm; 128 mm ] Rib -Skin U magn : [ 16 mm; 21 mm ] Spar -Spar U x : [ -1.2 mm; -0.35 mm ] Rib -Spar U x : [ -0.21 mm; -0.02 mm ]

Figure 11 :

 11 Figure 11: Detail views of the wing analysis: zooms are depicted in order to highlight the non-conforming interfaces between the different parts (skins, spars, ribs). The color range of the displacement field is scaled accordingly to each zones to show the displacement continuity at these weak coupling interfaces.

Table 1 :

 1 The Preconditioned Conjugate Projected Gradient algorithm with re-orthogonalization.

Table 2 :

 2 Performance study for the homogeneous cantilever beam problem described in figure4. The reference energy was taken as E mono = 5.568207e-2 to compute the relative energy error of the coupled solution.

  [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]

					D
	0 1 2 3 4	204 448 1224 3928 13 944 136 55 16 6 24 9 40 15 72 29	5 7 9 9 10	368 812 2228 7172 144 38 32 10 48 14 80 22 25 508 272 69	6 8 10 10 11

Table 3 :

 3 Performance study for the heterogeneous beam problem described in figure5.

  [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]

					D
	0 1 2 3 4	222 510 1446 4758 140 28 44 76 17 142 268 118 9 14 29 60	6 7 11 18 26	222 510 1446 4758 140 28 44 76 17 142 268 117 7 14 28 59	4 4 6 11 16

Table 4 :

 4 Performance study for the bending plate problem described in figure6.

  [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF]

						D
	0 1 2 3 4	527 1399 4295 14 695 53 927 1038 175 78 47 142 73 270 124 526 158	34 6.15e-3 35 1.19e-4 45 2.85e-6 59 1.01e-7 81 3.71e-9	708 1652 4692 15 380 55 188 1052 405 92 70 156 117 284 188 540 276	34 3.55e-3 38 6.44e-5 48 1.92e-6 64 8.54e-8 86 3.52e-9	1513 3017 7537 22 625 1620 276 468 842 76 993 3156 1783 121 167 44 317 52 558 64 949 87

Table 5 :

 5 Performance study for the Scordelis-Lo shell problem described in figure7. The reference energy was taken as E mono = 4826.577028 to compute the relative energy error of the coupled solution.

Table 6 :

 6

		5A
	ref 0 1 2 3 4	N tot 1227 3313 10 365 35 989 133 317 1827 1106 N λ I 147 109 259 193 483 341 931 613	F-1 D 50 68 85 103 125	ref 0 1 2	6A N λ 9246 2272 1874 N tot I 28 554 4096 -102 444 7744 -	F-1 D 258 299 378

Performance study (left) for the stiffened panel problem 5A described in figure

7

and (right) the wing problem 6A described in figures 9, 10, and 11.