Robin Bouclier
email: bouclier@insa-toulouse.fr

Jean-Charles Passieux
email: passieux@insa-toulouse.fr

A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis

Keywords: Isogeometric analysis, Non-symmetric Nitsche method, Non-intrusive coupling, Non-conforming geometries, Multiresolution

In this work, we propose a new non-intrusive coupling algorithm for global/local isogeometric structural analysis. In contrast to the existing non-intrusive strategies that rely on a Lagrange multiplier coupling, the algorithm makes use of the non-symmetric Nitsche method. It results in an accurate and efficient tool to compute any evolution of a local model within a fixed global NURBS one. The reason for this is the robustness and simplicity of the coupling (no auxiliary fields, no dual space approximation, no stabilization parameters), which enables to directly handle all the non-conforming coupling scenarios encountered through the global/local multiresolution process. The performance of the methodology is numerically demonstrated through a series of two-dimensional elastic benchmarks involving conforming and non-conforming couplings, along straight, curved, and bi-material interfaces. In all examined problems, the proposed Nitsche algorithm provides optimal accuracy. Finally, to illustrate both the efficiency in a multiple query context and the robustness of the method to arbitrary non-conforming scenarios, a simple structural optimization problem is carried out using the developed non-intrusive solver, which simplifies the process and ensures computational time saving.

Introduction

The IsoGeometric Analysis (IGA) concept, which was first introduced in Hughes et al. [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF] and later formalized in the book by Cottrell et al. [START_REF] Cottrell | Isogeometric analysis: Toward Integration of CAD and FEA[END_REF], relies on the use of the same basis functions for the representation of both the geometry in Computer-Aided Design (CAD) and the numerical approximations in Finite Element Analysis (FEA). As a consequence, typical Lagrange polynomials are replaced in the computations by smooth and higher-order functions such as Non-Uniform-Rational-B-Splines (NURBS) [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF][START_REF] Rogers | An introduction to NURBS With Historical Perspective[END_REF], which constitute the most commonly used technology in CAD. Other geometry descriptions include T-splines [START_REF] Bazilevs | Isogeometric analysis using T-splines[END_REF] and subdivision surfaces [START_REF] Cirak | Integrated modeling, finite-element analysis and engineering design for thin-shell structures using subdivision[END_REF]. Within this work NURBS are used. Beyond the reinforced link between CAD and analysis, IGA turned out to be a superior computational mechanics technology, which on a per-degree-of-freedom basis exhibits increased accuracy and robustness in comparison to standard FEA [START_REF] Evans | n-Widths, sup-infs, and optimality ratios for the k-version of the isogeometric finite element method[END_REF]. The reason for this is the higher order of regularity of spline-based functions, namely C (p-1) through the knot-span elements of the mesh for a polynomial degree p, whereas only C 0 continuity is available for Lagrange polynomials. However, in contrast to the standard nodal basis, a multivariate NURBS basis comprises a rigid tensor product structure which precludes the simple modeling of local phenomena. Indeed, we necessarily end up with a structured mesh in a NURBS patch. This unavoidably leads to the overlap of some global knot-span elements to allow for a truly global-mesh independent local region to be incorporated for the modeling of any specific local behaviors (e.g., introduction of a hole [START_REF] Guo | Global-local model coupling for composite shell structures in the framework of isogeometric analysis[END_REF], of an inclusion [START_REF] Nguyen | Nitsche's method for two and three dimensional NURBS patch coupling[END_REF], crack propagation [START_REF] Gupta | Analysis of three-dimensional fracture mechanics problems: A non-intrusive approach using a generalized finite element method[END_REF][START_REF] Passieux | Local/global non-intrusive crack propagation simulation using multigrid XFEM solver[END_REF], emergence of a plastic zone [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF], of local damage [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF], modeling of local contact [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF], etc). Fig. 1(a) illustrates this issue for the example of the integration of a geometrical detail (a hole) within an initial NURBS plate. The coupling of two domains along a non-conforming interface (i.e., an interface that intersects through the global mesh), is involved. The analysis of trimmed NURBS patches needs thus to be addressed, which is known not to be a trivial task in the IGA community [START_REF] Marussig | A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects[END_REF].

From the literature, it seems that three main strategies could be applied to numerically solve the resulting NURBS non-conforming global/local problem. The first strategy consists in recovering the usual situation of boundary fitted discretizations (see Fig. 1(b)). In order to do so, a re-parametrization of the whole global NURBS model is required, leading to the splitting of the new geometry into several patches with C 0 continuity at the boundaries (see, e.g., [START_REF] Wang | Polycube splines[END_REF][START_REF] Akhras | Isogeometric analysis-suitable trivariate NURBS models from standard B-Rep models[END_REF]). This may entail a considerable modeling and computational effort which is often as complex and time consuming as standard mesh generation and then, is opposed to the core idea of IGA.

Concurrently, a second approach initiated in Ruess et al. [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF] and based on the combination of a fictitious domain approach with a weak coupling, may constitute an interesting option to carry out global/local NURBS simulations. To this purpose, the authors made use of the so-called Finite Cell Method (FCM) (see, e.g., [START_REF] Schillinger | The Finite Cell Method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF] for a detailed review). Unlike the first approach, the idea here is to simply use the resulting unfitted structured mesh for the interpolation of the global fields (see Fig. 1(c)), while the trimmed geometry is accurately captured by means of suitable quadrature rules for cut knot-span elements [START_REF] Legrain | A NURBS enhanced extended finite element approach for unfitted CAD analysis[END_REF][START_REF] Nagy | On the numerical integration of trimmed isogeometric elements[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary methods[END_REF][START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF][START_REF] Stavrev | Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method[END_REF]. This strategy appears and has proved to be highly efficient in the context of immersed finite elements. However, let us notice that in the general and more industrially ubiquitous case of a local region that may evolve during the simulation (e.g., to carry out the shape optimization of local entities, or to model crack propagation, or the expansion of a plastic zone, etc), we expect several re-assemblies and re-factorizations of the resulting global/local stiffness operator during the multiresolution process. The left-hand side operators may also appear ill-conditioned depending on the encountered trimmed configuration [START_REF] Schillinger | The Finite Cell Method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF][START_REF] Schillinger | An isogeometric design-throughanalysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces[END_REF].

Relatively connected to the previous approach, the last strategy revolves around the concept of nonintrusiveness. This notion has originally emerged in the standard FE-based community through the so-called non-intrusive global/local coupling methods. Based on an idea put forward by Whitcomb [START_REF] Whitcomb | Iterative global/local finite element analysis[END_REF] and later developed in Gendre et al. [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF], this new class of methods has attracted large attention these last years and has now been successfully applied in many domains in FEA (see, e.g., [START_REF] Gupta | Analysis of three-dimensional fracture mechanics problems: A non-intrusive approach using a generalized finite element method[END_REF][START_REF] Passieux | Local/global non-intrusive crack propagation simulation using multigrid XFEM solver[END_REF] for crack propagation, [START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF] for taking localized uncertainties into account, [START_REF] Guguin | Nonintrusive coupling of 3D and 2D laminated composite models based on finite element 3D recovery[END_REF] for non-intrusive plate/3D coupling, [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF] for non-linear domain decomposition, [START_REF] Chantrait | Low intrusive coupling of implicit and explicit time integration schemes for structural dynamics: Application to low energy impacts on composite structures[END_REF][START_REF] Bettinotti | A fast weakly intrusive multiscale method in explicit dynamics[END_REF] for transient dynamics analysis and [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF] for real aeronautical strcutures). More recently, an extension for NURBS has been proposed in Bouclier et al. [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF] and has proved to be a good candidate for NURBS local enrichment. This strategy seeks to locally modify a global NURBS model without having an impact on its corresponding numerical operators. The method relies on an iterative process between global and local computations, thus involving the initial global non-trimmed NURBS mesh (see Fig. 1(d)). The replacement of part of this mesh locally is carried out exactly and non-intrusively: only interface data are transmitted from one model to the other and the initial global stiffness operator remains unchanged (independently of the shape of the local domain). We emphasize that the terminology "non-intrusive" used here only characterizes the numerical solver but not the geometrical construction (and parametrization) of the local model that still needs to be addressed. We believe that this method is of the greatest interest when addressing NURBS global/local analysis. Especially in case the local region evolves, thus leading to the resolution of a sequence of similar problems (in a multiresolution fashion [START_REF] Néron | Time-space PGD for the rapid solution of 3D nonlinear parametrized problems in the many-query context[END_REF][START_REF] Gosselet | Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems[END_REF]), the method meets its full potential: the global stiffness operator is assembled and factorized only once and the systems to be solved remain well-conditioned, which ensures both robustness and computational time saving.

