Marie Tirvaudey
email: tirvaude@insa-toulouse.fr

Ludovic Chamoin
email: chamoin@lmt.ens-cachan.fr

Robin Bouclier
email: robin.bouclier@math.univ-toulouse.fr

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

A posteriori error estimation and adaptivity in non-intrusive couplings between concurrent models

Keywords: Concurrent models, Non-intrusive coupling, A posteriori error estimation, Modeling error, Model adaptation, Goal-oriented approaches

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The fine analysis of localized complex phenomena has always been of major interest in simulation-based structural mechanics engineering. A typical case is aerospace engineering, in which local phenomena due to nonlinearities, heterogeneities, or geometric details are frequently studied over structures exhibiting various scales (from micrometer-sized composite fibers up to meter-sized aircraft components). For this purpose, and in order to avoid the use of numerically expensive physical models over the whole structure, a natural and usual practice is to introduce such high-fidelity models in a neighborhood of the regions of interest alone, at a local scale, while simpler and computationally cheaper surrogate models are used elsewhere at the global scale. These latter models may be obtained from homogenization or considering the linear range of the material behavior, for instance. The corresponding model reduction procedure leads to numerical simulations which can be performed with affordable computing resources in industrial applications. In this context, a wide variety of numerical methods, dedicated to multiscale and/or multi-model computing, have emerged. They can roughly be categorized in two main classes. The first class, mainly devoted to multiscale analysis with micro/macro models, consists in model enrichment by means of augmented approximation spaces (using finer meshes or specific enrichment functions) and superposition of micro/macro solutions. We may cite in this class :

-methods with enrichment based on a partition of unity (PUM) [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF] such as the Generalized Finite Element Method (GFEM) [START_REF] Strouboulis | The design and analysis of the generalized finite element method[END_REF][START_REF] Duarte | Analysis and applications of a generalized finite element method with global-local enrichment functions[END_REF] or the eXtended Finite Element Method (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] ; -other methods, such as adaptive localized Multiscale FEM (MsFEM) [START_REF] Hou | A multiscale finite element method for elliptic problems in composite materials and porous media[END_REF][START_REF] Efendiev | Multiscale Finite Element Methods: Theory and Applications[END_REF][START_REF] Chamoin | A posteriori error estimation and adaptive strategy for the control of MsFEM computations[END_REF], in which specific basis functions encode fine-scale details of the solution ; -methods with local correction, as performed in the Variational MultiScale method (VMS) [START_REF] Hughes | The variational multiscale method -A paradigm for computational mechanics[END_REF], the hierarchical modeling method [START_REF] Oden | Hierarchical modeling of heterogeneous solids[END_REF], multigrid methods [START_REF] Parsons | The multigrid method in solid mechanics -Part I: Algorithm description and behaviour[END_REF][START_REF] Rannou | A local multigrid XFEM strategy for 3D crack propagation[END_REF], the bridging scale method [START_REF] Wagner | Coupling of atomistic and continuum simulations using a bridging scale decomposition[END_REF], the chimera method [START_REF] Brezzi | Analysis of a chimera method[END_REF], numerical homogenization [START_REF] Engquist | Heterogeneous multiscale method: a general methodology for multiscale modeling[END_REF][START_REF] Feyel | A multilevel finite element method (FE2) to describe the response of highly non-linear structures using generalized continua[END_REF], or structural zooming with FE patches [START_REF] Glowinski | Finite element approximation of multi-scale elliptic problems using patches of elements[END_REF][START_REF] Picasso | Multiscale algorithm with patches of finite elements[END_REF]. Nevertheless, a major drawback of these model enrichment methods is that they can hardly be used in practical multiscale engineering activities due to their level of intrusiveness in existing commercial software. The second class of numerical methods, on which we particularly focus in this work, refers to model coupling methods with interface data transfers. These have received much interest with the emergence of new simulation trends in which several models, potentially coming from different software or physics, are used into parallel computations that are run on modern clusters. Among the wide list of coupling methods, and out of traditional sub-modeling (with one-way weak coupling) which is still a standard in industry [START_REF] Jara-Almonte | The specified boundary stiffness/force SBSF method for finite element subregion analysis[END_REF][START_REF] Voleti | Global-local analysis of large-scale composite structures using finite element methods[END_REF], we may refer to several advanced methods with strong coupling :

-improved iterative sub-modeling methods with global correction (taking into account the influence of local phenomena) or static condensation [START_REF] Hirai | An efficient zooming method for finite element analysis[END_REF][START_REF] Sun | A refined global-local finite element analysis method[END_REF][START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF] ; -the mortar method [START_REF] Belgacem | The mortar finite element method with Lagrange multipliers[END_REF][START_REF] Bernardi | Basics and some applications of the mortar element method[END_REF][START_REF] Brivadis | Isogeometric mortar methods[END_REF] enforcing weak equalities at the coupling interface by means of Lagrange multipliers ; -the Nitsche method [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF][START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF] ; -energy averaging methods with volume interface such as the Arlequin method [START_REF] Dhia | Problèmes mécaniques multi-échelles : la méthode Arlequin[END_REF][START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Prudhomme | Analysis of an averaging operator for atomic-tocontinuum coupling methods by the Arlequin approach[END_REF], the bridging domain method [START_REF] Belytschko | A bridging domain method for coupling continua with molecular dynamics[END_REF], or the MAAD method for atomic-to-continuum couplings [START_REF] Broughton | Concurrent coupling of length scales: methodology and application[END_REF]. Domain decomposition methods, such as well-known FETI [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF], BDD [START_REF] Mandel | Balancing domain decomposition[END_REF], FETI-DP [START_REF] Farhat | FETI-DP: a dual-primal unified FETI method -part I: A faster alternative to the two-level FETI method[END_REF], or mixed LA-TIN [START_REF] Ladevèze | A micro-macro and parallel computational strategy for highly heterogeneous structures[END_REF][START_REF] Daghia | A micro-meso computational strategy for the prediction of the damage and failure of laminates[END_REF], are also coupling methods based on Schwarz algorithms [START_REF] Lions | On the Schwarz method[END_REF] and are widely used in structural engineering [START_REF] Tallec | Domain decomposition methods in computational mechanics[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF]. Here again, all these model coupling approaches are intrusive as such, in the sense that they require quite deep modifications of FE solvers and software, as well as time-consuming meshing procedures, which is not always feasible in an industrial context.

More recently, a new and attractive class of model coupling methods referred to as non-intrusive localglobal coupling has emerged [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF], following pioneering ideas developed in [START_REF] Whitcomb | Iterative global-local finite element analysis[END_REF]. It consists in a substitution approach, with iterative solver, that enables local modifications of an existing finite element model (in terms of mesh refinement, introduction of local features related to the geometry or material behavior, . . .) while keeping the corresponding initial numerical operators unchanged at the global level. It defines a coarse global numerical model over the whole physical domain, in which geometry, connectivity, operators and solver are fixed (the initial factorized global matrix is thus conserved along the iterative coupling procedure), while local modeling evolutions are performed through a separate numerical model defined over local zones or patches. Interface data are then iteratively exchanged between these two models. Besides the increased flexibility with no global remeshing, this non-intrusive technique involving independent local-global solvers permits an easy merging of commercial software with any other specific simulation code dedicated to the modeling of complex phenomena of interest. Indeed, no modification of the commercial software is required and standard input/output specifications of such software can be fulfilled. Over the last decade, the non-intrusive local-global coupling method has been extensively applied and analyzed in many engineering situations exhibiting complex local phenomena. It was in particular implemented for problems with local plasticity [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Blanchard | Space/time global/local noninvasive coupling strategy: application to viscoplastic structures[END_REF], for crack propagation problems [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 a multigrid X-FEM solver[END_REF][START_REF] Gerasimov | A non-intrusive global/local approach applied to phase-field modeling of brittle fracture[END_REF], for the analysis of local uncertainties from a global deterministic operator [START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF][START_REF] Nouy | A multiscale method for semi-linear elliptic equations with localized uncertainties and non-linearities[END_REF], for 2D/3D couplings in thin composite panels with local stress concentration and debonding [START_REF] Guguin | Nonintrusive coupling of 3D ans 2D laminated composite models based on finite element 3D recovery[END_REF][START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF], for problems involving a NURBS definition of the domain shape and including local geometric details, fracture, or mesh refinement [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF][START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF], or for transient dynamics problems [START_REF] Bettinotti | A fast weakly intrusive multiscale method in explicit dynamics[END_REF][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]10]. It was also used in conjunction with domain decomposition techniques [START_REF] Duval | Non-intrusive coupling: recent advances and scalable domain decomposition[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF][START_REF] Gosselet | Non-invasive global-local coupling as a Schwarz domain decomposition method: acceleration and generalization[END_REF]. A global overview of the current capabilities of the non-intrusive local-global coupling method is available in [START_REF] Duval | Non-intrusive coupling: recent advances and scalable domain decomposition[END_REF].

In addition to model coupling, the present paper refers to a main challenge in simulation-based engineering, identified in the report of the NSF Simulation-Based Engineering Science panel [START_REF] Oden | Simulationbased Engineering Science: Revolutionizing Engineering Science through simulation[END_REF], and that deals with the certification of simulation models and methods. As any numerical method, the non-intrusive local-global coupling method is impacted by errors coming from various sources (discretization, modeling, iterative solution strategy) which need to be controlled in order to certify the numerical accuracy of the method and permit its transfer and intensive use in industrial activities related to robust design. In other words, several numerical parameters are inherent to the accuracy of the non-intrusive local-global coupling method (such as the number of iterations between local and global problems, the location and size of zones on which the local model should be implemented, or the mesh size used to perform computations with this latter model), and these parameters need to be carefully selected in order to get relevant simulation results in terms of output values used for decision-making. In the current literature, and contrary to other multiscale or multi-model methods [START_REF] Larson | Adaptive variational multiscale methods based on a posteriori error estimation: energy norm estimates for elliptic problems[END_REF][START_REF] Abdulle | A posteriori error analysis of the heterogeneous multiscale method for homogenization problems[END_REF][START_REF] Larsson | On two-scale adaptive FE analysis of micro-heterogeneous media with seamless scale-bridging[END_REF][START_REF] Jhurani | Multiscale modeling using goal-oriented adaptivity and numerical homogenization. Part 1: mathematical formulation and numerical results[END_REF][START_REF] Henning | An adaptive multiscale finite element method[END_REF][START_REF] Chamoin | A posteriori error estimation and adaptive strategy for the control of MsFEM computations[END_REF], there are very few works dealing with error estimation and adaptivity for non-intrusive local-global couplings. The very recent work detailed in [START_REF] Duval | Residual error based adaptive mesh refinement with the non-intrusive patch algorithm[END_REF] is probably the most advanced one in this context ; it constructs a cheap and global (i.e. in the energy norm) a posteriori error estimator based on an explicit residual technique. This estimator enables to control discretization and convergence (or algebraic) errors, and it may be used in practice to drive both mesh adaptation in the local model zone (supposed to have a fixed definition in [START_REF] Duval | Residual error based adaptive mesh refinement with the non-intrusive patch algorithm[END_REF]) and iteration stopping. Nevertheless, currently available verification tools for non-intrusive local-global couplings do not provide a quantitative error assessment with computable bounds. Moreover, they do not consider modeling and pollution errors which are major concerns in model coupling ; this is a drawback for robust design. In the present work, we wish to go much further in the certification of the non-intrusive local-global coupling method, so that the quality of simulation results can be fully controlled and the method can be confidently applied for industrial purposes [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF]. We thus develop advanced tools, in terms of fully computable a posteriori error estimator and indicators, in order to assess all error sources and drive effective adaptive procedures. These are constructed from verification approaches which have been extensively studied and implemented for a posteriori error estimation and mesh adaptation in the context of the finite element method [START_REF] Verfürth | A review of a posteriori error estimates and adaptive mesh-refinement techniques[END_REF][START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF][START_REF] Ladevèze | Mastering calculations in linear and nonlinear mechanics[END_REF]19] ; more specifically, we refer here to the residual functional. In this context, we develop verification tools which are dedicated to non-intrusive multiscale couplings and which require low implementation efforts. These tools can be seen as extensions of the pioneering ones developed in [START_REF] Oden | Estimation of modeling error in Computational Mechanics[END_REF] for modeling error estimation, and that were initially applied to heterogeneous materials in the context of the hierarchical modeling method [START_REF] Oden | Analysis and adaptive modeling of highly heterogeneous elastic structures[END_REF][START_REF] Oden | Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 1: Error estimates and adaptive algorithms[END_REF][START_REF] Vemaganti | Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 2: A computational environment for adaptive modeling of heterogeneous elastic solids[END_REF], before being implemented in other multiscale contexts [START_REF] Zohdi | A method of substructuring large-scale computational micromechanical problems[END_REF][START_REF] Oden | Multi-scale modeling of physical phenomena: adaptive control of models[END_REF][START_REF] Romkes | Multi-scale goal-oriented adaptive modeling of random heterogeneous materials[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF][START_REF] Chamoin | Control of modeling errors in the coupling of linear transport and diffusion models[END_REF][START_REF] Zaccardi | Error estimation and model adaptation for a stochastic-deterministic coupling method based on the Arlequin framework[END_REF]. The numerical strategy that we propose is also mainly devoted to the control of localized scalar quantities of interest defined from the local model, as such quantities represent critical outputs when resorting to non-intrusive local-global couplings. Consequently, we refer to the classical goal-oriented error estimation framework with the introduction of an adjoint problem [START_REF] Rannacher | A feedback approach to error control in finite element methods: application to linear elasticity[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF][START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF] whose solution acts as a filter to capture only part of the error that impacts the quantity of interest. In addition, error indicators are developed to separate contributions of each individual error source, including modeling error with pollution effects. They are computed at each step of the iterative and adaptive local-global coupling process. They thus feed a greedy adaptive algorithm that aims at automatically and iteratively meeting a given error tolerance with minimal computing effort, tuning at best the numerical model (mechanical model and mesh) as well as parameters of the coupling algorithm (number of local-global iterations). In particular, the local-global iterations are stopped when the convergence error (associated with goal-oriented unbalance at the coupling interface) becomes insignificant compared to other error contributions ; this is an alternative to classical stopping criteria based on the decrease of a norm of the interface residual, and it avoids useless and costly iterations which would not improve the quality of the solution outputs. Eventually, the verification procedure indicates where to put the final coupling interface (according to the level of modeling error), and which discretization should be used in the local model zone (according to the level of discretization error), so that a trade-off is obtained between solution accuracy and numerical cost. We emphasize that the non-intrusive feature of the local-global coupling substantially facilitates the implementation of the error estimation and adaptive procedures, as mesh refinement in the local model zone and modifications in the geometry of this zone can be performed independently of the global model. Moreover, it brings flexibility in the analysis of various scenarios for optimal and certified modeling. Throughout the paper, we assume that the continuous fine-scale model is free of error, that boundary conditions are perfectly known (no variability or uncertainty in their definition), that interface data are fully transmitted from one model to the other (using geometrically compatible meshes in particular), and that error sources coming from rounding of loading/geometry representation are negligible. For the sake of simplicity and clarity, the theoretical and subsequent numerical developments are conducted for linear models and local heterogeneities, even though extensions to nonlinear models are possible using linearized operators. Performance of the overall approach is shown on several two-dimensional numerical experiments with diversified coupling configurations and adaptivity sequences.

The remainder of the paper is organized as follows : in Section 2, we recall the basics of the non-intrusive coupling method and we introduce useful notations ; in Section 3, we develop goal-oriented verification tools based on the residual functional (i.e. weighted residual method) and devoted to the non-intrusive coupling context ; Section 4 focuses on the proposed adaptive strategy ; numerical results are reported in Section 5 ; finally, conclusions and prospects to this work are drawn in Section 6.

Non-intrusive coupling : context and basic implementation

Reference model

We consider a structural mechanics problem defined over a body occupying the closure of an open bounded domain Ω ⊂ R d (d = 1, 2 or 3 being the space dimension), with regular Lipschitz boundary ∂Ω (Fig. 1). We assume that a given displacement field u d is prescribed on a non-zero measured part ∂ u Ω ⊂ ∂Ω, while given traction forces F d are prescribed on the complementary part

∂ F Ω ⊂ ∂Ω, such that ∂ u Ω ∩ ∂ F Ω = ∅ and ∂ u Ω ∪ ∂ F Ω = ∂Ω.
A given body force field f d may also be active in Ω. In the following, and without loss of generality, we choose u d = 0 (homogeneous Dirichlet boundary conditions). Furthermore, we consider a quasi-static isothermal evolution with small perturbations regime. The material behavior is supposed to be described by a heterogeneous linear elasticity model, with possible fast variations of the material parameters. The mechanical problem then consists in finding the displacement-stress pair (u, σ σ) verifying :

u = 0 on ∂ u Ω (kinematic constraints)    div σ σ + f d = 0 in Ω σ σn = F d on ∂ F Ω (balance equations) σ σ = K (u) in Ω (constitutive relation) (1)
where n is the outward unit normal vector, (u) = 1 2 Grad(u) + Grad T (u) is the linearized strain tensor, and K is the heterogeneous linear Hooke operator. The weak form of this problem reads :

Find u ∈ V such that Ω K (u) : (v) = Ω f d • v + ∂ F Ω F d • v ∀v ∈ V (2)
where V = {v ∈ [H 1 (Ω)] d ; v = 0 on ∂ u Ω} is the appropriate functional space. We assume that in the previously considered model problem (2), phenomena of interest are localized in space. Consequently, a natural approach to reduce computational efforts consists in performing substructuring and restrict the use of a high-fidelity model to localized zones inside Ω, switching to a simpler model (in terms of material behavior, but also later in terms of mesh size) in the complementary part. We thus initially partition the physical domain Ω in two non-overlapping zones (see Fig. 2) :

Surrogate model and intrusive iterative coupling strategy

• a local zone Ω L ⊂ Ω, also denoted patch, that should encompass the support of the phenomena of interest to be analyzed. For the sake of simplicity, the zone Ω L is assumed here to be located strictly inside Ω. In this zone, the initial high-fidelity model (based on the constitutive operator K) is preserved ; • the complementary zone Ω 0 = Ω\Ω L in which a coarser concurrent model is implemented. It is defined by substituting the initial material behavior with an homogenized linear elastic behavior with Hooke's operator K 0 .

Remark 1 In the remainder of the paper, we assume that the initial position of the patch Ω L is a priori set from the support of phenomena of interest. An alternative, in the case where these phenomena are not identified, would consist in using a coarse model over the whole domain Ω and determine critical zones (by means of standard error estimates) from which the initial location of Ω L should be defined [START_REF] Picasso | Multiscale algorithm with patches of finite elements[END_REF].

The coupling problem is here formulated in a weak form using the Lagrange multipliers method. Introducing the interface Γ between zones Ω L and Ω 0 , the continuous coupling problem then consists in finding a global displacement field u G defined in Ω 0 , a local displacement field u L defined in Ω L , and a Lagrange multiplier field λ ∈ M (representing reaction forces on Γ), verifying :

• a global problem over Ω 0 :

Find u G ∈ V 0 such that Ω0 K 0 (u G ) : (v G ) = Ω0 f d • v G + ∂ F Ω F d • v G - Γ λ • v G ∀v G ∈ V 0 (3) with V 0 = {v ∈ [H 1 (Ω 0 )] d ; v = 0 on ∂ u Ω} ; • a local problem over Ω L : Find u L ∈ V L such that Ω L K (u L ) : (v L ) = Ω L f d • v L + Γ λ • v L ∀v L ∈ V L (4)
with

V L = {v ∈ [H 1 (Ω L )] d } ; • a continuity condition on Γ : Γ (u L -u G ) • µ = 0 ∀µ ∈ M (5)
This formulation naturally ensures the kinematic compatibility between global and local displacements and the balance of tractions [START_REF] Hansbo | A Lagrangian multiplier method for the FE solution of elliptic interface problems using non-matching grids[END_REF] on the interface Γ.

Using a discretization method with FE spaces V H 0 ⊂ V 0 (defined from a coarse partition τ H of Ω 0 ), V h L ⊂ V L (defined from an independent and usually finer partition τ h of Ω L ), and M h ⊂ M (i.e. the trace space defined from τ h ), the algebraic formulation of the above problem reads :

     K 0 0 C T G 0 K L -C T L C G -C L 0           U G U L Λ      =      F 0 F L 0      (6)
where U G , U L , and Λ are nodal value vectors of discretized fields u H G , u h L , and λ h , respectively, K 0 and K L are stiffness matrices in Ω 0 and Ω L , respectively, and C G and C L are coupling mortar operators. In practice, and in order to conform with domain decomposition techniques and parallel computing, the previous coupling problem is not solved in a monolithic way but rather by means of an iterative Dirichlet-Neumann solver. To do so, an asymmetric local-global algorithm with alternated interface data transfer is introduced. After initializing λ (0) = 0 (zero interface reaction), the continuous problem at iteration n consists in finding (u

(n) G , u (n) L , λ (n) ) ∈ V 0 × V L × M verifying
• a global problem over Ω 0 , with given Neumann boundary conditions on Γ, providing u

(n) G : Ω0 K 0 (u (n) G ) : (v G ) = Ω0 f d • v G + ∂ F Ω F d • v G - Γ λ (n-1) • v G ∀v G ∈ V 0 (7)
• a local problem over Ω L with given Dirichlet boundary conditions on Γ, providing (u

(n) L , λ (n) ) : u (n) L|Γ = u (n) G|Γ Ω L K (u (n) L ) : (v L ) - Γ λ (n) • v L = Ω L f d • v L ∀v L ∈ V L (8)
The corresponding algebraic formulation reads :

K 0 U (n) G = F 0 -C T G Λ (n-1) ;   K L -C T L -C L 0     U (n) L Λ (n)   =   F L -C G U (n) G   (9)
Remark 2 We assume here that meshes τ H and τ h are geometrically conforming even though they do not match on Γ, i.e. the interface is aligned with the edges of the local and global elements. As a result, the continuity of displacements can be enforced exactly (while traction equilibrium is enforced weakly on the interface approximation space) and a general mortar method [START_REF] Bernardi | Basics and some applications of the mortar element method[END_REF] is used to transfer interface data between local and global problems. The numerical experiments reported in Section 5 are performed in this context. In the more general case of a non-conforming interface, the transfer would require special attention in the implementation process, evaluating reaction forces with suitable quadrature rules, as performed in [START_REF] Bouclier | Local enrichment of NURBS patches using a non-intrusive coupling strategy: geometric details, local refinement, inclusion, fracture[END_REF] for NURBS geometry representations. Alternative matching conditions have also been introduced in the literature for non-intrusive couplings, such as these based on a more regular Mortar method [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], on a Nitsche method [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF], or on the use of a transition mesh to address topology changes between models [START_REF] Guguin | Nonintrusive coupling of 3D ans 2D laminated composite models based on finite element 3D recovery[END_REF][START_REF] Guguin | On the computation of plate assemblies using realistic 3D joint model: a non-intrusive approach[END_REF][START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF].

Non-intrusive coupling strategy

A drawback of the previous intrusive coupling technique is that the stiffness matrix K 0 , that depends on the geometrical definition of Ω 0 (and thus Ω L ), should be computed for each particular configuration of the local zone Ω L . Indeed, it requires the construction of a global mesh which is conforming with the potentially complex geometry of Ω 0 . Consequently, remeshing and new factorization of K 0 are necessary each time the location or shape of Ω L is changed (e.g. in case of crack propagation or optimization of local entities). This appears to be much time consuming, in particular for large domains with many dofs involved, and in a multi-query context. To circumvent this issue and enhance the numerical efficiency, the key idea of the non-intrusive local-global coupling strategy is to modify the global problem (3), defining the support of its solution u G over the whole domain Ω. In order to define the new global problem, the homogenized linear elasticity behavior is fictively prolongated to Ω L . Consequently, using additivity of the integral over

Ω 0 ∪ Ω L , the initial global problem (3) is recast as : Find u G ∈ V such that Ω K 0 (u G ) : (v G ) = Ω0 f d • v G + ∂ F Ω F d • v G - Γ λ • v G + Ω L K 0 (u G ) : (v G ) ∀v G ∈ V (10
) where n Ω L is the outward unit normal vector of Ω L . We emphasize that the corresponding solution u G , even though defined over the whole domain Ω, is usually non-physical in Ω L and irrelevant to analyze local phenomena of interest correctly (all the more so when u G is approximated using a coarse mesh). Furthermore, it is not unique in Ω L and a specific solution is in practice selected from the choice of the initialization in the coupling algorithm. Nevertheless, these do not represent issues as u G is eventually replaced by the local fine-scale solution u L in Ω L , so that there is no impact on the local-global solution.

Using the new discretization space V H , obtained from a coarse mesh τ H defined over the whole domain Ω, the non-intrusive procedure leads to the following change regarding the algebraic formulation of the global problem :

K 0 U G = F 0 -C T G Λ =⇒ K 0Ω U G = F 0 -C T G Λ + K 0L U G = F 0Ω -C T G Λ + R LG (11)
where K 0Ω (resp. K 0L ) is the invertible stiffness matrix computed over the whole domain Ω (resp. over the subdomain Ω L ) using the smooth linear operator K 0 , while

R LG = K 0L U G -F 0L
is the discrete interface reaction forces coming from the fictitious part of the global model, computed in practice from volume integrals.

Introducing again an iterative Dirichlet-Neumann solution scheme (fixed point algorithm), the nonintrusive local-global coupling method consists in finding, at each iteration n of the process and after initializing u

(0) G = 0 and λ (0) = 0, the set (u (n) G , u (n) L , λ (n) ) ∈ V × V L × M verifying • a global problem over Ω, with given internal reaction forces on Γ, providing u (n) G : Ω K 0 (u (n) G ) : (v G ) = Ω0 f d • v G + ∂ F Ω F d • v G - Γ λ (n-1) • v G + Ω L K 0 (u (n-1) G ) : (v G ) ∀v G ∈ V (12)
• a local problem over Ω L with given Dirichlet boundary conditions on Γ, providing (u

(n) L , λ (n) ) : u (n) L|Γ = u (n) G|Γ Ω L K (u (n) L ) : (v L ) - Γ λ (n) • v L = Ω L f d • v L ∀v L ∈ V L (13)
The corresponding algebraic formulation reads :

K 0Ω U (n) G = F 0 -C T G Λ (n-1) + K 0L U (n-1) G ;   K L -C T L -C L 0     U (n) L Λ (n)   =   F L -C G U (n) G   = F 0Ω -C T G Λ (n-1) + R (n-1)
LG [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF] We point out that the global stiffness matrix K 0Ω , as well as the global force vector F 0Ω , are fixed independently of the local model parameters (position and shape of Ω L , mesh size used in τ h ). They correspond to quantities that would be initially computed considering a smooth behavior over the whole structure, i.e. without any analysis of local complex phenomena, and using a global coarse mesh. The global stiffness operator is therefore assembled and factorized only once.

Essentially, the non-intrusive coupling technique thus consists in alternating between local calculations over Ω L with prescribed displacements on the coupling interface, and global correction calculations over the whole domain Ω which include inner corrective loads (in terms of equilibrium residual, i.e. reaction forces mismatch) in order to reduce the imbalance between concurrent models. A sketch of the associated local-global algorithm is given in Fig. 3.

Two independent numerical codes may be used to perform the local and global calculations. It can be shown that, under some conditions (e.g. multiscale elliptic problem, or local model operator not stiffer than the global model one, which is the usual case in practical applications), the solution to the fixed point ( 12)-( 13) converges to the solution to the initial coupling problem (3)-( 5). We refer to [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable domain decomposition[END_REF] for a review on these aspects, based on a global reformulation of the non-intrusive local-global coupling strategy that can be interpreted as a quasi-Newton algorithm on reaction force equilibrium.

Remark 3 The number of solver iterations is the price to pay in the non-intrusive coupling method compared with the intrusive version. Classically, a norm on the interface residual is used as a convergence indicator and stopping criterion. Nevertheless, we mention that convergence acceleration techniques can be used in this framework, such as the dynamic Aitken relaxation [START_REF] Irons | A version of the Aitken accelerator for computer iteration[END_REF][START_REF] Duval | Non-intrusive coupling: recent advances and scalable domain decomposition[END_REF], or the update of the global operator (without factorizing it again) using the symmetric rank one (SR1) update [START_REF] Conn | Convergence of quasi-Newton matrices generated by the symmetric rank one update[END_REF] and/or the Shermann-Morison and Woodbury formulas [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF]. Moreover, mixed interface conditions may be considered between local and global models [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF][START_REF] Oumaziz | A parallel non-invasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF]. All these techniques will not be implemented in the present work.

Remark 4 Basically, the non-intrusive local-global coupling framework is merely seen as a behavior substitution in Ω L (numerical zoom), starting from an initial smooth behavior (with operator K 0 ) defined over the whole domain Ω. We choose here to adopt another vision, deriving the coupling problem from an initial reference model in which the complex material behavior (with operator K) is introduced everywhere in Ω. This enables to have a consistent definition of the reference solution, from which error measures will be later defined.

Error definition

From the previous non-intrusive approach, and using discretization with meshes τ H and τ h , an approximate continuous local-global displacement field u hH(n) LG ∈ V can be recovered at each iteration n of the process. It is constructed as :

u hH(n) LG =    u h(n) L in Ω L u H(n) G in Ω 0 (15)
However, it should be noticed that the corresponding local-global stress field σ σ hH(n) LG , defined as :

σ σ hH(n) LG =    σ σ h(n) L = K (u h(n) L ) in Ω L σ σ H(n) G = K 0 (u H(n) G
) in Ω 0 [START_REF] Brivadis | Isogeometric mortar methods[END_REF] does not respect equilibrium in any weak sense (before convergence) across the interface Γ.

From the local-global displacement field u

hH(n) LG ∈ V, an error field e hH(n) LG = u -u hH(n)
LG ∈ V can be defined. Such an error field describes the discrepancy between the exact solution to the reference problem (2) and the approximate local-global solution at hand. Several scalar measures of the error field may then be used. In this paper, we will consider the measure Q(e hH(n) LG

) defined from a given linear quantity of interest Q. Nevertheless, global measures could also be used, such as the measure in the energy norm

• = Ω K (•) : (•).
In the non-intrusive coupling framework, error sources are of three types :

• modeling error due to the use of a surrogate model in Ω 0 , associated with a smooth material operator K 0 and a fixed (i.e., not adaptive) coarse mesh τ H . It may generate pollution effects when dealing with the accuracy of quantities of interest defined inside Ω L . The amplitude of this error source can be reduced by increasing the size of the critical zone Ω L , and it vanishes when Ω L = Ω ; • discretization error due to the use of a mesh τ h in order to approximate the solution of the local problem [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]. The amplitude of this error source can be reduced by decreasing the mesh size h in τ h , and it vanishes when h goes to zero ; • convergence error due to the use of an iterative local-global algorithm. The amplitude of this error source can be reduced by increasing the number of local-global iterations, and it vanishes when n tends to +∞.

In practical applications of the non-intrusive coupling method, Ω L and τ h are usually defined empirically from the a priori user experience, without any quantitative assessment of associated modeling and discretization errors. In addition, the convergence of the local-global iterative algorithm is classically controlled using stopping criteria (or convergence indicators) based on the magnitude of a norm on the interface residual. This procedure may be very pessimistic and may mobilize unnecessary computing resource, as : (i) the error tolerance on outputs of interest may be fulfilled even though the full local-global solution has not converged, so that the iterative algorithm could be stopped earlier without sacrificing the accuracy on these outputs ; (ii) the convergence error, even large, may rapidly become negligible compared to other error sources, so that further iterations become useless to decrease the overall error. Consequently, it is of interest to design tools that provide for a quantitative assessment of error measures as well as individual error contributions coming from various sources. Such tools could then be effectively used to drive an automated adaptive algorithm that optimally defines Ω L , τ h , and the required number of iterations (for a prescribed error tolerance), so that numerical performance in terms of computational cost is substantially enhanced. This is the topic of the following sections, in which we develop fully computable error estimators and indicators.

Verification tools based on weighted equilibrium residuals

Weak forms and residual functional

The reference problem (2) can be recast as :

find u ∈ V such that a(u, v) = l(v) ∀v ∈ V (17) with a(u, v) = Ω K (u) : (v) ; l(v) = Ω f d • v + ∂ F Ω F d • v (18)
Based on this weak form, we introduce the residual functional R :

V × V → R : R(w, v) = l(v) -a(w, v) (19)
which will be used for error estimation. The property [START_REF] Broughton | Concurrent coupling of length scales: methodology and application[END_REF] directly yields R(u, v) = 0 for any v ∈ V.

Further introducing the following notations :

a L (u, v) = Ω L K (u) : (v) ; a 0Ω (u, v) = Ω K 0 (u) : (v) ; a 0L (u, v) = Ω L K 0 (u) : (v) b Γ (λ, u) = Γ λ • u ; l L (v) = Ω L f d • v ; l 0 (v) = Ω0 f d • v + ∂ F Ω F d • v
(20) we also define additional weak forms for problems introduced in Section 2 :

• the continuous weak form of the local-global non-intrusive coupling (without iterative fixed-point solution scheme at this stage), coming from ( 4), [START_REF] Belgacem | The mortar finite element method with Lagrange multipliers[END_REF], and (10), reads :

find (u G , u L , λ) ∈ V × V L × M such that a 0Ω (u G , v G ) -a 0L (u G , v G ) + a L (u L , v L ) -b Γ (λ, v L -v G ) + b Γ (µ, u L -u G ) = l 0 (v G ) + l L (v L ) ∀(v G , v L , µ) ∈ V × V L × M (21)
or in a more condensed writing :

a LG ((u G , u L , λ), (v G , v L , µ)) = l LG (v G , v L , µ) ∀(v G , v L , µ) ∈ V × V L × M (22)
This provides the approximate solution u LG ∈ V defined as

u LG =    u L in Ω L u G in Ω 0 ;
• introducing the FE space V H , associated with coarse mesh τ H over Ω, the partially discretized version of [START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF] reads : find

(u H G , u L , λ) ∈ V H × V L × M such that a LG (u H G , u L , λ), (v H G , v L , µ) = l LG (v H G , v L , µ) ∀(v H G , v L , µ) ∈ V H × V L × M (23)
This provides the approximate solution u H

LG ∈ V defined as

u H LG =    u L in Ω L u H G in Ω 0 ;
• introducing the FE spaces V h L , and M h , associated with the mesh τ h used in Ω L , the fully discretized version of [START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF]

reads : find (u H G , u h L , λ h ) ∈ V H × V h L × M h such that a LG (u H G , u h L , λ h ), (v H G , v h L , µ h ) = l LG (v H G , v h L , µ h ) ∀(v H G , v h L , µ h ) ∈ V H × V h L × M h (24)
This provides the approximate solution u hH LG ∈ V defined as

u hH LG =    u h L in Ω L u H G in Ω 0 ;
• eventually, introducing the fixed-point scheme, the weak form at iteration n stemming from ( 12)-( 13) reads : find (u

H(n) G , u h(n) L , λ h(n) ) ∈ V H × V h L × M h such that a 0Ω (u H(n) G , v H G ) + a L (u h(n) L , v h L ) -b Γ (λ h(n) , v h L ) + b Γ (µ h , u h(n) L -u H(n) G ) = l 0 (v H G ) + l L (v h L ) + a 0L (u H(n-1) G , v H G ) -b Γ (λ h(n-1) , v H G ) ∀(v H G , v h L , µ h ) ∈ V H × V h L × M h (25)
or in a more condensed writing :

a (n)
LG (u [START_REF] Brezzi | Analysis of a chimera method[END_REF], and which is the available computed field when resorting to the non-intrusive local-global coupling framework.

H(n) G , u h(n) L , λ h(n) ), (v H G , v h L , µ h ) = l (n) LG (v H G , v h L , µ h ) ∀(v H G , v h L , µ h ) ∈ V H × V h L × M h (26) This provides the approximate solution u hH(n) LG ∈ V defined in
We emphasize again that global solutions u G (in ( 21)-( 22)) and u H G (in ( 23)-( 24)) are not unique in Ω L , but u 25)-( 26)) is. However, non-uniqueness is not an issue as : (i) global solutions are eventually replaced by fine-scale local solutions Ω L in the definition of the overall local-global solutions ; (ii) contributions of global solutions in Ω L vanish in the term a 0Ω (•, •) -a 0L (•, •) of the residual.