If in FEA the situation of conforming interfaces appears sufficient for the global/local coupling since the basis functions do not suffer from a global tensor product structure, the challenge when addressing NURBS is to handle non-conforming coupling interfaces (see Fig. 1(a) again). In its standard FEA and first NURBS version, the non-intrusive strategy comes with a Lagrange multiplier based coupling method. Indeed, thanks to the introduction of additional unknowns that represent the interface traction, it appears, at first glance, easier to make the two models communicate through interface exchanges while keeping the initial structure of the global and local stiffness operators. Regarding the more general context of NURBS coupling, such methods involving the definition of dual basis functions have been successfully applied (see, e.g., [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] for standard elasticity, [START_REF] Dornisch | The weak substitution method -An application of the mortar method for patch coupling in NURBS-based isogeometric analysis[END_REF][START_REF] Coox | A robust patch coupling method for NURBS-based isogeometric analysis of non-conforming multipatch surfaces[END_REF] for the weak substitution method and [START_REF] Apostolatos | Domain decomposition methods and Kirchhoff-Love shell multipatch coupling in isogeometric analysis. Isogeometric analysis and applications[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] Sommerwerk | Reissner-Mindlin shell implementation and energy conserving isogeometric multi-patch coupling[END_REF] for the coupling of shells). However, it has to be noted that the performance of these techniques crucially depends on the compatibility between the approximation subspaces considered for the primal and the dual variables. In our case of arbitrary non-conforming coupling configurations, this rules out many choices for the subspaces. As a result, Nitsche-based coupling have been preferred in the field of immersed boundary methods. Nitsche formulations ensure variational consistency by introducing flux terms along the coupling interface expressed by primal unknowns of the coupling domains. Thus, the method is free of auxiliary fields, which simplifies the theory. Nevertheless, in its widespread symmetric version (see, e.g., [START_REF] Annavarapu | A robust Nitsche's formulation for interface problems[END_REF] for the origin and then [START_REF] Nguyen | Nitsche's method for two and three dimensional NURBS patch coupling[END_REF][START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF] for standard elasticity), additional stabilization terms are required to ensure coercivity of the formulation.