H(n) G (in (

The weighted residual method

In this section, we focus on the local error measure

Q(u)-Q(u hH(n) LG
) defined according to a given scalar quantity of interest Q(u) that is a specific (and usually fine-scale) feature of the solution u. We assume here that Q : V → R is linear, even though nonlinear quantities of interest could also be considered with minor changes (see [START_REF] Oden | Estimation of modeling error in Computational Mechanics[END_REF]). We also naturally assume that the quantity Q refers to features of u located in the initial configuration of Ω L . We then develop goal-oriented error estimation in a similar way as in [START_REF] Oden | Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 1: Error estimates and adaptive algorithms[END_REF][START_REF] Oden | Estimation of modeling error in Computational Mechanics[END_REF], and based on the definition of an adjoint problem.

Adjoint problem and error representation

We introduce the adjoint problem of ( 17), associated with Q, referring to [START_REF] Paraschivoiu | A posteriori finite element bounds for linear functional outputs of elliptic partial differential equations[END_REF][START_REF] Prudhomme | On goal-oriented error estimation for elliptic problems : application to the control of pointwise errors[END_REF][START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF] for an overview on the theoretical bases on this problem. It consists in finding u ∈ V such that

a(v, u) = a * ( u, v) = Q(v) ∀v ∈ V (27)
a * being constructed from the adjoint model operator. In the present case, the model operator is selfadjoint so that a * = a.

From the adjoint solution u, it is straightforward that for any approximation u app ∈ V of u, the error Q(u) -Q(u app ) can be represented as :

Q(u) -Q(u app ) = Q(u -u app ) = a(u -u app , u) = R(u app , u) (28)
Then introducing any approximation u app ∈ V of the adjoint solution u, the error representation also reads :

Q(u) -Q(u app ) = R(u app , u app ) + R(u app , u -u app ) (29)
Remark 5 The quantity of interest is usually defined in a global way by means of extraction functions. It is written under the form :

Q(u) = Ω σ σ Σ : (u) + Ω f Σ • u + ∂ F Ω F Σ • u + Ω K (u Σ ) : (u) - Ω u Σ • f ( 30
)
where σ σ Σ , f Σ , F Σ , and u Σ are extractors. These are defined explicitly or implicitly (depending on the quantity Q), and they can be mechanically interpreted as pre-stress, body force, traction force, and predisplacement, respectively, in the loading of the adjoint problem. The field u Σ , vanishing on ∂ F Ω, enables to extract components of the stress vector σ σ(u)n on ∂ u Ω (reaction forces).

Residual-based error estimator in the non-intrusive coupling method

Using the previous error representation, we now develop a computable error estimate on Q(u) when using the non-intrusive local-global coupling strategy. From [START_REF] Duval | Residual error based adaptive mesh refinement with the non-intrusive patch algorithm[END_REF], and noticing that u hH(n) LG ∈ V, we first write :

Q(u) -Q(u hH(n) LG ) = R(u hH(n) LG , u) (31)
For the term in the right-hand side to be computable, the adjoint solution u should be replaced by an approximate solution u app as described in [START_REF] Engquist | Heterogeneous multiscale method: a general methodology for multiscale modeling[END_REF]. Nevertheless, a relevant approximation should be computed so that R(u

hH(n) LG , u) ≈ R(u hH(n) LG , u app ) (i.e. R(u hH(n) LG
, uu app ) ≈ 0 can then be neglected in this case). This is the saturation assumption. In practice, this means that in order to catch the various error sources accurately, the approximation space used to compute u app should be richer than that used for u hH(n) LG

. Considering u app in the same approximation space as u hH(n) LG would lead to a poor error estimate of the error on Q.

Remark 6 A typical and well-known case illustrating the previous statement is the mere finite element approximation u f em of the solution u of (17) in a subspace

V f em ⊂ V. The Galerkin orthogonality R(u f em , v) = 0 for all v ∈ V f em indicates that the discretization error estimate R(u f em , u app ) is mea- ningless when u app is searched in V f em . Considering a richer space V + f em ⊂ V (with finer mesh size) to compute u app , the result R(u f em , u app ) = Q(u + f em ) -Q(u f em ) with u + f em ∈ V + f em also shows that the estimate catches all the error on Q(u) except the part Q(u) -Q(u + f em ).
For the considered non-intrusive local-global coupling method, enriching the approximation space for the solution of the adjoint problem means : (i) sufficiently enlarging the zone Ω L in which the original high-fidelity model is preserved (this enrichment is referred to with subscript "L + " in the following) ; (ii) sufficiently refining the mesh τ h used in this zone (this enrichment is referred to with superscript "h + " in the following) ; (iii) be sufficiently close to convergence in the iterative algorithm (referred to with superscript "∞" in the following). Consequently, after computing u h + H(∞) L + G ∈ V using an enriched nonintrusive local-global coupling method, with local part u

h + (∞) L + ∈ V h +
L + , an overall and fully computable error estimate of the error on Q reads :

η tot Q = R(u hH(n) LG , u h + H(∞) L + G ) ( 32
)
Remark 7 Due to the specific loading of the adjoint problem, which is concentrated inside Ω L , it is expected that the iterative local-global algorithm converges very fast when computing u

h + H(∞) L + G
.

Remark 8 In order to further reduce the computational cost without sacrificing too much the quality of the error estimate, it would be possible to approximate the residual functional R (initially defined from the reference model) considering the enriched approximation space used to solve the adjoint problem. Nevertheless, such an approximation does not prevent from projections between meshes for the computation of R(u

hH(n) LG , u h + H(∞) L + G
). This alternative will not be investigated here.

Adaptive strategy

Residual-based error indicators

The estimate (32) comprises all error sources. As described in Section 2.4, these are threefold : modeling, discretization, convergence. Indeed, introducing solution fields defined in Section 3.1, the error on Q can be split as :

Q(u) -Q(u hH(n) LG ) = Q(u) -Q(u H
LG )

∆ mod Q + Q(u H LG ) -Q(u hH LG ) ∆ dis Q + Q(u hH LG ) -Q(u hH(n) LG ) ∆ conv Q ( 33
)
where ∆ mod Q , ∆ dis Q , and ∆ conv Q correspond to modeling, discretization, and convergence parts of the error, respectively. We develop below some error indicators on each of these parts.

They are defined as follows :

• the indicator on convergence error, denoted by η conv Q , is constructed from a converged approximate adjoint solution u hH(∞) LG ∈ V with no enrichment in terms of mesh τ h and local zone Ω L used. It reads :

η conv Q = R LG (u hH(n) LG , u hH(∞) LG ) ( 34
)
where the residual R LG is defined from operators a LG and l LG associated with the local-global coupling problem (i.e., reference problem providing for u H

LG and u hH LG ) : R LG (u

hH(n) LG , u hH(∞) LG ) =l 0 ( u H(∞) G ) + l L ( u h(∞) L ) -a 0Ω (u H(n) G , u H(∞) G ) + a 0L (u H(n) G , u H(∞) G ) -a L (u h(n) L , u h(∞) L ) ( 35
)
The indicator is such that η conv • the indicator on discretization error, denoted by η dis Q , is constructed from a converged approximate solution u h + H(∞) LG ∈ V computed with a finer local mesh τ h + alone, while the shape of Ω L remains unchanged compared to that used for the computation of u hH(n) LG . It reads :

η dis Q = R LG (u hH(n) LG , u h + H(∞) LG ) -η conv Q (36)
and is such that η dis Q -→ h→h + ≈ 0. It should provide a relevant quantitative indication on the discretization error ∆ dis Q provided h + is small enough. • eventually, the indicator on modeling error, denoted by η mod Q , is constructed from an approximate solution u hH(∞) L + G ∈ V computed with a larger zone Ω L + alone, while the mesh τ h is unchanged compared to that used for the computation of u hH(n) LG . It reads :

η mod Q = R(u hH(n) LG , u hH(∞) L + G ) -η conv Q (37)
and is such that

η mod Q -→ Ω L →Ω L + ≈ 0.
It should provide a relevant quantitative indication on the modeling error provided Ω L + is large enough. An alternative construction of the indicator η mod Q , giving in practice slightly different values but decreasing the number of adjoint solutions, stems from the following (and still empirical) definition :

η mod Q = η tot Q -η conv Q -η dis Q (38)
Remark 9 It is worth noticing that numerical strategies which have to be implemented for the computation of the estimator η tot Q , as well as indicators η conv Q , η dis Q , and η mod Q , are in accordance with the non-intrusive framework. Indeed, the definitions of the enriched spaces which are used to compute the approximate adjoint solutions u

h + H(∞) L + G , u h + H(∞) LG
, and u hH(∞) L + G require modifications of V h L alone, while V H is kept unchanged. This can be easily performed using the non-intrusive coupling methodology. In addition, the non-intrusive framework applied to the solution of the adjoint problem enables to select specific error sources and analyze various modeling configurations in a suitable manner. By a flexible introduction of additional patches to Ω L + , located on some preselected zones (e.g. in the vicinity of geometrical details such as holes), the corresponding adjoint solution automatically filters targeted error sources due to orthogonality properties described in Section 3.2.2. Therefore, the critical phenomena that affect the accuracy on the quantity of interest, even though located far from the region over which the quantity of interest is defined (pollution effects), are easily detected. These phenomena would then need to be further modeled accurately, i.e. at the fine-scale level.

Adaptive algorithm

From the previously defined error estimator η tot Q and indicators η conv Q , η dis Q , and η mod Q , it is possible to set up a relevant adaptive algorithm in order to drive the non-intrusive coupling algorithm. We propose to use here a greedy algorithm closely related to those proposed in [START_REF] Oden | Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 1: Error estimates and adaptive algorithms[END_REF][START_REF] Vemaganti | Estimation of local modeling error and goal-oriented modeling of heterogeneous materials; Part 2: A computational environment for adaptive modeling of heterogeneous elastic solids[END_REF][START_REF] Oden | Multi-scale modeling of physical phenomena: adaptive control of models[END_REF][START_REF] Romkes | Multi-scale goal-oriented adaptive modeling of random heterogeneous materials[END_REF][START_REF] Bauman | Adaptive multiscale modeling of polymeric materials: Arlequin coupling and goals algorithms[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF][START_REF] Zaccardi | Error estimation and model adaptation for a stochastic-deterministic coupling method based on the Arlequin framework[END_REF] (so-called Goals algorithms). The approach, which refers to goal-oriented adaptivity, aims at automatically tuning the parameters of the local-global coupling method (shape of Ω L , mesh size in τ h , number of local-global iterations) in order to predict the quantity of interest Q within a preset error tolerance γ tol while optimizing the computational cost. This is achieved by generating a sequence of approximate solutions u (k) app so that for some integer k 0 , the overall error on Q satisfies :

|Q(u) -Q(u (k0) app )| ≤ γ tol |Q(u (k0) app )| (39)
At each iteration of the adaptive process, and before stopping the full adaptive algorithm when the error tolerance is met (quantitative information given by η tot Q ≤ γ tol |Q(u app )|), the goal is to reduce the major error source which is identified comparing indicators η conv Q , η dis Q , and η mod Q . Adaptations in discretization and modeling are conducted locally after decomposing the indicators over predefined subdomains in Ω L and Ω 0 , respectively. In practice, subdomains in Ω L are chosen as elements of τ h , while subdomains in Ω 0 are defined from elements of the coarse mesh τ H (even though larger subdomains could be used). This decomposition is possible by observing that indicators η dis Q and η mod Q correspond to residual terms defined from space integrals.

After initializing Ω L (as a neighborhood of the region over which the quantity of interest is defined) and τ h (with similar mesh size as for τ H ), and after specifying the error tolerance γ tol for the quantity of interest, the proposed adaptive algorithm reads as follows :

Compute the adjoint solution

u h + H(∞) L + G
(using an appropriate enriched space) ; 1. Set n = 1 ; 2. Solve the primal surrogate problem for u

hH(n) LG ; 3. Compute the estimate η tot Q ; 4. If |η tot Q /Q(u hH(n) LG
)| ≤ γ tol then STOP. Otherwise proceed to Step 5 ;

5. Compute solutions u hH(∞) L + G , u h + H(∞)
LG

, and indicators

η conv Q , η dis Q , and η mod Q : • if max(|η conv Q |, |η dis Q |, |η mod Q |) = |η conv Q |, increment n + 1 → n and go to Step 2 ; • if max(|η conv Q |, |η dis Q |, |η mod Q |) = |η dis Q |, decompose η dis Q
and locally refine τ h up to reaching

|η dis Q /Q(u hH(n) LG )| ≤ γ tol /3, then go to Step 0 ; • if max(|η conv Q |, |η dis Q |, |η mod Q |) = |η mod Q |, decompose η mod Q
and locally enlarge Ω L up to reaching

|η mod Q /Q(u hH(n) LG
)| ≤ γ tol /3, then go to Step 0.

This adaptive algorithm prevents from useless local-global iterations for the primal problem (when discretization or modeling error is larger than convergence error). It also indicates, at the end of the adaptive process, a suitable definition of Ω L and τ h for reaching the error tolerance.

Numerical results

In this section, we perform a series of 2D numerical experiments to show the performance of the proposed verification strategy. Several linear elasticity problems (with plane stress assumption) and quantities of interest are considered. In all cases, the local (fine-scale) model is initially preserved in a critical region of the structure that corresponds to an encompassing neighborhood of the region over which the quantity of interest is defined.

Quasi-1D slender structure in tension

In this first example, a simple slender structure under traction loading is considered ; the objective is to illustrate the splitting of error sources and the adaptive procedure. We assume that the beam-like structure, of unit length (L = 1) and width b = 0.1, is clamped on its left end and that a uniform traction load F = 1 is applied on its right end (Fig. 4). The beam is made of a linear isotropic elastic material, and we assume that the Young modulus E varies along the longitudinal coordinate. Indeed, its nominal value E 0 = 1 is weakened in some local regions in order to enforce the occurence of local phenomena that may need to be captured using a local-global coupling. The local weakening of the Young modulus is modeled in terms of Gaussian functions, and two configurations are considered for the evolution of E(x) as described in Fig. 5. In order to simplify notations in the remainder of this section, Case 1 and Case 2 will refer to the configurations shown in Fig. 5(a) and Fig. 5(b), respectively. The Poisson ratio is assumed to be constant and is set to ν = 0. Regarding the numerical solution, the global mesh τ H is made of 16 quadratic macro elements which are placed along the structure and the thickness of each element corresponds to that of the beam (i.e. there is only one element across the thickness for the mesh associated with the global model). The same transverse mesh size (one element across the thickness) applies for the local mesh τ h . The problem is solved using a non-intrusive coupling algorithm as detailed in Section 2 with :

-a global problem defined all along the structure Ω and in which a constant Young modulus E 0 = 1 is considered ; -a local problem defined in a sub-zone Ω L that corresponds to a subset of macro elements of τ H , and in which the variation of the Young modulus E is considered. A finer mesh τ h may also be considered in this zone if need be. An example of the associated coupling configuration is shown in Fig. 4.

The quantity of interest Q which is chosen to illustrate the adaptive procedure is the average of the longitudinal displacement u x in the subregion ω Q =]13/16, 15/16[×]0, b[ of the structure, that is

Q(u) = 1 |ω Q | ω Q u x .
The associated adjoint problem is also solved using a non-intrusive method so that the global adjoint model is not changed when the local adjoint model is modified (and consequently local error analysis can be easily performed). Its loading consists of a uniform horizontal body force f Σ = 8/b applied in ω Q . The initial local patch for the primal problem is placed in the subregion 13/16 ≤ x ≤ 15/16 (i.e. Ω L = ω Q ) ; it has the size of two macro elements (Fig. 6). The relative error tolerance γ tol is set to 5%. In order to perform Step 3 of the adaptive procedure and find the global error estimate η tot Q (as defined in (32)), a larger patch is needed to solve the adjoint problem and catch error sources. In order to reduce the computational time, this adjoint problem is in practice replaced with a set of adjoint problems with smaller patches that can be solved in parallel thanks to the non-intrusive coupling method. Even though the resulting combination of adjoint solutions does not coincide with the one that would be obtained with a monolithic numerical approach, it leads to similar results in terms of error estimation due to orthogonality properties mentioned in Section 3.2.2. To that extent, six adjoint problems are defined with a local zone Ω L + that contains a patch in the region ]13/16, 15/16[ and a second patch with varying locations depending on the adjoint problem number, as explained in Fig. 6(a). After solving each adjoint problem (with enough iterations to ensure the convergence of the local-global coupling, the global error estimate η tot Q after convergence of the primal local-global coupling can be evaluated for Cases 1 and 2 ; this estimate is defined as a sum of error contributions computed from the individual adjoint problems defined above. The distribution of the estimate per macro element is shown in Fig. 6(b) and Fig. 6(c). The main corresponding error source can be here interpreted as the modeling error source (as the mesh in the local zone is quite fine and the coupling algorithm has converged here). We observe that for Case 1 (Fig. 6(b)), the error is concentrated in the region ]11/16, 13/16[ whereas for Case 2 (Fig. 6(c)), it is rather concentrated in the region ]3/16, 7/16[. These results are in good agreement with the respective evolutions of the Young modulus E in both cases. Indeed the local zone of the primal problem being initially placed in the region ]13/16, 15/16[ (while taking E 0 = 1 in the remainder of the structure), modeling errors are large in the area where the exact value of E is far from E 0 , and they impact the quantity of interest. )|) are reported for each adaptive step. The number n of iterations in the local-global coupling algorithm is indicated on the graphs (red digits), and the evolution of the size of the local zone Ω L in the primal problem is also given at the bottom of the plots. This representation illustrates again that when the indicator on modeling error η mod Q is the highest, the local zone Ω L is extended (e.g. an additional patch is added at Step 1, and then enlarged at Step 3, for Case 2). For these examples, results show that discretization error in Ω L is low and has almost no influence on the quantity of interest. It is also interesting to note that, in both cases, there is no need to model the small decrease of the Young modulus for a tolerance of 5% ; only the highest peak in Fig. 5(a) and Fig. 5(b) needs to be taken into account. Nevertheless, further computations for a tolerance of 1% would require to model this small decrease and to refine the mesh in Ω L . Additionally, the comparison with the reference error (computed from an overkill solution and represented by the black curves in Fig. 7) enables to verify that η tot Q is a relevant estimate of the overall error on Q.

The conclusion of this study is that ideal coupling configurations (with respect to the given error tolerance γ tol = 5% on the quantity of interest) are :

-for Case 1, a local zone Ω L defined by a unique patch in the region ]11/16, 15/16[ with no specific refinement of τ h and 4 iterations (n = 4) in the coupling algorithm, according to Fig. 7(a) ; -for Case 2, a local zone Ω L defined by patches in regions ]3/16, 7/16[ and ]13/16, 15/16[, with no specific refinement of τ h and 3 iterations (n = 3) in the coupling algorithm, according to Fig. 7(b).

Remark 10 Keeping this simple example, other quantities of interest can further benefit from the adaptive procedure. Considering for instance the sub-zone ω Q =]13/16, 15/16[, we may consider :

-the average of the traction stress component σ xx in ω Q , the adjoint loading being then a pre-stress -the average of the traction strain component xx in ω Q , the adjoint loading being then a pre-stress

σ σ Σ = 1 |ω Q | K   1 0 0 0   = 1 |ω Q | . E 1 -ν 2   1 0 0 ν   in ω Q , or equivalently a body force f Σ = -∇ • σ σ Σ in ω Q and tractions σ σ Σ n on the boundary ∂ω Q of ω Q ;
σ σ Σ = 1 |ω Q |   1 0 0 0   in ω Q , or equivalently tractions σ σ Σ n on the boundary ∂ω Q of ω Q .
In these two cases, the algorithm indicates that no adaptive step is needed as stress and strain fields do not depend much on the upstream configuration. Consequently, the quantities of interest can be accurately computed with very little numerical effort.

2D plate in traction with local weakening inclusions

We now consider a square plate (size L × L with L = 1) in which localized weakening of the material stiffness is considered. The structure, represented in Fig. 8, is clamped on its left side and subjected to a uniform traction on its right side ; other boundaries are free. The global mesh τ H is made of 100 (10×10) first-order quadrangular elements. Local variations of the Young modulus E(x, y) take the form of Gaussian functions in five zones which are smaller than a macro element. They act as inclusions inside the material where the Young modulus is lower than its nominal value E 0 = 1. Two cases are here considered : (i) low weakening with minimal Young's modulus value E min = 0.45 ; (ii) high weakening with minimal Young's modulus value E min = 0.0026. The impacted macro elements are shown in Fig. 8(b), while the variation of E inside the structure is shown in Fig. 8(c) for a large contrast (i.e. high weakening). The Poisson ratio is fixed and set to ν = 0.3. The reference solutions (where the exact material behavior and a very fine mesh are considered over the whole domain Ω) are given in Fig. 9 for these two cases, in terms of component xx of the strain field. The quantity of interest is the average longitudinal displacement on the right edge x = L where the traction loading is applied. The goal of the adaptation procedure is to find the optimal configuration for the coupled problem regarding this quantity of interest, and with respect to a given error tolerance. This tolerance is set to γ tol = 0.5% (this value enables to detect the small impact of the modified Young modulus on the predicted value of the quantity of interest).

Analysis of the adaptive algorithm

We show here the adaptivity results when considering a small or large contrast in the local change of the Young modulus. The initial configuration of the coupled problem is such that the local zone Ω L is composed of elements 91 to 100 on which the traction loading is applied (see Fig. 8(b)) ; the discretization in this local zone (mesh τ h ) is similar to that used for the global mesh τ H . We first show results for a small modification of the Young modulus (E min = 0.45) in Fig. 10. Starting from the initial solution given in Fig. 10(a), where we observe effects of the clamping on left corners when using a coarse mesh, the adaptive procedure is performed. In Fig. 10(c), the values of the different relative estimator and indicators (i.e. normalized by the approximate value of the quantity of interest) are given at each adaptation step. These are

|η tot Q /Q(u hH(n) LG )|, |η conv Q /Q(u hH(n) LG )|, |η dis Q /Q(u hH(n) LG
)|, and

|η mod Q /Q(u hH(n) LG
)|. We observe that 3 adaptation steps are performed : the first two deal with model adaptation by adding two macro elements (first element 83 then element 89) in the primal local zone Ω L , and the last one increases the number n of iterations in the coupling algorithm. The error tolerance is then reached, and the obtained approximate local-global solution is given in Fig. 10(b). The discretization error is again very low, and the corresponding normalized indicator converges to zero very rapidly even though no refinement of τ h is performed. In this final coupling configuration : (i) only elements 83 and 89 (close to the right-hand side of the structure and relying on the fine-scale Young modulus) are added to Ω L ; (ii) two iterations are performed in the non-intrusive local-global coupling algorithm (i.e. n = 2) ; (iii) no mesh refinement is necessary in Ω L . Even though the approximate solution shown in Fig. 10(b) is quite different from that given in Fig. 9 (left), it is accurate enough for the prediction of the quantity of interest. In particular, some of the local variations of the Young modulus have very little impact on the quantity of interest and do not need to be represented. Remark 11 A similar example was run using an even smaller contrast (E min = 0.8). In this case, the adaptive procedure shows that the major error source, which is modeling error, is located in elements 1 and 10 which are corner elements where the Dirichlet boundary condition are applied and where large gradients occur (stress concentration). Two additional patches with refined mesh should thus be incorporated in Ω L in order to decrease the error on the quantity of interest effectively.

We now show results for a large modification of the Young modulus (E min = 0.0026). Starting from the same initial configuration of the coupled problem, the adaptive procedure is performed and results are shown in Fig. 11. The values of the different relative indicators are given at each adaptation step in Fig. 11(c). We observe that 7 adaptation steps are now required to reach the error tolerance : five are related to model adaptation by adding macro elements inside the local zone Ω L , and two steps increase the number of iterations in the coupling algorithm. In the first three steps, the error is mainly due to modeling and local evaluations indicate that elements 89, 83, and 58 should be successively included in the local model. The other macro elements (44 and 25) which are the support of a change in Young's modulus are added to the local zone Ω L in steps 5 and 7. Consequently, the algorithm indicates that the local zone should here cover all elements which contain an inclusion in order to reach the error tolerance on the quantity of interest. Nevertheless, the local weakening of the Young modulus do not need to be described precisely as no mesh refinement is required for τ h . We also notice (looking at Step 8 on Fig. 11(c)) that mesh refinement would be the next adaptation step if a lower error tolerance was chosen. The approximate solution obtained at the end of the adaptive process is reported in Fig. 11(b) ; it is interesting to see that it is sufficient to predict the quantity of interest accurately even though it remains a coarse approximation of the exact strain field shown in Fig. 9 (right).

Heterogeneous local variation of the Young modulus

We now consider a case where a zone on which the Young modulus is decreased has a larger area and impacts more than one macro element of the global mesh τ H . The associated configuration is shown in Fig. 12(a), and the contrast is such that E min = 0.45.

For this example, eleven adaptive steps are required to reach a tolerance γ tol = 1% on the quantity of interest, as shown in Fig. 12(b), and these are mostly related to model adaptation. In order to detail the adaptive process, we represent in Fig. 13 the spatial distribution of the indicator η mod Q on modeling error per macro element of τ H , the position of the local patches constituting Ω L (grey zones), and the macro element (in black) that was included in the local zone Ω L at the current adaptation step. At the end of the adaptive process, the configuration of the coupling problem is such that Ω L is made of the element set {91 -100, 58, 59, 57, 25, 83, 68, 48, 89} (elements are listed in the order they are included in Ω L ), 3 iterations are performed in the local-global coupling algorithm, and no refinement is needed.

Eventually, for this last configuration of the Young modulus distribution, we consider the control of the error on another quantity of interest Q. It is the average of the strain component xx in the macro element 68 (which is in the neighborhood of the large weakened zone). The local zone Ω L initially consists of macro elements 58 and 68. Applying the adaptive process for this quantity indicates that the main error sources are initially due to coupling iterations and local discretization so that the mesh τ h in the local zone Ω L needs to be refined in order to reach the tolerance γ tol = 2%. This tolerance, obtained after 4 iterations of the adaptive algorithm, also requires n = 3 local-global iterations but no extension of Ω L . We show in Fig. 14 several features of the goal-oriented adaptation strategy : the adjoint solution (that exhibits large localized gradients in the vicinity of the region of interest) is shown in Fig. 14(a), the evolution of error estimator and indicators along the adaptive process are given in Fig. 14

Bending plate with holes

In this final application, we consider a plate with a regular (periodic) distribution of 160 holes with constant radius r = 0.15, and submitted to a bending loading. The dimensions of the plate and boundary conditions are detailed in Fig. 15. The Young modulus is E = 1 and the Poisson ratio is ν = 0.3. The reference solution, in terms of yy component of the strain field, is given in Fig. 16(a) ; all the 160 holes are considered in this case, and an overkill computation is performed.

Using the local-global coupling framework, the solution is approximated considering : -a global model made of the plate without any hole and with an homogenized Young modulus E 0 = (1 -πr 2 )E (effective modulus obtained from a weighted average). The global mesh τ H used for this model is composed of 8×20 first-order quadrangle elements ; -a local model in a zone Ω L ⊂ Ω made of a set of patches, each patch representing a squared domain including a hole (Fig. 15(b)). The size of one local patch is 1×1 (that is, the size of a macro element), the Young modulus is E = 1, and the unstructured mesh is composed of first-order triangular elements.

In the adaptive process, we consider that the local mesh τ h is fine enough so that discretization error in Ω L is neglected. Consequently, the adaptive procedure aims at setting the optimal number n of iterations in the coupling procedure as well as selecting the holes which need to be represented using a patch (definition of the size of Ω L ). The quantity of interest is the average of the vertical displacement on the right edge (x = 20) of the structure. Naturally, the initial local zone Ω L is placed in the vicinity of this edge (see Fig. 18(a)). It is made of the layer of 8 macro elements of τ H which are connected to the right edge of Ω ; out of the 8 associated holes, the other holes are not represented. We show in Fig. 16(b) the approximate solution obtained when considering this coupling configuration.

Setting the error tolerance to γ tol = 2%, the adaptive process is performed, starting from the previous coupling configuration with one iteration in the local-global algorithm. Using the non-intrusive framework, potential critical zones are analyzed by placing an additional patch in each macro element of Ω 0 when solving the adjoint problem, in order to catch the associated error sources. The obtained set of adjoint problems is then solved in parallel and enables to detect zones where the modeling error is the highest.

Remark 12 For the computation of residuals, a fine projection grid is used and elements of this grid which are inside holes are marked and discarded in the computation.

The adaptation results are shown in Fig. 17, where the evolutions of the relative error estimator and indicators (in terms of iteration and modeling error sources) are shown along the adaptive process (Fig. 17(a)), as well as the final approximate coupled solution verifying the tolerance on the quantity of interest (Fig. 17(b)). This solution is obtained after 63 iterations of the adaptive process ; it requires n = 3 iterations in the coupling algorithm and an enlarged local zone Ω L taking into account pollution effects from the coarse global model. The adaptive process can be detailed as follows :

-from Step 1 to Step 16, modeling error is predominant so that patches are added to Ω L . The size of this zone at Step 16 is shown in Fig. 18(b) ; -at Step 17, an additional iteration is performed in the coupling algorithm ; -from Step 18 to Step 62, news patches are added to Ω L . The size of this zone at Step 62 is shown in Fig. 18(c) ; -at Step 63, an additional iteration is performed in the coupling algorithm and the error reaches the preset tolerance γ tol . We emphasize that we chose here to add only one patch to Ω L at each iteration of the adaptive process, so that the required number of iterations to reach the error tolerance is quite large. An alternative to decrease the number of iterations would be to add several patches in the same time (selecting macro elements of τ H in which modeling error is larger than a threshold). Nevertheless, optimality of the final coupling configuration may be lost with this procedure, and this is why we chose not to apply it in the paper.

Eventually we now consider as the quantity of interest the average of the stress component σ xx in a local zone ω Q ⊂ Ω that corresponds to the upper-left macro element of τ H . For this quantity, the loading of the adjoint problem consists of a pre-strain inside ω Q . Starting from an initial configuration where Ω L = ω Q , and setting an error tolerance γ tol = 2%, adaptation results are reported in Fig. 19. The evolutions of the relative error estimator and indicators along the adaptive process are shown in Fig. 19(a) ; they indicate that the preset error tolerance is reached after 4 iterations of the adaptive algorithm, with n = 2 local-global iterations, and with an enlarged local zone Ω L . The final configuration of Ω L is displayed in Fig. 19(b), and the final approximate local-global solution is shown in Fig. 19(c). It is interesting to notice that for this last case, only 3 holes out of 160 really need to be represented in order to reach the error tolerance on the quantity of interest, so that much computing resources and meshing effort can be saved.

Conclusions and prospects

We presented in this paper a new strategy to control the accuracy of local-global coupling strategies with respect to some quantities of interest. This strategy, based on residual functionals and adjointbased techniques, enables to compute right at the right cost. It defines a fully computable error estimate (quantitatively certifying the quality of the approximation) as well as error indicators which are used in an adaptation process. These indicators enable to split error sources between iterations (i.e. lack of convergence at the coupling interface), modeling, and discretization so that useless over-computations are avoided (e.g. the iterative solver is usually stopped before reaching convergence in terms of the usual interface equilibrium). It is important to notice that the strategy is made consistent with the nonintrusive framework of the coupling ; it can thus be performed when coupling two different codes, and local analyses for error sources (by adding local patches when solving the adjoint problem) can advantageously benefit from this non-intrusive framework. Consequently, the adjoint solution does not require prohibitive computing resources but is rather conducted by defining individual and manageable problems (that differ by the position of local patches) which can be all solved in parallel at a global cost similar to that of the primal coupling problem. Also, the proposed strategy can be extended to nonlinear problems using linearized operators, even though the error estimator and indicators may not be fully robust in some cases. It thus appears as an attractive numerical tool for practical engineering applications. Further developments should address the application to large structures with complex nonlinear material behaviors (such as damage), as well as the computation of robust (e.g. mathematically guaranteed) error bounds on quantities of interest. Another interesting study would consist in assessing the impact of an interface error (e.g. approximate data transfer coming from incompatible meshes) on outputs of interest. All these points will be the topics of forthcoming research works.

Figure 1 .

 1 Figure 1. The reference problem and its environment.

Figure 2 .

 2 Figure 2. Sub-structuring of the physical domain.

Figure 3 .

 3 Figure 3. Illustration of the non-intrusive strategy.

  should provide a quantitative indication on the convergence error ∆ conv Q , enabling to define a relevant stopping criterion for the local-global iterative solver.

Figure 4 .

 4 Figure 4. Initial configuration of the local-global coupling.

Figure 5 .

 5 Figure 5. Two examples of evolution of the Young modulus E along the longitudinal direction.

  (a) Different locations of a second patch in the adjoint solution.

  Error map for Case 1.

  Error map for Case 2.

Figure 6 .

 6 Figure 6. Adjoint solution strategy (left), and initial global error map for Case 1 (center) and Case 2 (right).

  (a) Local patch locations. (b) Local patch locations.

Figure 7 .

 7 Figure 7. Evolution of the error estimator and indicators at each adaptive step for Case 1 (left) and Case 2 (right). The reference error is computed using an overkill solution (no coupling) of the primal problem.

  (a) Considered structure with boundary conditions. (b) Position of the zones with degenerated Young modulus. The impacted macro elements, numbered {25, 44, 58, 83, 89}, are the same for the highly and slightly degenerated cases.(c) Evolution of the Young modulus inside the structure for the highly degenerated case.

Figure 8 .

 8 Figure 8. Considered problem with local variations of the Young modulus.

Figure 9 .

 9 Figure 9. Reference solution (component xx of the strain field) for low (left) and high (right) contrast in the Young modulus field.

  (a) Primal coupled solution at initial step. (b) Primal coupled solution after adaptation. (c) Evolution of error indicators along the adaptive process.

Figure 10 .

 10 Figure 10. Results for a low local variation of the Young modulus : (a) and (b) represent the field of the xx strain component at steps 0 and 4, respectively, of the adaptive process ; (c) gives the evolution of the different error indicators at each adaptation step.

  (b), while the final local mesh τ h and the final approximate local-global solution (requiring n = 3 local-global iterations) are shown in Fig. 14(c) and Fig. 14(d), respectively.

  (a) Primal coupled solution at initial step. (b) Primal coupled solution after adaptation. (c) Evolution of error indicators along the adaptive process.

Figure 11 .

 11 Figure 11. Results for a large local variation of the Young modulus : (a) and (b) represent the field of the xx strain component at steps 0 and 7, respectively, of the adaptive process ; (c) gives the evolution of the different error indicators at each adaptation step.

( a )

 a Larger weakening zone : the local decrease of the Young modulus also impacts the neighborhood of element 58. (b) Evolution of error indicators along the adaptive process.

Figure 12 .

 12 Figure 12. Local weakened zone inside the structure (a) and adaptivity results (b).

Figure 13 .

 13 Figure 13. Distribution of the indicator on modeling error at different steps of the adaptation procedure. In the top figures of each step, the local zone Ω L is in grey and the newly added elements in the local zone are in black.

  (a) Stress component σxx of the adjoint problem. (b) Evolution of error indicators along the adaptive process. (c) Final mesh used to approximate the solution in the local-global coupling process. (d) Primal coupled solution after adaptation.

Figure 14 .

 14 Figure 14. Results when considering as a quantity of interest the average of the strain component xx in the vicinity of the large weakened zone : influence on the adjoint solution (a), and final stress field in the structure after applying the adaptive algorithm.

  (a) Configuration of the studied problem. (b) Zoom on the local mesh used to represent each hole.

Figure 15 .

 15 Figure 15. Description of the bending plate problem. The reference geometry (a) is composed of 160 holes that may be each represented by a patch (b) in the numerical approximation.

  Approximate solution (initial coupling configuration).

Figure 16 .

 16 Figure 16. Map of the εyy strain component for the considered plate problem : (a) when all the holes are considered (reference solution, no coupling) ; (b) when a local-global coupling strategy is used with a local zone made of one layer of macro elements (holes are represented in this zone alone).

  Map of the yy strain component for the final coupling solution.

Figure 17 .

 17 Figure 17. Evolution of error quantities and final coupling configuration for the control on the plate with holes.

Figure 18 .

 18 Figure 18. Evolution of the local zone Ω L along the adaptation procedure. [10] Bettinotti O, Allix O, Perego U, Oancea V, Malherbe B. Simulation of delamination under impact using a global-local method in explicit dynamics. Finite Elements in Analysis and Design 2017; 125:1-13.

  (a) Evolution of error estimator and indicators. (b) Final location of Ω L and Ω 0 , with associated meshes, inside the global structure.

  Map of the yy strain component for the final local-global solution.

Figure 19 .

 19 Figure 19. Adaptation results when considering the average of σxx in the local zone ω Q as the quantity of interest. [19] Chamoin L., Díez P. (ed). Verifying calculations, forty years on: an overview of classical verification techniques for FEM simulations. SpringerBriefs, 2016.