Appropriate estimates of the stabilization parameters involving the resolution of an additionnal eigenvalue problem (whose resolution can be performed locally per element [START_REF] Jiang | A robust Nitsche's formulation for interface problems with spline-based finite elements[END_REF]) are thus involved, which may appear delicate from an algorithmic viewpoint. Therefore, there has been an increasing interest in developing methods without such mesh-dependent stabilization parameters [START_REF] Kollmannsberger | Parameter-free, weak imposition of Dirichlet boundary conditions and coupling of trimmed and non-conforming patches[END_REF][START_REF] Baiges | A symmetric method for weakly imposing Dirichlet boundary conditions in embedded finite element meshes[END_REF].

In this context, the non-symmetric version of the Nitsche method has recently reemerged in Schillinger et al. [START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF] for isogeometric immersed methods. Since then, few applications can already be found (for instance, see [START_REF] Guo | A parameter-free variational coupling approach for trimmed isogeometric thin shells[END_REF] for the coupling of shells and [START_REF] Dettmer | A stabilised immersed boundary method on hierarchical b-spline grids[END_REF] which is related to the incompressible Navier-Stokes equations). Originally introduced as a discontinuous Galerkin method [START_REF] Oden | A discontinuous hp finite element method for diffusion problems[END_REF][START_REF] Rivière | A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems[END_REF][START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF], it is based on variationally consistent numerical flux conditions that are introduced in such a way that the criterion for stability is (weakly) satisfied. Therefore, it does not require the introduction of additional stabilization terms in contrast to its symmetric counterpart and thus, its performance does not rely on an appropriate estimation of additional sensitive parameters.

The object of this work is to further improve the NURBS non-intrusive approach [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF] by extending it to incorporate the non-symmetric Nitsche coupling. The goal is to offer simplicity, accuracy and robustness, regardless of the non-conforming coupling scenarios encountered through the global/local multi-resolution process. Indeed, we believe that the combination of the non-symmetric Nitsche coupling with the global/local iterative non-intrusive strategy will result in an accurate and efficient tool to compute any evolution of a local model within a fixed NURBS one. The paper is organized as follows: following this introduction, Section 2 specifies the modeling and the derivation of the non-symmetric Nitsche coupling in the context of linear elasticity. Then, in Section 3, we develop a non-intrusive iterative algorithm that makes use of the non-symmetric Nitsche method for the coupling. It results in an original iterative procedure whose mechanical interpretation is rather different from that of the existing non-intrusive strategies. In section 4, the performance of the proposed methodology over the existing NURBS one is numerically demonstrated through a series of benchmarks involving local mesh refinement, introduction of a geometrical detail, and coupling at a bi-material interface. In addition, the ability of the method to adapt to any evolution of the local model is highlighted by carrying out a simple structural optimization example. Finally, Section 5 concludes on this work by summarizing our most important points and motivating future research based on the proposed algorithm.

Reference mechanical coupling problem and formulation

The framework of this study draws on research dealing with IGA [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric analysis: Toward Integration of CAD and FEA[END_REF] based on NUBRS [START_REF] Piegl | The NURBS Book (Monographs in Visual Communication)[END_REF][START_REF] Rogers | An introduction to NURBS With Historical Perspective[END_REF], nonintrusive global/local analysis [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF], and non-symmetric Nitsche coupling [START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free 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] Dettmer | A stabilised immersed boundary method on hierarchical b-spline grids[END_REF]. More precisely, the interaction of two non-conforming linear elastic NURBS meshes is considered: one domain is the foreground and defines the local model, the other is the background and constitutes the global NURBS model to be enriched. Special care is taken to handle the most general situations we may encounter; that is, for instance, when the shape of the local region is expected to evolve during the simulation. This results in a local mesh that is truly independent from the underlying global NURBS discretization.

Governing equations

Let us start by considering a global NURBS (possibly multi-patch) model. Its mesh is characterized by a physical domain Ω 1 ⊂ R d (d = 2 or 3 being the dimension of the ), which is divided into two disjoint, open and bounded subsets Ω 11 and Ω 12 . Those two non-overlapping subdomains share a common interface denoted by Γ such that Ω 1 = Ω 11 ∪ Ω 12 ∪ Γ and Ω 11 ∩ Ω 12 = ∅ (see Fig. 2(a) (left) for illustration). We assume that this coarse NURBS mesh is sufficient for the modeling except in the small region Ω 12 where a local phenomena is to be introduced. As a consequence, a local (more detailed) "sub-model" characterized by domain Ω 2 is constructed to replace the global model in Ω 12 (see Fig. We assume that the two non-overlapping subdomains Ω and, displacements u g 11 and u g 2 are prescribed over boundaries Γ u11 and Γ u2 (see Fig. 2 again). In each subdomain, the kinematic constraints, the equilibrium equations and the constitutive relations have to be verified. Using subscript m to denote a quantity that is valid over region Ω m , with m = 11 and 2, the corresponding governing equations read:

u m = u g m over Γ um ; (1a) div(σ m ) + f g m = 0 in Ω m ; (1b) σ m n m = F g m over Γ Fm ; (1c)
σ m = C m ε (u m ) in Ω m . (1d)
In the above equations, ε (u m ) denote the infinitesimal strain tensors, σ m the Cauchy stress tensors and C m the Hooke tensors. n 11 and n 2 represent the outward unit normals to Ω 11 and Ω 2 , respectively. . To complete the formulation of the boundary value problem, the interface condition has to be added. The writing of this condition depends on the considered domain Ω m :

u + -u -= 0 over Γ ; (2a) σ + n + + σ -n -= 0 over Γ ; (2b)
where subscript + is used to denote a variable of subdomain Ω m whereas subscript -refers to a variable of the neighboring subdomain. It ensures kinematic compatibility and equilibrium of the tractions, respectively, along the coupling interface Γ between the subdomains.

Non-symmetric Nitsche coupling formulation

As stated in the introduction, the non-symmetric Nitsche method has recently reemerged for the weak enforcement of boundary and coupling conditions in isogeometric immersed methods (see, e.g., [START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF] that provides a detailed review of the non-symmetric Nitsche method applied to the Poisson problem). In this part, we rewrite the derivation of the non-symmetric Nitsche formulation for the weak coupling of multiple domains in the context of linear elasticity. The obtained formulation constitutes the reference coupling of the present work.

In order to do so, let us start by defining the functional spaces U m and V m over subdomain Ω m that will contain the displacement solution and test functions respectively:

U m = u m ∈ H 1 (Ω m ) d , u m | Γu m = u g m ; V m = v m ∈ H 1 (Ω m ) d , v m | Γu m = 0 . (3)
We also introduce the space for the stress solution and for the virtual stress field:

Σ m = τ m ∈ H (, Ω m ) , τ m = τ T m , (4)
where:

H (, Ω m ) = τ m ∈ L 2 (Ω m ) d×d such that div (τ m ) ∈ L 2 (Ω m ) d . (5)
We note that Eq. ( 1a) is enforced strongly by considering space U m . To derive a Nitsche-type method, the first step consists in imposing weakly the equilibrium and behavior of the structure. In order to do so, we multiply Eqs. ( 1d) and (1b) by C -1 m τ m with τ m ∈ Σ m and by v m ∈ V m , respectively, and we perform integration by parts on each subdomain Ω m . Making also use of the Neumann boundary condition (1c), we find:

Ωm τ m : C -1 m σ m dΩ = - Ωm div (τ m ) • u m dΩ + ∂Ωm τ m n m • u m dΓ ; (6a) Ωm ε (v m ) : σ m dΩ = Ωm v m • f g m dΩ + Γ Fm v m • F g m dΓ + Γ v m • σ m n m dΓ. (6b)
We then discretize Eqs. (6a) and (6b) by means of the finite element subspaces

U h m ⊂ U m , V h m ⊂ V m and Σ h m ⊂ Σ m .
We end up with the flux formulation [START_REF] Arnold | Unified analysis of discontinuous Galerkin methods for elliptic problems[END_REF]: find u h m and σ h m such that:

Ωm τ h m : C -1 m σ h m dΩ = - Ωm div τ h m • u h m dΩ + Γu m ∪Γ Fm τ h m n m • u h m dΓ + Γ τ h m n m • dΓ ; (7a) Ωm ε v h m : σ h m dΩ = Ωm v h m • f g m dΩ + Γ Fm v h m • F g m dΓ + Γ v h m • n m dΓ ; ( 7b
)
where the numerical fluxes and are the approximations of the displacement and stress, respectively, on the interface Γ. To obtain the non-symmetric Nitsche formulation, the following expressions for the numerical fluxes are considered:

= 3 2 u h + - 1 2 u h -; (8a) = 1 2 σ h + + σ h -. (8b)
To make use of these expressions, we now consider the following integration by parts formula:

- Ωm div τ h m • u h m dΩ = Ωm τ h m : ε u h m dΩ - ∂Ωm τ h m n m • u h m dΓ. (9)
Considering Eq. (7a), replacing its second term by the right terms of Eq. ( 9), choosing τ h m = C m ε v h m on its left term, and inserting the flux approximation (8a), we then get:

Ωm ε v h m : σ h m dΩ = Ωm ε v h m : C m ε u h m dΩ + Γ C m ε v h m n m • 1 2 u h + -u h -dΓ. ( 10
)
Finally, by making use of the flux approximation (8b) into (7b), relating the result to the left-hand side of [START_REF] Gupta | Analysis of three-dimensional fracture mechanics problems: A non-intrusive approach using a generalized finite element method[END_REF] and summing over the two subdomains Ω 11 and Ω 2 , we end up with the non-symmetric Nitsche coupling formulation: find

u h 11 , u h 2 ∈ U h 11 × U h 2 such that: m a m u h m , v h m + Γ {} n 11 • u h dΓ - Γ v h • σ h n 11 dΓ = m l m v h m , ∀ v h 11 , v h 2 ∈ V h 11 × V h 2 . ( 11
)
In the above formulation, the standard elastic bilinear forms a m and linear forms l m are recovered:

a m (u m , v m ) = Ωm ε (v m ) : C m ε (u m ) dΩ ; (12a) l m (v m ) = Ωm v m • f g m dΩ m + Γ Fm v m • F g m dΓ ; (12b)
and the jump and average operators are introduced, respectively, as follows:

u h = u h 11 -u h 2 ; ( 13a
)
σ h = 1 2 σ h 11 + σ h 2 | Γ = 1 2 C 11 ε u h 11 + C 2 ε u h 2 . ( 13b
) - m Ωm v m • (div(σ m ) + f g m ) dΩ + m Γ Fm v m • (σ m n m -F g m ) dΓ (14a) + Γ 1 2 n 11 • (u 11 -u 2 ) dΓ + Γ v 11 • 1 2 (σ 11 -σ 2 ) n 11 dΓ (14b) + Γ 1 2 n 11 • (u 11 -u 2 ) dΓ + Γ v 2 • 1 2 (σ 11 -σ 2 ) n 11 dΓ = 0. (14c)
As a consequence, there is no need for additional stabilization terms with associated parameters in contrast to symmetric Nitsche methods [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Nguyen | Nitsche's method for two and three dimensional NURBS patch coupling[END_REF][START_REF] Annavarapu | A robust Nitsche's formulation for interface problems[END_REF][START_REF] Jiang | A robust Nitsche's formulation for interface problems with spline-based finite elements[END_REF]. Moreover, the formulation is free of auxiliary fields, thereby circumventing the difficulty of choosing a suitable dual space (particularly in case of non-conforming couplings) when considering mixed-type methods [START_REF] Brivadis | Isogeometric Mortar methods[END_REF]. Reminding the case of an evolving local model within an underlying NURBS mesh, the attractiveness of the non-symmetric Nitsche method becomes clear: it provides an easy and robust coupling in every non-conforming situations that may be encountered.

Non-intrusive resolution

We now address the non-intrusive resolution of the non-symmetric Nitsche coupling. In order to do so, we need to extend the strategy proposed in the framework of Lagrange multiplier coupling [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF]. After clarifying the classical monolithic approach which essentially consists in the approach illustrated in Fig. 1(c), the proposed non-intrusive algorithm is constructed. Finally, a reformulation of the iterative procedure is performed to highlight the mechanical and numerical properties of the algorithm.

Monolithic resolution

For ease of understanding of the next developments, let us first note that Eq. ( 11) reads: find u h 11 and u h 2 such that:

a 11 u h 11 , v h 11 - Γ v h 11 • σ h n 11 dΓ + Γ 1 2 n 11 • u h dΓ = l 11 v h 11 , ; (15a)
a 2 u h 2 , v h 2 - Γ v h 2 • σ h n 2 dΓ - Γ 1 2 n 2 • u h dΓ = l 2 v h 2 , . (15b)
Resolution [START_REF] Marussig | A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects[END_REF] of the coupling problem constitutes the classical monolithic approach: the resulting multiscale model of Fig. 1(c) is computed directly using a single direct solver. This strategy can be characterized as intrusive. Indeed, it is important to note that the global stiffness operator a 11 depends at this stage on the interface Γ or, in other words, on the shape of the local domain Ω 2 . Therefore, if domain Ω 2 is expected to evolve, not only the local operator a 2 but also the global operator a 11 have to be fully re-built and re-factorized.

Construction of a non-intrusive resolution algorithm

Rather than directly solving system (15), we start by considering the following asymmetric algorithm, where problems over Ω 11 and problems over Ω 2 are alternatively solved until convergence. For the n th iteration, starting with initial guesses u h (0)

11 and u h (0)

2

, we look for u h (n)

11 and u h (n) 2 such that:

1. Resolution of a problem over Ω 11 :

a 11 u h (n) 11 , v h 11 = l 11 v h 11 + Γ v h 11 • σ h (n-1) n 11 dΓ - Γ 1 2 n 11 • u h (n-1) dΓ. ( 16
)
2. Resolution of a problem over Ω 2 :

a 2 u h (n) 2 , v h 2 - Γ v h 2 • 1 2 σ h (n) 2 n 2 dΓ + Γ 1 2 n 2 • u h (n) 2 dΓ = l 2 v h 2 (17a) + Γ v h 2 • 1 2 σ h (n) 11 n 2 dΓ + Γ 1 2 n 2 • u h (n) 11 dΓ. ( 17b
)
Remark 1. A mechanical interpretation of the resulting iterative strategy will be given later in the paper (see section 3.3). Nevertheless, one may already notice that formulation (17) corresponds to the non-symmetric Nitsche weak imposition of the coupling conditions for subdomain Ω 2 only. Indeed, formulation (17) can be obtained by applying the procedure of section 2.2 to subdomain Ω 2 (only the sum over the two subdomains is not performed at the end). Another way of realizing this is to investigate the consistency of the formulation.

Considering the continuous version of (17), integrating by parts, and bringing all terms on the left-hand side, we find:

- Ω2 v 2 • div(σ (n) 2 ) + f g 2 dΩ + Γ F 2 v 2 • σ (n) 2 n 2 -F g 2 dΓ (18a) + Γ 1 2 n 2 • u (n) 2 -u (n) 11 dΓ + Γ v 2 • 1 2 σ (n) 2 -σ (n) 11 n 2 dΓ = 0. ( 18b
)
The above equation clearly shows that, at each iteration of the proposed algorithm, the interface term (14c) is recovered from the initial coupling formulation.

We now switch to the matrix form of problems ( 16) and ( 17) for ease of reading with respect to the implementation of the method. To this end, let us introduce the finite element basis functions N 3)), respectively. Denoting by u 11 and u 2 the associated degrees of freedom vector, the fixed point ( 16)-( 17) reads: for the n th iteration, starting with initial guesses u

K 11 u (n) 11 = f 11 -K N 11 T -K N 11 u (n-1) 11 -K N 21 T -K N 12 u (n-1) 2 . ( 19
)
2. Resolution over Ω 2 :

K 2 + K N 22 T -K N 22 u (n) 2 = f 2 -K N 12 T -K N 21 u (n) 11 . (20)
Operators K 11 (respectively F 11 ) and K 2 (resp. F 2 ) are the classical stiffness matrices (resp. vector forces) associated to subdomains Ω 11 and Ω 2 . K N is the standard Nitsche coupling operator (see, e.g., [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Nguyen | Nitsche's method for two and three dimensional NURBS patch coupling[END_REF]) defined as follows:

K N = K N 11 K N 12 K N 21 K N 22 , (21)
where the different blocks read:

K N 11 = Γ - 1 2 D 11 B 11 T N 11 dΓ ; (22a) K N 12 = Γ 1 2 D 11 D 11 T N 2 dΓ ; (22b) K N 21 = Γ - 1 2 D 2 D 2 T N 11 dΓ ; (22c) K N 22 = Γ 1 2 D 2 D 2 T N 2 dΓ. ( 22d
)
In the above equations, N 11 and N 2 represent the standard shape function matrices, B 11 and B 2 are the standard strain-displacement operators and, D 11 and D 2 model the Hooke constitutive law. In addition, is introduced to perform the product between the stress tensor and the outward unit normal to Ω 11 (see [START_REF] Nguyen | Nitsche's method for two and three dimensional NURBS patch coupling[END_REF][START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF] for more details regarding the construction of such operators).

Even if two problems are now solved separately, procedure ( 19)-( 20) is still intrusive because only the contributions over Ω 11 of the global model are involved. To go further, we make use of the continuous prolongation of the global solution from Ω 11 to Ω 12 and we apply the additivity of the integral with respect to domain Ω 1 = Ω 11 ∪ Ω 12 ∪ Γ. Denoting by u 1 the degrees of freedom vector associated to domain Ω 1 , this gives us:

K 1 u 1 = K 11 u 1 + K 12 u 1 . (23)
In the above equation, K 11 and K 12 are the extensions to Ω 1 of the classical stiffness matrices K 11 and K 12 related to subdomains Ω 11 and Ω 2 , respectively. They formally contain the classical stiffness operators and are padded with zeros to make them the same dimension of u 1 . As well, we define f 1 = f 11 + f 12 the load vector associated to domain Ω 1 . Equality ( 23) is used to expand problem [START_REF] Schillinger | The Finite Cell Method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF] from Ω 11 to Ω 1 . We finally end up with the following algorithm: for the n th iteration, starting with initial guesses u such that:

1. Resolution over Ω 1 :

K 1 u (n) 1 = f 1 -K N 11 T -K N 11 u (n-1) 1 -K N 21 T -K N 12 u (n-1) 2 + r (n-1) 12 . ( 24
)
2. Resolution over Ω 2 :

K 2 + K N 22 T -K N 22 u (n) 2 = f 2 -K N 12 T -K N 21 u (n) 1 . (25)
We note that the Nitsche operator K N simply consists of the prolongation of former operator K N from Ω 11

to Ω 1 . r 12 is introduced to denote the discrete reaction forces at Γ of the global model in the covered part Ω 12 . It reads at iteration n -1:

r (n-1) 12 = K 12 u (n-1) 1 -f 12 . (26)
Note that in expression [START_REF] Whitcomb | Iterative global/local finite element analysis[END_REF], only the degrees of freedom concerned with interface Γ are not zero. . Thanks to this prolongation, the whole stiffness matrix of the global NURBS model is now considered without any modification. In this sense, the iterative strategy ( 24)-( 25) is said to be non-intrusive. We emphasize that the global stiffness operator can be assembled and factorized only once (in a pre-processing step), and system (24) remains symmetric and well-conditioned (regardless of the evolution of the local shape). Futhermore, the global and local problem being solved alternatively and the interaction between the two models being restricted to interface Γ, the formalism offers the possibility to couple different numerical codes with few implementation effort. The price to pay is the number of iterations of the fixed point, that may appear quite important in some challenging situations. For more efficiency in such situations, we will see in forthcoming section 3.3 that acceleration techniques, such as based on an Aitken's Delta Squared method or a Quasi-Newton method [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF], can be applied to the present coupling. Numerical experiments to account for this point will also be carried out in Section 4.

Remark 2. It has to be noticed that the computation of the reaction forces r 12 (see Eq. ( 26)) requires to perform a volume integration to be consistent with the discrete form of the problem (see, e.g., [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF] for more information regarding this point). As a consequence, specific quadrature rules need to be implemented to integrate over pieces of global knot-span elements. In order to do so, a large amount of techniques, mostly taken from immersed boundary methods, may be used: for instance, the standard sub-triangulation technique in the context of X-FEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], the technique used in the NURBS Enhanced FEM [START_REF] Sevilla | NURBS enhanced finite element method (NEFEM)[END_REF], the recursive quadrature approach applied in the FCM [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Schillinger | The Finite Cell Method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF][START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF], or more recent geometrically faithfull quadratures [START_REF] Legrain | A NURBS enhanced extended finite element approach for unfitted CAD analysis[END_REF][START_REF] Nagy | On the numerical integration of trimmed isogeometric elements[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary methods[END_REF][START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF][START_REF] Stavrev | Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method[END_REF][START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF]. In this work, for simplicity and robustness, we perform the recursive quadrature approach of FCM.

Remark 3. We emphasize that the terminology "non-intrusive" used here only characterizes the numerical solver but not the geometrical construction (and parameterization) of the local model, which may not be a trivial task.

Incremental formulation, mechanical interpretation and convergence acceleration

In this section, another incremental formulation of the proposed non-intrusive algorithm is derived for a better understanding of the numerical procedure. It may be noted that the emphasis is made on the mechanical interpretation of the formulation, but not on the rigorous theoretical analysis of the algorithm.

To start with, let us denote by r N the Nitsche-based discrete reaction forces at Γ applied to the global model in Eq. [START_REF] Stavrev | Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method[END_REF]. We have at iteration n -1:

r N (n-1) = K N 11 T -K N 11 u (n-1) 1 -K N 21 T -K N 12 u (n-1) 2 . ( 27
)
Then, let us notice the following equality that corresponds to the discrete equilibrium of the global model at iteration n -1:

K 1 u (n-1) 1 = f 1 + r (n-1) 11 + r (n-1) 12 , (28)
where r 11 characterizes the discrete reaction forces at Γ produced by the global model in subdomain Ω 11 :

r (n-1) 11 = K 11 u (n-1) 1 -f 11 . (29)
Now, performing the subtraction between Eqs. ( 24) and ( 28), we arrive at:

K 1 u (n) 1 -u (n-1) 1 = -r N (n-1) + r (n-1) 11 . ( 30
)
As a result, procedure ( 24)-( 25) can be rewritten as follows: for the n th iteration, starting with initial guess u (0)

1 , we look for u

(n) 1
such that:

u (n) 1 = u (n-1) 1 -K -1 1 g u (n-1) 1 , (31)
where the application g is defined as:

g u (n-1) 1 = r N (n-1) + r (n-1) 11 , (32)
provided the local solution u

(n-1) 2

in r N (n-1) is computed from u

(n-1) 1

through system ( 25) (expressed at iteration n -1).

With formulation [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF] in hand, one can realize that the iterative procedure simply consists in a fixed point aiming at solving g = 0. The following residual is thus considered to stop the algorithm:

η = ||r N + r 11 || ||f 1 || 2 + ||f 2 || 2 . ( 33
)
In the following, we denote this residual by the "non-intrusive residual". Moreover, writing the continuous version of Eq. ( 32), we literally find:

g u (n-1) 1 = Γ 1 2 n 11 • u (n-1) 1 -u (n-1) 2 dΓ - Γ v 1 • 1 2 σ (n-1) 1 + σ (n-1) 2 n 11 dΓ + Γ v 1 • σ (n-1) 1 n 11 dΓ, ( 34
)
which corresponds to the consistency term (14b) of the initial coupling formulation. Therefore, the proposed iterative strategy aims at ensuring the coupling conditions for subdomain Ω 11 in a non-symmetric Nitsche weak sense. In addition, reminding that the other interface consistency term (14c) is enforced at each iteration of the non-intrusive procedure through the local problem [START_REF] Schillinger | An isogeometric design-throughanalysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and T-spline CAD surfaces[END_REF] (see term (18b) that corresponds the coupling conditions related to domain Ω 2 ), we ensure that if the proposed algorithm is convergent, all the interface terms of the initial coupling formulation ( 14) are recovered. In other words, if the developed non-intrusive algorithm is convergent, it leads to the same solution as the one computed from the monolithic approach (15). This numerical procedure is original and its mechanical interpretation is quite different from the existing non-intrusive strategies that aims at recovering the interface forces equilibrium provided that the Dirichlet conditions (over Γ) are transmitted, at each iteration, form the global to the local problem through the use of a Lagrange field (see Duval et al. [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF]).

Remark 4. Note that we need to compute r 11 (involving K 11 ) to evaluate residual [START_REF] Gosselet | Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems[END_REF]. The calculation is performed from the already computed stiffness K 12 (see Eq. ( 26)), i.e.:

K 11 = K 1 -K 12 .
Although different from the existing non-intrusive strategies in its mechanical interpretation, the algorithm can still be interpreted as a modified Newton method on g = 0 (see Eq. ( 31)). As a consequence, the application of the same acceleration techniques as the ones usually performed in the context of non-intrusive Lagrange multiplier coupling appears straightforward. In this work, we consider two acceleration methods: the Aitken dynamic relaxation and the Quasi-Newton update. The first technique, based upon the Aitken's Delta Squared formulation [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF][START_REF] Küttler | Fixed-point fluid-structure interaction solvers with dynamic relaxation[END_REF], proved to be very simple, cheap and of relatively good acceleration performance, especially when the global and local models are rather close (e.g., local mesh refinement, local plasticity). With a slightly more computational overhead, the second method enables to further speed up the algorithm and to make it convergent even in challenging situation. The inverse of the Symmetric Rank One update [START_REF] Kelley | Local Convergence of the Symmetric Rank-One Iteration[END_REF][START_REF] Khalfan | A Theoretical and Experimental Study of the Symmetric Rank-One Update[END_REF] (which converges towards ∇g -1 ) is computed in a non-intrusive fashion by means of a Sherman-Morrison formula. For more information regarding these accelerations techniques, the interested reader is advised to consult [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF].

Remark 5.

Numerical experiments

To assess the performance of the developed algorithm ( 24)-( 25), we now present a series of numerical experiments in two-dimensional linear elasticity. The plane stress assumption is considered for all test cases. The examples include different discretization scenarios that cover conforming and non-conforming couplings along straight and curved interfaces. The presentation is divided into two parts: first, the strategy is validated through three different benchmarks involving local mesh refinement, introduction of a geometrical detail, and coupling at a bi-material interface. In this part, the results are compared with those obtained from the existing standard approach based on Lagrange multipliers [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF], which demonstrates the performance of the proposed methodology. Then, the situation of an evolving local model is addressed. More precisely, a simple shape optimization problem is solved: the position of an inclusion within a global plate is optimized so as to provide maximal rigidity for the structure. The optimization algorithm calls upon the non-intrusive solver each time the objective function needs to be evaluated.

Validation of the non-intrusive coupling strategy

From here on, the mesh composed of N elements along the first length and M elements along the second length will be denoted N × M . The k-refinement strategy [START_REF] Cottrell | Studies of refinement and continuity in isogeometric structural analysis[END_REF] is performed so that the maximum level of continuity (i.e., C (p-1) for a polynomial degree p) is ensured across the interior knots. In the illustrations, we keep the notations introduced in Section 2; in particular, domain Ω 1 = Ω 11 ∪ Ω 12 ∪ Γ characterizes the global NURBS model and the local model of domain Ω 2 is expected to replace the global model in subdomain Ω 12 .

Curved beam subjected to end shear

The first example consists of a curved beam subjected to end shear. The problem, together with its global/local discretization, is illustrated in Fig. 3(a). A constant horizontal displacement of u 0 = 0.01 units is prescribed over the lower beam boundary. An analytical solution is available for the problem in [START_REF] Zienckewicz | The Finite Element Method -The Basis[END_REF]. In the upper half of the structure, the global NURBS model is meant to be replaced by a more refined (along the radial direction) NURBS model. Quadratic NURBS functions are considered for the discretization of both models. The situation of conforming meshes to be coupled is investigated here: the global/local interface is aligned with the edges of the global knot-span elements. As a result, no specific quadrature rules need to be implemented for the computation of r 12 (see remark 2 and Eq. ( 26)). However, the two meshes can be non-matching: the global and local mesh may not be aligned along the interface. To integrate over the interface, we divide the straight line into 1D sub-elements by combining the trace on the interface of the two meshes.

The results obtained by performing algorithm ( 24)-( 25) with the discretization of Fig. 3(a) are grouped in Figs. 3(b)-3(e). More precisely, the solution obtained once the iterative algorithm has converged is plotted in Figs. 3(b) and 3(c) in terms of displacement and Von Mises stress, respectively. We note that it is the coupled solution in Ω 11 ∪ Ω 2 that is mapped (the fictitious prolongation of the global solution over Ω 12 is not represented). The solution appears to be in a good agreement with [START_REF] Zienckewicz | The Finite Element Method -The Basis[END_REF]. In particular, there is no visible error concentration around the coupling interface Γ. Figs. 3(d) and 3(e) enable the convergence of the iterative strategy to be appreciated: first, in terms of the non-intrusive residual (Eq. ( 33)) and then in terms of the error on the displacements in the energy norm . No acceleration techniques are implemented here since the standard fixed point already appears efficient. The non-intrusive residual goes to zero and the error on the displacements reaches an asymptotic value (corresponding to the NURBS finite element approximation). This accounts for the convergence of the algorithm. More importantly, it may be actually noticed from Fig. 3(e) that it takes only 3 iterations to reach the stagnation of the NURBS finite element interpolation error. This means that it takes only 3 iterations for the coupling error to be negligible compared to that of the discretization. It seems that a criterion based on a coupling error of 10 -4 is a bit too demanding. To go further, a suitable error estimator should be implemented to stop the iterative algorithm at the right time [START_REF] Duval | Residual error based adaptive mesh refinement with the non-intrusive patch algorithm[END_REF]. Such an estimator, which should be extended to the case of NURBS, goes beyond the scope of this work, and we chose to stop the algorithm at 10 -4 in the following.

Then, the convergence of the method with the refinement of the mesh is studied. For this purpose, Figs. 4 and5 show the convergence curves for the error in the L 2 norm and in the energy norm, respectively. To refine the coupled solution, we consider the meshes indicated in Tab. 1. The situations of matching and non-matching meshes are investigated. For each approximation, the first mesh discretizes domain Ω 11 (this is the global mesh divided in half along the circumferential direction) and the second mesh is used for domain Ω 2 (this is the local mesh). For comparison purpose, the convergence curves of equivalent singlepatch solutions have been plotted in Figs. 4 and5. Also, the equivalent single-patch meshes are reported in Tab. 1. The convergence curves are plotted with respect to the equivalent mesh size normalized by the mesh size of the equivalent coarsest mesh. Finally, the standard non-intrusive solutions [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF] based on a simple Lagrange multiplier coupling have been added to the graphics. A pragmatic strategy is adopted to build the space for the Lagrange multiplier field: it is the trace along the interface of the local shape functions.

Total number of elements Single-patch mesh Matching meshes Non-matching meshes The present results confirm the behavior of the non-symmetric Nitsche coupling reported for the Poisson problem in [START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free imposition of weak boundary and coupling conditions in immersed finite elements[END_REF]: although optimal rates of convergence are achieved, the level of L 2 accuracy is reduced by a constant factor. Conversely, the corresponding convergence curves for the error in the energy norm are optimal: the rate of convergence and the error constant are equivalent to those of the single-patch discretization. These results indicate that the error regarding the derivatives of the primal variable is significantly smaller than the error due to the NURBS finite element approximation. This means that the accuracy of the derivative quantities is preserved with the proposed methodology. In the context of solid mechanics where, from an engineering point of view, the accuracy of the stresses is crucial for the structural design and optimization, this property is of primary importance while the reduced L 2 accuracy is acceptable. This encourages us to pursue in this direction.

(Ω 11 ∪ Ω 2 ) (Ω 11 ∪ Ω 2 ) 24 6 × 4 3 × 4 ∪ 3 × 4 3 × 3 ∪ 3 ×
Unlike its non-symmetric Nitsche counterpart, the standard non-intrusive procedure based on a Lagrange multiplier coupling exhibits optimal results in the energy norm as well as in the L 2 norm on this test case. We emphasize here that we are considering a simple coupling situation of conforming meshes. As a consequence, the choice adopted for the approximation of the dual variable appears suitable in this particular situation, which gives the full accuracy of the resulting mixed method. As will be illustrated in the next examples, this property of the Lagrange multiplier method will be deteriorated when addressing truly non-conforming interfaces.

Infinite plate with a circular hole

In the next step, we study the performance of the proposed methodology for the introduction of a geometrical detail in a given NURBS model. As an example, we consider the popular test case of an infinite plate with a circular hole under in-plane tension. The geometry, material, boundary conditions and the analytical stress solution [START_REF] Sadd | Elasticity, Theory, Applications, and Numerics[END_REF] are given in Fig. 6. The hole constitutes the geometrical detail to be incorporated in the initial NURBS square plate. The discretization of the problem following the developed non-intrusive strategy is illustrated in Fig. 7(a). A regular B-Spline mesh is used for the plate without the hole (domain Ω 1 ) and a circular refined ring-shaped NURBS domain is constructed, as the local model, for an accurate representation of the stress concentration around the hole. Therefore, the case of non-conforming coupling along a curved interface is this time investigated. The recursive quadrature approach [START_REF] Schillinger | The Finite Cell Method: A review in the context of higher-order structural analysis of CAD and image-based geometric models[END_REF] applied in the FCM is implemented to integrate over intersected elements for the evaluation of r 12 (see Figs. The results obtained with algorithm ( 24)-( 25) and the discretization of and 7(e) in terms of stress and non-intrusive convergence, respectively. Removing the smooth non-physical fictitious prolongation in Ω 12 , and replacing part of it by the local solution in Ω 2 , the resulting global/local solution appears smooth and to be in good agreement with results obtained using NURBS fitted discretizations [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry, and mesh refinement[END_REF][START_REF] Cottrell | Isogeometric analysis: Toward Integration of CAD and FEA[END_REF]. Regarding the convergence of the proposed algorithm, the Quasi-Newton acceleration technique seems to be necessary on this test case. This may be expected here since the stiffness gap between domains Ω 12 and Ω 2 is significant. Taking the advantage of the Quasi-Newton update, we are able to make the number of iterations relatively low: a residual of 10 -4 is obtained in about 20 iterations.

The finite element convergence of the method when applied to this non-conforming scenario is then investigated (see Figs. 8 and9). 2 nd , 3 rd and 4 th -order NURBS shape functions are considered for both of the models. The different meshes of Fig. 8 Unlike what was experienced for the curved beam problem of the previous section, the Lagrange multiplier approach appears to fail to reproduce the optimal convergence rates, especially when increasing the polynomial order of the NURBS shape functions. We remind that the performance of such mixed methods highly depends on the compatibility between the approximation subspaces considered for the primal and the dual variables. As a result, the reason for this poor accuracy may be due to the non-conforming situation encountered which rules out many choices for the dual subspace. Conversely, without any additional specific treatment, the proposed non-symmetric Nitsche version seems to achieve the optimal rates of convergence, which demonstrates that the methodology does not interfere with the accuracy of the NURBS functions regardless of the coupling scenario encountered. Note that for the finest quartic coupling discretization, the error level is so low that it may be deteriorated by rounding errors.

For completeness, the jump of the traction forces across Γ with respect to the parametric coordinate along the interface is plotted in Fig. 10 in case of a Nistche-based as well as a Lagrange multiplier based non-intrusive strategy. The second quartic discretization is considered. We observe that the Nitsche based strategy enables to significantly reduce the numerical jump in order to reproduce Eq. (2b) whereas important oscillations appear in case of a Lagrange multiplier based procedure, which confirms the superiority of the developed strategy when addressing non-intrusive non-conforming coupling.

Plate with a central inclusion

With the next example, the situation of a non-conforming, curved, and bi-material interface is investigated, considering the modeling of a plate with a central inclusion subjected to constant in-plane tension (see Fig. 11(a)). The Young's modulus is chosen to be a hundred times larger for the inclusion than for the plate (E i = 100 × E p ). The Poisson's ratios are the same. Regarding the discretization, a regular quadratic B-Spline grid composed of 28×40 elements is used for the plate while a circular quadratic NURBS mesh with 48 (circumferential direction) ×12 (radial direction) elements is constructed to model the inclusion. For the integration over intersected entitites, we keep the quadrature approach and we refine the NURBS interface irrespective of the underlying models. The interest of such a test case involving an overlapping mesh both finer and stiffer than the underlying grid, is that it may display mesh locking, as shown in the context of classical finite elements (see, e.g., [START_REF] Bechét | A stable lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF][START_REF] Lew | A discontinuous-Galerkin based immersed boundary method[END_REF]). For comparison purpose, a refined conforming finite element solution is computed using the classical 4-node element implemented in the industrial software ABAQUS (see Fig. 11(b)).

The results in terms of Von Mises stress are given in Fig. 12. The Quasi-Newton acceleration technique was used to make the algorithm converge and about 50 iterations were necessary to handle such a challenging situation. Due to the contrast in Young's moduli at the coupling interface, the same treatment as proposed in [START_REF] Sanders | A Nitsche embedded mesh method[END_REF] for the symmetric Nitsche coupling in the context of standard immersed finite elements is applied to the proposed non-intrusive non-symmetric Nitsche coupling. The average operator (13b) is thus modified as follows to counterbalance the material gap at the interface:

σ h = γσ h 1 + (1 -γ)σ h 2 with γ = E 2 E 1 + E 2 . ( 35
)
Thanks to this precaution, we are able to accurately recover the physical stress state with the developed non-intrusive strategy (the stiffer behavior of the inclusion is well captured, see Fig. 12(a)) whereas spurious oscillatory stresses are encountered at the interface with the standard Lagrange multiplier based procedure due to mesh locking (see Fig. 12(b)). To better appreciate the phenomenon, the jump of the traction forces across Γ with respect to the parametric coordinate along the interface is also plotted in Fig. 13. These results confirm the attractiveness of the proposed non-intrusive method over the standard approach: it can be easily extended with a suitable average operator to properly handle severe bi-material interfaces.

Remark 6.

Application to a simple optimization problem

With the last example, we seek to demonstrate the ability of the proposed strategy to accurately and efficiently compute any evolution of a local model within a fixed global NURBS one. Such a feature opens the door to tackle various applications. One of them concerns the optimization of local entities that can arise in different contexts: for instance, for the shape optimization of holes [START_REF] Wall | Isogeometric structural shape optimization[END_REF][START_REF] Fubeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF], or inclusions, or stiffeners within stiffened panels [START_REF] Mulani | EBF3PanelOpt: An optimization framework for curvilinear blade-stiffened panels[END_REF], or even, to perform local topology optimization [START_REF] Seo | Isogeometric topology optimization using trimmed spline surfaces[END_REF][START_REF] Parvizian | Topology optimization using the finite cell method[END_REF][START_REF] Groen | Higher-order multi-resolution topology optimization using the finite cell method[END_REF]. In this section, we compute, as a preliminary example, a simple optimization problem involving the positioning of an inclusion within a global plate. The objective is to show the potential of the developed method for carrying out structural optimization and to motivate future research in this direction.

The considered structural shape optimization problem is illustrated in Fig. 14. For the numerical model, the same parameters as for the previous test case (see Fig. 11(a)) are taken expect that the Young's modulus of the inclusion is chosen to be, this time, a hundred times smaller than for the plate (E p = 10 3 and E i = 10). The design variables are the horizontal and vertical coordinate of the center of the inclusion: x c and y c , respectively. The objective function is the compliance. As an optimization algorithm, we make use of the black box FMINSEARCH available in MATLAB. This routine uses a gradient-free algorithm: it is the Nelder-Mead simplex algorithm as described in Lagarias et al. [START_REF] Lagarias | Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions[END_REF]. We admit that more sophisticated algorithms may be probably used for improving the optimization performance. Nevertheless, our interest here being to prove the efficiency of the developed non-intrusive solver rather than building an advanced The optimization algorithm calls upon the developed non-intrusive solver to get the displacement solution at each configuration encountered through the optimization process. A maximum number of iterations of 20 is prescribed for the non-intrusive iterative algorithm. The interest of using algorithm (24)-( 25) is clear in this situation: the stiffness operators K 1 and K

2 + K N 22 T -K N 22
are assembled and factorized once for all in a pre-processing step. Then, the only things to do to compute the successive coupled solutions encountered at each iteration of the optimization process are:

1. build the integration rule for the small immersed region Ω 12 ; Conversely, a direct solver would require to fully rebuild and recompute the resulting global/local problem [START_REF] Marussig | A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects[END_REF] (that may also appear ill-conditioned) at each of these steps. The results of the optimization using the non-intrusive solver are presented in Fig. 16. More precisely, Figs. 16(a) and 16(b) show the convergence of the objective function and of the design variables, respectively. The convergence criterion for the optimization algorithm relies on the stagnation of the objective function as well as of the design variables: it stops when it gets below 10 -4 . It can be observed that we are able to recover the expected optimized configuration through the optimization process, which accounts for the ability of the non-intrusive solver to adapt to any arbitrary non-conforming situations, in particular when the edge of an element of the local model approaches that of an element of the global model.

Conclusion

In this paper, we proposed a new non-intrusive global/local coupling algorithm. In contrast to the existing non-intrusive strategies that rely on a Lagrange multiplier coupling, the developed iterative procedure makes use of the non-symmetric Nitsche approach. Recenty, the non-symmetric variant of the Nitsche method was successfully applied to the weak enforcement of boundary and interface conditions in non-boundaryfitted discretizations [START_REF] Schillinger | The non-symmetric Nitsche method for the parameter-free 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]. This work further consolidated this technique by establishing its application in the context of global/local isogeometric structural analysis. In this field, the non-symmetric Nitsche coupling gave the full potential of the non-intrusive approach: it offered the opportunity to address simply and efficiently all the non-conforming global/local coupling configurations that were encountered when integrating a local region within an underlying global NURBS model. The reason for this is the robustness and simplicity of the coupling. Indeed, the coupling formulation is (1) free of auxiliary fields which enables to circumvent the difficulty of choosing a suitable dual space (particularly in case of non-conforming couplings), and ( 2) is intrinsically stable (in a weak sense), thus eliminating the need for additional stabilization terms with appropriate parameters (whose determination can be delicate depending on the cutting scenarios).

The extension of the non-intrusive technique resulted in an original iterative process that aims at recovering the interface conditions (in a non-symmetric Nitsche weak sense) for the global model, provided that the coupling conditions viewed by the local model are enforced at each iteration of the algorithm. Although its mechanical interpretation is rather different form that of the existing non-intrusive strategies, the numerical procedure could still be interpreted as a modified Newton method so that, it was possible to apply efficient acceleration procedures in order to maintain a limited number of iterations (a few dozen) even in challenging situations. More precisely, the performance of the developed methodology over the existing NURBS one was numerically demonstrated through a range of two-dimensional elastic benchmarks involving conforming and non-conforming couplings along straight, curved, and possibly bi-material interfaces. Our results clearly indicated the superiority of the proposed approach when addressing any arbitrary non-conforming situation: the non-symmetric Nitsche version always achieved optimal accuracy while, due to a choice for the dual space, the existing Lagrange version appeared to fail to reproduce optimal behavior in some situations. In addition, mesh locking was reported when coupling at a severe bi-material interface with the Lagrange multiplier based method whereas the Nitsche counterpart handled successfully such interfaces thanks to the use of a suitable averaging operator.

Given its attractive properties, the developed method appears particularly suitable to compute any evolution of a local model within a fixed global NURBS one. From an engineering point of view, this appears of great interest since in many applications it is necessary to solve a sequence of similar problems, in a multiresolution process, where only the behavior at a local scale needs to be updated (e.g., to carry out the shape optimization of local entities [START_REF] Fubeder | Fundamental aspects of shape optimization in the context of isogeometric analysis[END_REF][START_REF] Mulani | EBF3PanelOpt: An optimization framework for curvilinear blade-stiffened panels[END_REF], or to perform the identification of local mechanical parameters [START_REF] Passieux | Multiscale displacement field measurement using digital image correlation: Application to the identification of elastic properties[END_REF][START_REF] Réthoré | Automatic crack tip detection and stress intensity factors estimation of curved cracks from digital images[END_REF], or to model crack propagation [START_REF] Passieux | Local/global non-intrusive crack propagation simulation using multigrid XFEM solver[END_REF], or the expansion of a plastic zone [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF], a damage zone [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF], etc). To show the potential of the method in this context, a simple optimization problem involving the positioning of an inclusion within a global plate was carried out. The optimization algorithm called upon the nonintrusive solver each time the objective function needed to be evaluated, which facilitated the process and ensured computational time saving. The success of this study opens the door for the development of effective optimization tools, based on the proposed non-intrusive solver, to handle more realistic applications (e.g., for the shape optimization of stiffeners within stiffened panels [START_REF] Mulani | EBF3PanelOpt: An optimization framework for curvilinear blade-stiffened panels[END_REF]). This motivates future research in this direction and, more generally, encourages to make use of such an algorithm when performing multiresolution processes.

  Interfarce exchanges(d) Approach 3: non-intrusive coupling.

Figure 1 :

 1 Figure 1: Example of a global/local NURBS problem. The global model over subdomain Ω 12 is replaced by the finer local model of domain Ω 2 through interface Γ, which enables to integrate a geometrical detail (a hole) within the initial NURBS model (top). Three different strategies may be applied to numerically solve the global/local problem (bottom).

  2(a) (right)). The substitution of the local model within the global one is achieved through interface Γ. The resulting global/local problem to be solved concerns the coupling of two overlapping NURBS meshes (multi-domain problem in Ω 11 ∪ Ω 2 ∪ Γ), the global solution in Ω 12 being discarded (see Fig. 2(b)).

F 2 (

 2 a) Description of the initial global/local NURBS probem.

11 2 (

 112 b) Resulting coupling problem to be solved.

Figure 2 :

 2 Figure 2: Reference mechanical coupling problem.

1 .

 1 Resolution over Ω 11 :

  Displacement field (magnitude).

  Convergence of the algorithm: non-intrusive residual.

  Convergence of the algorithm: error in energy norm.

Figure 3 :

 3 Figure 3: Global/local non-intrusive Nitsche analysis of the curved beam problem (NURBS mesh of quadratic 24 (circumferential direction) ×12 (radial direction) elements for Ω 1 (12 × 12 for Ω 11 ), and NURBS mesh of quadratic 12 × 20 elements for Ω 2 ).

Figure 4 :

 4 Figure 4: Convergence of the error in the L 2 norm for the curved beam problem.

  Non-matching meshes.

Figure 5 :

 5 Figure 5: Convergence of the error in the energy norm for the curved beam problem.

  7(b) and 7(c)). To ensure accuracy, we employ 6 levels of quadrature sub-cells. This precaution holds for all remaining examples of the paper. To integrate over the non-conforming interface, we refine the NURBS curve irrespective of the underlying global and local meshes.

Figure 6 :

 6 Figure 6: Infinite plate with a hole: description and data of the problem.

2 (

 2 Fig 7(a) are shown in Figs. 7(d)12 1 a) Discretization of the coupled problem. (b) Adaptive sub-cells for integration over Ω 12 . Normal stress σxx.

  acceleration (e) Convergence of the iterative strategy: non-intrusive residual.

Figure 7 :

 7 Figure 7: Global/local non-intrusive Nitsche analysis of the infinite plate with a hole (rectangular B-spline mesh of cubic 6 × 6 elements for Ω 1 , recursive Gaussian quadrature with 6 levels of sub-cells for Ω 12 , and circular NURBS mesh of cubic 4 × 4 elements for Ω 2 ).

  are built for the numerical study. They correspond to several uniform refinements starting from a mesh of 3 × 3 B-spline elements for the global model and 2 × 2 NURBS elements for the local model. The corresponding convergence curves obtained with the error in the energy norm are plotted in Figs. 9(a) and 9(b) for the proposed Nitsche-based strategy and for the standard Lagrange multiplier based strategy, respectively. For the Lagrange multiplier coupling, we keep the choice of a dual space constructed by means of the trace along the interface of the local shape functions.

Figure 8 :

 8 Figure 8: Sequence of meshes considered for the convergence study.

  Non-intrusive Lagrange strategy.

Figure 9 :

 9 Figure 9: Convergence curves in energy norm for the infinite plate with a hole.

  Figure 10:

  NURBS non-intrusive problem. (b) Reference FE mesh (from ABAQUS).

Figure 11 :

 11 Figure 11: Plate with a central inclusion: description and discretization of the problem.

  (a) Non-intrusive Nitsche strategy. (b) Non-intrusive Lagrange strategy. (c) Reference Abaqus solution.

Figure 12 :Figure 13 := 5 .

 12135 Figure 12: Plate with a central inclusion: Von Mises stress distribution.

Figure 14 :

 14 Figure 14: Simple optimization problem: the position of the inclusion (xc, yc) in the plate is optimized so as to provide maximal rigidity for the structure. The optimization algorithm calls upon the non-intrusive solver each time the objective function needs to be evaluated.

Figure 15 :

 15 Figure 15: Simple optimization problem: Von Mises stress distribution.

Figure 16 :

 16 Figure 16: Simple optimization problem: results.

  [START_REF] Passieux | Local/global non-intrusive crack propagation simulation using multigrid XFEM solver[END_REF] and Ω 2 are subjected to body forces f g

	Global model
			Local model
	u 11
	F 11	12	Substitution
			2	u 2
	11
	11 2 are associated to boundaries Γ F11 and Γ F2 11 and F g 2 , respectively. Furthermore, surface forces F g and f g
			4

  ,..,n2} that discretize the global and local model, respectively. In addition, we denote by N 11 C C∈{1,2,..,n11} the restricted part to subdomain Ω 11 of the shape functions of the global model. Following the principle of isoparametric elements, the basis N 11 C C∈{1,2,..,n11} and N 2 B B∈{1,2,..,n2} are used to build the above-mentioned finite element spaces U h 11 and U h 2 (see Eq. (

	1
	A A∈{1,2,..,n1}
	and N 2 B B∈{1,2

Table 1 :

 1 Meshes considered to study the convergence behavior of the curved beam problem.

	5