Mickaël Duval

Alexei Lozinski

Jean-Charles C Passieux

Michel Salaün

Residual error based adaptive mesh refinement with the non-intrusive patch algorithm

Keywords: Multiscale, Finite Elements, Non-intrusive coupling, Error estimation, Mesh refinement

teaching and research institutions in France or abroad, or from public or private research centers.

Residual error based adaptive mesh refinement with the non-intrusive patch algorithm

Introduction

It exists a wide variety of applications in which a given finite element mesh may not be locally fine enough in order to take into account some localized phenomena. In such a case, when pre-processing with mesh adaptation is either not possible or hardly feasible, several numerical methods still allow to carry on an analysis guaranteeing a sufficient accuracy. For instance, if the vectorial space generated by the finite element discretization is too poor to well take into account local scale phenomena (crack for instance), then the model can be enriched using the eXtended Finite Element Method [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] or, in a more general context, the Generalized Finite Element Method [START_REF] Strouboulis | The generalized finite element method[END_REF][START_REF] Kim | A new generalized finite element method for two-scale simulations of propagating cohesive fractures in 3-D[END_REF]. Another efficient way for bridging different scales is to rely on multigrid algorithms [START_REF] Rannou | A local multigrid XFEM strategy for 3D crack propagation[END_REF][START_REF] Passieux | Direct estimation of generalized stress intensity factors using a three-scale concurrent multigrid XFEM[END_REF], which allow for relevant computations while keeping a reasonable computational cost. Nevertheless, setting up an enriched finite element model or making use of a multigrid solver requires to use an ad-hoc software. While remaining very efficient, such methods may not be suited to all situations, especially if one uses a software which does not support such features.

Then a flexible and efficient solution is to rely on an iterative algorithm with patches of finite elements [START_REF] Glowinski | Finite element approximation of multi-scale elliptic problems using patches of elements[END_REF][START_REF] Rezzonico | Multiscale algorithm with patches of finite elements[END_REF][START_REF] Pironneau | Numerical Zoom for localized Multiscales[END_REF], which is based upon the wider class of Schwarz algorithms [START_REF] Gander | Schwarz Methods over the Course of Time[END_REF]]. Such methods rely on separate finite elements models and solvers. A patch of finite elements is used to "zoom" the solution provided by a global scale model, without modifying it. Then, in its non-intrusive version [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], the patch algorithm can be used regardless of the used software and the underlying class of solver for each scale.

The main idea of this paper is to join together non-intrusive patch algorithm and mesh refinement. A major advantage is that it allows to confine mesh refinement to a restricted area, defined at the beginning of the procedure. Then, this area, our "patch", may be refined classically, in a h and/or p manner and several times, fully independently from the mesh, model and operators of the global domain, its coupling to the patch being driven by the non-intrusive algorithm.

In a rather similar way, [START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF]Bauman et al. 2008;Bernardi et al. 2005] present local/global coupling where the local model size is driven by a posteriori error estimates. Nevertheless, the underlying coupling methods are intrusive as they require to modify locally the global model. It requires for instance to remove some coarse elements, which implies renumbering, reassembly and re-factorization of the stiffness operator of the global model, each time the patch shape and mesh evolve. More precisely, the full scale FE model of complex industrial products takes often long time to build. In addition, it is developed in a particular FEA tool and, most of the time, it is very difficult to transfer the model to another software. So, the full scale FE model is generally the result of a trade-off between precision and computation time. In that respect, it cannot be easily remeshed or refined to improve its ability to analyse a complex localized phenomena if it was not initially intended to. The goal of non-intrusive coupling consists in improving, at the lowest cost, an existing and unmodifiable global FE model, the remeshing efforts being concentrated on the patch (see Fig. 1).

Finally, this paper is a first contribution to the analysis of a posteriori finite element error in the context of local/global non-intrusive coupling. Even though non-intrusive coupling is mostly relevant in a nonlinear context [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], the paper is restricted to the case of linear elasticity, where the theory can be developed rigorously.

Remark 1. Local/global coupling being relevant only in the case of localized complexities, an implicit assumption is also that the error source is localized. In this respect, another problem, which may arise, concerns pollution effects, when localized fine scale phenomena are affected by effects coming from large scale distances. This point will not be addressed here as it does not seem adapted to a local/global strategy.

So the paper is organized as follows. First is introduced the patch method and it is shown how it can be used in a non-intrusive way with the coupling algorithm. Then, an explicit residual based error estimator is proposed, which is classical in its construction, but however must be adapted to the multi-scale nature of the coupled problem. Such an estimator is used to build error map which is used to drive the local mesh refinement. In this respect, illustrations are given, arising from a structural analysis benchmark. Finally, we show how the convergence of the coupling algorithm can be monitored using error estimation.

Non-intrusive coupling and patch algorithm

The two-dimensional linearized elasticity problem

Let Ω be a bounded polygonal domain in R 2 and ∂Ω its boundary. We split this boundary into two disjoint parts: ∂Ω = ∂Ω D ∪ ∂Ω N . Furthermore, the solid body, which occupies domain Ω, is assumed to be fixed along ∂Ω D . Moreover, it is subject to body forces, say f = ( f 1 , f 2 ), on Ω and to boundary forces, say g = (g 1 , g 2 ), along ∂Ω N . Then, due to these forces, the solid is deformed and the corresponding displacement field, say u = (u 1 , u 2 ), is solution of the following boundary value problem

         -div σ(u) = f in Ω , u = 0 on ∂Ω D , σ(u) • n = g on ∂Ω N , ( 1
)
where n is the unit outward normal to ∂Ω and σ(u) is the stress tensor, which is linked to the displacement field thanks to Hooke's law, which reads (with summation convention on repeated indices, those indices i, j, k, l taking on values 1 and 2)

σ i j (u) = R i jkl γ kl (u) ,
where R i jkl is the stiffness tensor and

γ kl (u) = 1 2 (∂ k u l + ∂ l u k ) is the linearized strain tensor.
Now, let us introduce the space of admissible displacements

V = v = (v 1 , v 2 )/v i ∈ H 1 (Ω); v i = 0 on ∂Ω D . ( 2
)
Then u is solution of the following variational formulation

∀v ∈ V , Ω σ i j (u) γ i j (v) = Ω f i v i + ∂Ω N g i v i .
or else, by introducing Hooke's law

         Find u ∈ V such that for all v ∈ V a(u, v) ≡ Ω R i jkl γ kl (u) γ i j (v) = Ω f i v i + ∂Ω N g i v i . (3)
As the stiffness tensor is assumed to satisfy the following symmetry properties: R i jkl = R kli j , the bilinear form a is symmetric :

a(u, v) = a(v, u).
Finally, we recall that, if Ω is a bounded connected open subset and if ∂Ω D has a strictly positive measure, if the stiffness tensor satisfies the two following classical properties (1) R i jkl ∈ L ∞ (Ω) for all indices i, j, k and l, (2) there exists a strictly positive constant C such that, for all symmetric tensor τ: R i jkl τ i j τ kl ≥ C τ i j τ i j , then, there exists a unique u ∈ V solution of the variational problem (3). Moreover, v -→ √ a(v, v) is a norm on V, called the energy norm, which is equivalent to the usual norm ||.|| 1,Ω , and semi-norm |.| 1,Ω thanks to generalized Poincaré inequality, results which are gathered in the following inequalities

α |v| 1,Ω ≤ α ||v|| 1,Ω ≤ a(v, v) ≤ β |v| 1,Ω ≤ β ||v|| 1,Ω , for all v ∈ V , (4)
where 0 < α ≤ β.

Patch method

The "patch" method is now applied to this problem. It means that two meshes are defined: T H is a triangulation of Ω and T h a triangulation of Λ which is a small subdomain of Ω, also assumed to be polygonal (see Fig. 1). Then, we introduce the finite element spaces

           V H = v H ∈ [C(Ω)] 2 / v H|T ∈ (P 1 ) 2 , ∀T ∈ T H ; v H|∂Ω D = 0 V h = v h ∈ [C(Λ)] 2 / v h|t ∈ (P 1 ) 2 , ∀t ∈ T h ; v h|∂Λ∩∂Ω D = 0 M h = V h | Γ
where C(Ω) and C(Λ) are the spaces of continuous functions on Ω and Λ respectively, while P 1 stands for the space of degree 1 polynomial functions. Moreover, Γ = ∂Λ \ ∂Ω is the interface between Λ and Ω \ Λ, and M h is the space for the Lagrange multipliers on Γ. Let us remark that, all along this paper, we shall consider the case of nested triangulations, i.e. we suppose that T h is built up by a subdivision of the elements of T H lying inside Λ. The definitions above are valid for meshes T H and T h composed of triangles. Quadrilateral meshes can be equally well considered by replacing the linear polynomials by the bilinear ones, in the definitions of spaces V H and V h . The patch algorithm will be constructed so that it gives (upon convergence) a solution

u H ∈ V H , u h ∈ V h and λ h ∈ M h to the problem Ω\Λ σ(u H ) : γ(v H ) + Γ λ h • v H = Ω\Λ f • v H , ∀v H ∈ V H , (5)
Λ σ(u h ) : γ(v h ) - Γ λ h • v h = Λ f • v h , ∀v h ∈ V h , (6)
Γ µ h • (u H -u h ) = 0 , ∀µ h ∈ M h . (7)
Here, v • w denotes the usual scalar product in R 2 and we set: σ : γ = σ i j γ i j .

Remark 2. We assume from now on g = 0. This is done only to simplify the expressions. The general case of non zero boundary forces g can be easily covered by adding ∂Ω N \ ∂Λ g • v H and ∂Λ ∩ ∂Ω N g • v h to the right-hand sides of, respectively, (5) and ( 6).

As we are working with nested meshes and given our choice for M h , it is easy to see that the Lagrange multiplier can be completely eliminated. Indeed, (7) is equivalent to u H = u h on Γ. So, if we gather u H and u h into u Hh , such that u Hh = u H on Ω \ Λ and u Hh = u h on Λ, we observe that u Hh belongs to space V Hh defined by

V Hh = v ∈ [C(Ω)] 2 / v = v H on Ω \ Λ , v = v h on Λ with some v H ∈ V H and v h ∈ V h
Adding (5) and ( 6), the problem (5)-( 7) can be rewritten in terms of u Hh alone: Find u Hh ∈ V Hh such that

a(u Hh , v Hh ) = Ω f • v Hh , ∀v Hh ∈ V Hh (8) where a(u, v) = Ω σ(u) : γ(v).
Since our finite element space is conforming, i.e. V Hh ⊂ V (see (2)), under usual properties on meshes (regularity of the meshes, affine family of triangulations), the following a priori error estimate can be derived with Céa's lemma Theorem 3. Problem (8) admits a unique solution u Hh ∈ V Hh . If the exact solution u to (1) belongs to H 2 (Ω), we have moreover the following a priori error estimate

||u Hh -u|| 1,Ω ≤ C H |u| 2,Ω\ Λ + h |u| 2,Λ (9)
with a constant C independent of H, h and u.

We emphasize that equations (5-7) define a coupled problem involving both u H and u h at the same time. For instance, modifying the geometrical definition of Λ will result in a modification of both ( 5) and ( 6). So, when Ω is a domain which is very large and/or finely meshed, the stiffness matrix assembly and factorization will be very CPU time consuming. Moreover, if the patch location Λ has to be modified several times (for instance in shape optimization, crack propagation,...), so will be the stiffness matrix, which is not convenient for numerical efficiency. In the next section, in order to keep unchanged the global numerical operator on Ω, non intrusive coupling is used.

Non-intrusive resolution

Instead of defining equation ( 5) on Ω \ Λ, the idea of non-intrusive coupling is to solve both equations ( 5) and ( 6) on Ω and Λ respectively, in such a manner that u H and u h are at equilibrium on Γ. For this, the term Λ σ(u H ) : γ(v H ) is added on both sides of (5), which allows to use Chasles relation to extend the equation to the whole domain Ω. Then, freshly modified equations (5-6-7) are solved in an iterative manner, starting from a given u 0 H , in accordance with the following relations

Ω σ(u n+1 H ) : γ(v H ) = Ω\Λ f • v H - Γ λ n h • v H + Λ σ(u n H ) : γ(v H ) , ∀v H ∈ V H (10) Λ σ(u n h ) : γ(v h ) - Γ λ n h • v h = Λ f • v h , ∀v h ∈ V h (11) Γ µ h • (u n H -u n h ) = 0 , ∀µ h ∈ M h (12)
It can be convenient to rewrite (10) in an incremental form by adding -Ω σ(u n H ) : γ(v H ) on both sides, which leads to

Ω σ(u n+1 H -u n H ) : γ(v H ) = Ω\Λ f • v H - Γ λ n h • v H - Ω\Λ σ(u n H ) : γ(v H ) , ∀v H ∈ V H or equivalently a(u n+1 H -u n H , v H ) = -R(u n H , v H ) , ∀v H ∈ V H (13) with the following definition of R R(u n H , v H ) = Ω\Λ σ(u n H ) : γ(v H ) - Ω\Λ f • v H + Γ λ n h • v H
In a non-intrusive way, the patch algorithm consists in solving (13) iteratively, which requires of course to solve (11 -12) at each iteration in order to get λ n h . As this problem is solved using numerical methods (in this paper, the finite element method is considered), we introduce K and R, which are the discrete operators corresponding to a and R respectively. Then, the non-intrusive patch algorithm can be seen as a fixed point algorithm to solve equation R(u H ) = 0, according to

u n+1 H = u n H -K -1 R(u n H ) (14)
In the linear elastic structural analysis context, operator K is the stiffness matrix of the global model on Ω, and operator R corresponds to the generalized stress residual between global and local models on Γ. Then, from a practical point of view, algorithm (14) requires only the evaluation of reaction stress on the interface, and implies only linear system resolution at each iteration. Such algorithm is said to be non-intrusive as the global model (stiffness operator K) remains unchanged during all the iterative process. Under certain conditions (the patch model must not be stiffer than the global model on Λ, which will be always verified for mesh refinement in practical applications), this algorithm is shown to converge towards the solution (u H , u h ) of problem (5-6-7), see [START_REF] Chevreuil | A multiscale method with patch for the solution of stochastic partial differential equations with localized uncertainties[END_REF]. More details about this non-intrusive patch algorithm can be found in [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 posteriori estimates for the patch algorithm

As just shown, non-intrusive coupling allows for using a local refined model lying inside a global coarse one. Nevertheless, we still need a decision criterion on the location and size of the patch to define. We recall that our goal is to perform a structural analysis on an adapted mesh so that we can provide a given precision on the whole model. We will use a posteriori error estimation in order to evaluate the solution quality, and decide which parts of the mesh have to be adapted.

The aim of this paper is not to give a detailed review on error estimation. The references are numerous on such a subject. For example, a detailed review of a posteriori error estimation can be found in [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF]Ainsworth & Oden 2000;[START_REF] Grätsch | A posteriori error estimation techniques in practical finite element analysis[END_REF]. So, let us introduce the true error e Hh = u Hhu 1,Ω . Then, there exists three main kinds of global error estimators:

• Recovery based estimators are based upon the comparison between a computed field and a smoothed one [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineering analysis[END_REF], 1992a,b;[START_REF] Carstensen | Averaging technique for FE -a posteriori error control in elasticity. Part I: conforming FEM[END_REF][START_REF] Zang | A new finite element gradient recovery method: superconvergent property[END_REF]. They take advantage of the discontinuous nature of the stress fields arising from the finite element method.

• Residual based estimators are based upon the error representation equation, which can provide upper and lower bounds of the true error. Such estimators can be computed explicitly from the values of the elements internal residual and the normal strain jump across elements boundaries [Babuška & Rheinboldt 1978b;Babuška & Miller 1978;[START_REF] Kelly | A posteriori error analysis and adaptive processes in the finite element method[END_REF] or implicitly from the solution of auxiliary local boundary problems prescribed on elements [START_REF] Demkowicz | Adaptive finite elements for flow problems with moving boundaries. Part I: Variational principles and a posteriori estimates[END_REF]Bank & Weiser 1985] or sub-domains of elements [Babuška & Rheinboldt 1978a;Ainsworth 1993;[START_REF] Díez | Recovering lower bounds of the error by postprocessing implicit residual a posteriori error estimates[END_REF][START_REF] Díez | Accurate upper and lower error bounds by solving flux-free local problems in stars[END_REF][START_REF] Cottereau | Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method[END_REF]].

• Constitutive relation based estimators are built upon the computation of admissible fields (displacement or strain) which satisfies the equilibrium equation div σ + f = 0 [START_REF] Ladevèze | Error Estimate Procedure in the Finite Element Method and Applications[END_REF][START_REF] Ladevèze | Constitutive Relation Error Estimators and Adaptivity in Structural Engineering[END_REF][START_REF] Gallimard | A constitutive relation error estimator based on traction-free recovery of the equilibrated stress[END_REF]]. Those constitutive relation error can also be used to provide bounds of the true error e Hh .

A different approach in error estimates deals with goal oriented error estimators [START_REF] Chamoin | Strict and practical bounds through a non-intrusive and goal-oriented error estimation method for linear viscoelasticity problems[END_REF][START_REF] Ladevèze | Local error estimators for finite element linear analysis[END_REF][START_REF] Rannacher | A posteriori error control in finite element methods via duality techniques: application to perfect plasticity[END_REF][START_REF] Paraschivoiu | A hierarchical duality approach to bounds for the outputs of partial differential equations[END_REF][START_REF] Oden | Goal-oriented error estimation and adaptivity for the finite element method[END_REF]], which provide informations on the quality of the solution regarding local quantities (e.g. a local mean strain) and often need the resolution of a dual problem.

Indeed, there are many more sophisticated error estimators that can provide a guaranteed and precise estimation of the absolute error level, for example [START_REF] Gerasimov | An explicit residual-type error estimator for Q1-quadrilateral extended finite element method in two-dimensional linear elastic fracture mechanics[END_REF][START_REF] Kelly | A posteriori error analysis and adaptive processes in the finite element method[END_REF]. However, in this study, the error estimation is used for two goals:

• the error map is used to define the position and size of the next level of the local patch. For that, only the error distribution is required and a precise estimation of the error level is almost useless;

• at the end of the paper, a rather coarse approximation of the absolute error level is used to stop the iterative local/global coupling solver. At this stage, the error estimator has to be very very cheap as it is evaluated at each iteration. It has to be noted that a relatively large error on this approximation will lead to an over or under-estimation of only few iterations, because the convergence of the iterative solver is linear.

So, finally, we chose here a relatively cheap and easy to implement explicit residual based error estimator. Then, following the non-intrusive coupling strategy, this error estimation should not imply modification of the global existing model, and should not require the resolution of a "large" additional problem. So it will not bound the error as precisely as more accurate estimators (also more computationally expensive) but will allow to build an error map, accurate enough to be used for error driven mesh refinement.

Notations and basic recalls on interpolation

In order to make the paper self consistent, let us begin by introducing the notations and main results on which the remainder of the paper lies.

To derive an a posteriori error estimate, we need an appropriate interpolation operator on space V Hh , that works under certain hypotheses on the mesh, which we properly introduce now.

So, let T be an element of T H or T h . We denote by h T the diameter of T and by ρ T the diameter of the largest ball that can be inscribed into T . Then the mesh is assumed to be regular in the classical sense: it exists a strictly positive constant C such as h T ρ T ≤ C for all T . When element are all triangles or form an affine equivalent family of triangulation, no more hypothesis on the mesh is necessary. Nevertheless, when quadrilateral elements are used, the triangulation is generally not affine equivalent. In that case, the so-called shape regularity assumption (see [START_REF] Matthies | The inf-sup condition for the mapped Q k -P disc k-1 element in arbitrary space dimensions[END_REF]) is done, which imposes that the distortion of the quadrilateral elements from a parallelogram is uniformly bounded, and guarantees that the mapping between the reference element and the original one is bijective.

Now, let us introduce several notations and conventions. Capital letters (H, T , E...) will be used for the global initial "coarse" mesh, namely T H , and lower-case letters (h, t, e...) for the local "refined" mesh, namely T h . Moreover, symbol means that the corresponding inequalities holds up to a multiplicative constant which depends only on the regularity of the mesh.

Then, recall that we use lowest order polynomial functions in the numerical scheme, more precisely, polynomials which total degree in all coordinates is less or equal to 1 for triangles, and polynomials where the maximum power in each coordinate is less or equal to 1 for quadrilaterals.

We shall use the Scott-Zhang interpolation operator [START_REF] Scott | Finite Element Interpolation of Nonsmooth Functions Satisfying Boundary Conditions[END_REF]] on triangular meshes or its extension to quadrilaterals as in [START_REF] Heuveline | H 1 -interpolation on quadrilateral and hexahedral meshes with hanging nodes[END_REF]]. We thus introduce the operator I H : V → V H that, under the abovementioned regularity assumptions, satisfies for any v ∈ V

||v -I H v|| 0,T H T |v| 1, ω T , |v -I H v| 1,T |v| 1, ω T , ||v -I H v|| 0,E H E |v| 1, ω E (15)
on every element T ∈ T H and on every edge E of mesh T H , while H T is the diameter of element T and H E those of edge E. Finally, ω T (resp. ω E ) denotes the union of neighbours in T H of element T (resp. edge E). If we follow the construction of [START_REF] Scott | Finite Element Interpolation of Nonsmooth Functions Satisfying Boundary Conditions[END_REF], it can be seen that there is a certain freedom in the actual construction of I H and in the choice of ω T and ω E . In particular, for any node x of T H lying on the interface Γ, I H v(x) can be constructed by averaging v over an edge of T H lying on Γ. Accordingly, we can assume that if T ⊂ Ω \ Λ, then ω T ⊂ Ω \ Λ, and the same for ω E for any edge E inside Ω \ Λ. It means that when T (resp. E) touches the interface Γ, we can choose ω T (resp. ω E ) so that it only contains neighbours of T (resp. E) that are outside Λ. Otherwise, ω T is defined as the union of all the neighbouring elements of T . We recall also that I H is a projector, i.e.

I H v H = v H for any v H ∈ V H .
Similarly, it exists an interpolation operator

I h : V(Λ) → V h , where V(Λ) = {v ∈ H 1 (Λ) / v| ∂Ω D ∩∂Λ = 0}. And we have for any v ∈ V(Λ) ||v -I h v|| 0,t h t |v| 1, ω t , |v -I h v| 1,t |v| 1, ω t , ||v -I h v|| 0,e h e |v| 1, ω e (16)
on every element t ∈ T h and on every edge e of mesh T h , while h t is the diameter of element t and h e those of edge e. Again, there is a certain freedom in the construction of I h . In particular, it can be chosen so that (I h v)| Γ only depends on the trace of v on Γ, so that for any boundary edge e of T h lying on Γ

||v -I h v|| 0,e ||v|| 0, γ e
where γ e is the union of boundary edges of T h on Γ that share an endpoint with e.

Upper bound for the error

Using notations introduced in the previous section, we define now the combined interpolation operator I Hh : V → V Hh as follows. For all v ∈ V, I Hh v is the element of V Hh that takes the following values at the nodes of the two meshes

I Hh v(x) = I H v(x) if x is a node of mesh T H outside Λ or on Γ I h v(x) if x is a node of mesh T h inside Λ or on ∂Λ ∩ ∂Ω
It can be observed that relations ( 15) and ( 16) give error estimates for functions which support is contained respectively in Ω \ Λ and Λ. However, the definition of the interpolation operator I Hh makes it necessary to study more precisely what happens for a function whose support has an intersection with Γ. It is the aim of the following proposition.

Proposition 4. Let v ∈ H 1 (Ω) and t be an element of mesh T h that touches Γ, which means that at least one vertex of t is on Γ. Then, we have ||v -

I Hh v|| 0,t h t |v| 1, ω t + h 1/2 t ||v -I H v|| 0, ∂ω t ∩ Γ (17)
where ω t is the union of all elements of T h that share at least one node with t.

Similarly, let e be an internal edge of mesh T h that touches Γ, which means that one endpoint of e is on Γ. Then, we have ||v -

I Hh v|| 0,e h 1/2 e |v| 1, ω e + ||v -I H v|| 0, ∂ω e ∩ Γ ( 18
)
where ω e is the union of all elements of T h that share at least one node with e.

Proof. Let us denote v h = I Hh v| Λ .

• Let us begin with inequality (17). When an element t touches Γ, two cases may occur:

1. t has a side, say ẽ, that lies on Γ. So all the nodes of t, which are not on ẽ, are internal nodes of mesh T h and, consequently, v h -I h v vanishes at all those nodes. Then, on t, v h -I h v is a first order polynomial function (or a bilinear polynomial) which is fully determined by its values at the nodes on ẽ. So a simple argument involving scaling and equivalence of norms on a finite dimensional space allows us to write

||v h -I h v|| 0,t h 1/2 t ||v h -I h v|| 0,ẽ .
Then, using the fact that I h is a projector on V h and the interpolation estimate on the boundary edges, we obtain

||v h -I h v|| 0,t h 1/2 t ||I h (v h -v)|| 0,ẽ h 1/2 t ||v h -v|| 0, γ ẽ .
2. t has only a vertex, say x, that lies on Γ. Since v h -I h v vanishes at all the vertices of t except x, we get again by scaling

||v h -I h v|| 0,t h t |(v h -I h v)(x)| .
Let ẽ be an edge of T h which contains x and lies on Γ. For all w h ∈ V h , the inverse inequality leads to

|w h (x)| h -1/2 t ||w h || 0,ẽ .
Taking now w h = v h -I h v and following as in the previous case, we obtain

||v h -I h v|| 0,t h t |(v h -I h v)(x)| h 1/2 t ||v h -I h v|| 0,ẽ = h 1/2 t ||I h (v h -v)|| 0,ẽ h 1/2 t ||v h -v|| 0, γ ẽ .
Finally, as γ ẽ ⊂ (∂ω t ∩ Γ), in both cases, we have

||v h -I h v|| 0,t h 1/2 t ||v h -v|| 0, ∂ω t ∩Γ
and we conclude by using triangular inequality and above mentioned interpolation estimate

||v -v h || 0,t ≤ ||v -I h v|| 0,t + ||v h -I h v|| 0,t h t |v| 1, ω t + h 1/2 t ||v h -v|| 0, ∂ω t ∩Γ ,
which is equivalent to (17) since v h = I Hh v on Λ and v h = I H v on Γ.

• Let us now turn to the proof of ( 18). Let e be an internal edge of mesh T h that touches Γ. It means that one endpoint of e, say x, lies on Γ while the other is inside Λ. Then, we get again by scaling

||v h -I h v|| 0,e h 1/2 e |(v h -I h v)(x)| .
Using inverse inequality as in the previous part of the proof, we conclude

||v h -I h v|| 0,e ||v h -I h v|| 0,ẽ = ||I h (v h -v)|| 0,ẽ ||v h -v|| 0, γ ẽ ,
where ẽ is an edge of T h which contains x and lies on Γ. Finally, triangular inequality gives

||v -v h || 0,e ≤ ||v -I h v|| 0,e + ||v h -I h v|| 0,e h 1/2 e |v| 1, ω e + ||v h -v|| 0, γ ẽ , which implies (18) since v h = I Hh v on Λ and v h = I H v on Γ.
Before giving the main result of this section, we need some additional notations and definitions.

• Let us recall that two meshes are defined: T H on Ω and T h on Λ, which is a subdomain of Ω. As we shall need a partition of the whole domain Ω, it is convenient to introduce T H (Ω \ Λ) that is the union of elements from T H that are outside Λ.

• As the edges of elements belonging to T H (Ω \ Λ) and T h will be used, we split the corresponding set into three disjoint ones: E H (Ω \ Λ) is the union of edges of T H lying outside Λ (which excludes those on interface Γ), E h the union of edges of T h lying inside Λ except those on Γ, and E H (Γ) the union of edges of T H lying on Γ.

• Moreover, we introduce the two subsets T b h ⊂ T h and E b h ⊂ E h , which are respectively the sets of elements and internal edges of the fine mesh that touch Γ.

• As integrations by parts will be used in the following, we have to precise what we call jumps along the edges.

We thus assign a unit normal vector n ε to every edge ε of both meshes T H and T h . Then, for any quantity w defined on both sides of ε, we set [w] = w -w + , with w ± (x) = lim t→0 + w(x ± tn ε ), for any x on ε. It means in particular that the jump [w • n ε ] of a vector field w is defined independently from the arbitrary choice of the direction of n ε . If an edge ε lies on ∂Ω, we set [w] = w when ε is contained in ∂Ω N and [w] = 0 when ε is contained in the Dirichlet boundary ∂Ω D .

Then, the following error estimate holds.

Proposition 5. Under the regularity assumptions on the mesh, recalled at the beginning of this section, let u be the solution to (1) and u Hh the solution to (8). We have

||u -u Hh || E Θ(u Hh ) α , with Θ(u Hh ) 2 = T ∈T H (Ω\Λ) H 2 T || f + div σ(u H )|| 2 0,T + E∈E H (Ω\Λ) H E || [σ(u H ) n E ] || 2 0,E + E∈E H (Γ) H E || [σ(u Hh ) n E ] || 2 0,E + t∈T h h 2 t || f + div σ(u h )|| 2 0,t + e∈E h h e || [σ(u h ) n e ] || 2 0,e + t∈T b h H E,t h t || f + div σ(u h )|| 2 0,t + e∈E b h H E,e || [σ(u h ) n e ] || 2 0,e (19)
where

|| • || E = √ a(•, •)
is the energy norm while H E,t and H E,e stand for the maximum size of the elements of T H to which element t and edge e are attached. We finally recall that symbol means that the corresponding inequality holds up to a multiplicative constant which only depends on the mesh regularity.

Proof. Let us set w = uu Hh . Then, by Galerkin orthogonality, we have, for any w Hh ∈ V Hh

||u -u Hh || 2 E = a(u -u Hh , w) = a(u -u Hh , w -w Hh ) = Ω σ(u -u Hh ) : γ(w -w Hh ) = T ∈T H (Ω\Λ) T σ(u -u H ) : γ(w -w H ) + t∈T h t σ(u -u h ) : γ(w -w h )
An integration by parts over all elements yields to

||u -u Hh || 2 E = T ∈T H (Ω\Λ) T ( f + div σ(u H )) • (w -w H ) - E∈E H (Ω\ Λ) E [σ(u H ) n E ] • (w -w H ) - E∈E H (Γ) E [σ(u Hh ) n E ] • (w -w H ) + t∈T h t ( f + div σ(u h )) • (w -w h ) - e∈E h e [σ(u h ) n e ] • (w -w h )
So, taking w Hh = I Hh w and using interpolation estimates ( 15), ( 16), ( 17) and ( 18) for the appropriate elements/edges, we obtain Let us recall that T b h is the set of elements of the fine mesh T h that touch Γ. We shall decompose it into a collection of disjoint sets T b h (E) with E running over E H (Γ). Each set T b h (E) will contain elements of T h that have at least one vertex on E. The attribution of elements of T b h to the edges from E H (Γ) can be done as follows: we remind first of all that each t ∈ T b h has an edge or a vertex on Γ. If it is an edge, say e, as we suppose that T h is built up by a subdivision of the elements of T H lying inside Λ, it exists a unique edge E 0 ∈ E H (Γ) such as e ⊂ E 0 and we put t in T b h (E 0 ). If t has a vertex, say x, on Γ, two cases may occur. First, it exists a unique edge E 0 ∈ E H (Γ) such as x ∈ E 0 and we put t in T b h (E 0 ). Second, x is an endpoint of two edges of E H (Γ): {x} = E 1 ∩ E 2 ; then, we choose the longest edge, say E i , and put t in T b h (E i ) (if the two edges have the same length, we choose the first one). Then, using again (15) leads to

||u -u Hh || 2 E T ∈T H (Ω\Λ) || f + div σ(u H )|| 0,T H T |w| 1 , ω T + E∈E H (Ω\ Λ) || [σ(u H ) n E ] || 0,E H E |w| 1 , ω E + E∈E H (Γ) || [σ(u H ) n E ] || 0,E H E |w| 1 , ω E + t∈T h || f + div σ(u h )|| 0,t h t |w| 1 , ω t + t∈T b h || f + div σ(u h )|| 0,
t∈T b h || f + div σ(u h )|| 0,t h t ||w -I H w|| 0 , ∂ω t ∩Γ = E∈E H (Γ) t∈T b h (E) || f + div σ(u h )|| 0,t h t ||w -I H w|| 0 , ∂ω t ∩Γ E∈E H (Γ)           t∈T b h (E) h t || f + div σ(u h )|| 2 0,t           1/2 ||w -I H w|| 0,E E∈E H (Γ)           t∈T b h (E) H E h t || f + div σ(u h )|| 2 0,t           1/2 |w| 1 , ω E           t∈T b h H E,t h t || f + div σ(u h )|| 2 0,t           1/2         E∈E H (Γ) |w| 2 1 , ω E         1/2
Similarly, we may write

E b h = ∪ E∈E H (Γ) E b h (E)
, where E b h (E) are disjoint sets of edges of E h that have at least one vertex on E. Following the same idea as above, with (15) again, one has

e∈E b h || [σ(u h ) n e ] || 0,e ||w -I H w|| 0 , ∂ω e ∩Γ           e∈E b h H E,e [σ(u h ) n e ] || 2 0,e           1/2         E∈E H (Γ) |w| 2 1 , ω E         1/2
Getting all this together leads to

||u -u Hh || 2 E         T ∈T H (Ω\Λ) H 2 T || f + div σ(u H )|| 2 0,T + E∈E H (Ω\ Λ) H E || [σ(u H ) n E ] || 2 0,E + E∈E H (Γ) H E || [σ(u Hh ) n E ] || 2 0,E + t∈T h h 2 t || f + div σ(u h )|| 2 0,t + e∈E h h e || [σ(u h ) n e ] || 2 0,e + t∈T b h H E,t h t || f + div σ(u h )|| 2 0,t + e∈E b h H E,e || [σ(u h ) n e ] || 2 0,e           1/2 |w| 1 , Ω
and the result since

| • | 1 , Ω is a norm equivalent to || • || E (see(4)).

Application to adaptive mesh refinement

As seen in previous sections, the main goal of this work is to provide an non-intrusive mesh refinement procedure. By the way, even if very efficient remeshing algorithms have been developed [Alliez et al. 2003;[START_REF] Bois | A fully optimal anisotropic mesh adaptation method based on a hierarchical error estimator[END_REF]], such a procedure remains complex. Moreover, when dealing with commercial software, local mesh adaptation for quadrangular elements is an even more complex task [START_REF] Sandhu | Hierarchical mesh adaptation of 2D quadrilateral elements[END_REF][START_REF] Woodbury | Localized coarsening of conforming all-hexahedral meshes[END_REF], which is not as efficient as triangular mesh adaptation. A way to avoid main remeshing difficulties is to rely on non-conforming finite elements, allowing for hanging nodes. Then it is possible to locally refine a mesh without taking care of the conformity of elements, at the cost of additional complexity, for example the use of a mortar method [Bernardi & Maday 2000]. In this section, we will describe how the patch algorithm can be used for local remeshing, in a non-intrusive way.

So, it is considered as a test case an infinite plate with a hole (radius a), under tensile loading (magnitude T ∞ ). For numerical simulations, we only study a finite quarter of the plate (see Fig. 2) with symmetrical boundary conditions and prescribed displacements, given by the analytical solution of the infinite plate problem, which reads

               u x (r, θ) = T ∞ a 8µ r a (κ + 1) cos θ + 2a r (1 + κ) cos θ + cos 3θ - 2a 3 r 3 cos 3θ u y (r, θ) = T ∞ a 8µ r a (κ -3) sin θ + 2a r (1 -κ) sin θ + sin 3θ - 2a 3 r 3 sin 3θ with κ = 3 -ν 1 + ν and µ = E 2(1 + ν)
, E being Young's modulus and ν Poisson's ratio.

A local error indicator

To be able to use the previously introduced error estimate for mesh refinement, we have to build a local error indicator, i.e. an error map providing a scalar value on each element of the finite elements mesh. In order to remain consistent, the quadratic sum of all elementary indicators θ T must be equal to the error estimate Θ(u Hh ), given in (19).

Θ(u

Hh ) =         τ∈T H (Ω\Λ) ∪ T h θ 2 τ         1/2 .
So, the elementary indicators will be set up following the above definitions, for mesh T H (Ω \ Λ)

θ T =        H 2 T || f + div σ(u H )|| 2 0,T + 1 2 E∈∂T H E || [σ(u Hh ) ] || 2 0,E        1/2 if T ∈ T H (Ω \ Λ) ( 20
)
where E is a generic edge of the boundary ∂T of element T . For mesh T h , we set

                                   θ t =        h 2 t || f + div σ(u h )|| 2 0,t + 1 2 e∈∂t h e || [σ(u h ) n e ] || 2 0,e        1/2 if t ∈ T h \ T b h θ t =        h 2 t || f + div σ(u h )|| 2 0,t + 1 2 e∈∂t h e || [σ(u Hh ) n e ] || 2 0,e + H E,t h t || f + div σ(u h )|| 2 0,t + 1 2 e∈E b h ∩∂t H E,e || [σ(u h ) n e ] || 2 0,e           1/2 if t ∈ T b h ( 21
)
where e is a generic edge of the boundary ∂t of element t.

Remark 6. A coefficient 1 2 appears in front of the edge terms because most of the edges belong to two elements. However, we recall that, if an edge is contained in ∂Ω N , the normal jump has no meaning and has to be taken between the field defined in Ω (or Λ) and the Neumann boundary condition. In such a case, coefficient 1 2 has to be omitted. Furthermore, when edge is contained in the Dirichlet boundary, there is no normal jump and the corresponding term has to be ignored.

Local error visualization

Let us emphasize that, as far as mesh refinement is concerned, the local error indicators given by (20-21) will be used. Nevertheless, when error map visualization is concerned, those indicators appear to be not well suited, because the weight due to the element size on the local error indicator leads to very small values for refined elements, compared to the initial ones. So, we propose to use the following relative error indicator in visualized error maps

θ T = θ T √ meas(T ) (22)
where meas(T ) is the measure (here the surface) of element T . We then get a relative error estimator which will appear to be more convenient for such a multi-scale problem.

Another possibility to take into account different scales is to collect together local error indicators to the global mesh, according to

             ∀T ∈ T H (Ω \ Λ) , θ T = θ T , ∀T ∈ T H (Λ) , θ T =         t ∈T h , t ⊂T θ 2 t         1/2 . ( 23
)
Eventually, indicator θ T is designed to assess the error on each refined mesh (and then the quality of the mesh), whereas indicator θ T is better to provide comparison over a unique initial mesh. Now, in order to evaluate our residual based error indicator, we need to make comparison with the true error |u-u Hh | E,T . Then, it is possible to compute the corresponding absolute and relative error indicators θ re f T and

θ re f T            θ re f T = |u -u Hh | E,T , θ re f T = θ re f T √ meas(T )
However, in order to give a relevant comparison between the two error maps, corresponding to θ T and θ re f T , we scale them with their infinite norm

Θ ∞ = max T ∈T H ∪T h θ T , Θ re f ∞ = max T ∈T H ∪T h θ T re f .
So we introduce a normalized error ratio, which is the ratio between the scaled residual based error and the scaled reference error:

θ T θ re f T × Θ re f ∞ Θ ∞ .
Remark 7. As we are considering triangulations of a polygonal domain, we naturally induce an error on the geometry approximation at the hole edge: It will be a part of the discretization error.

Mesh refinement with patches of finite elements

Starting from an initial mesh on Ω, the classical residual based error estimate, defined by ( 20), is used to build an error map, which is then used to define the patch location [START_REF] Rezzonico | Multiscale algorithm with patches of finite elements[END_REF]]. A tolerance being given, the patch Λ will be the union of elements of T H which local error indicator θ T is greater than .

Λ =

T ∈T H , θ T > T .

The patch mesh is then uniformly refined by elements subdivision (see Fig. 5).

OMEGA_S OMEGA_P

Figure 5: Uniformly refined patch

In fact, we only have to take care of the interface, in order to guarantee that the meshes are nested on Γ, i.e. V H | Γ ⊂ V h | Γ . Let us remark that using such a patch for mesh refinement also allows us to locally modify the model geometry, for instance by moving some nodes in order to better fit the hole geometry (see Fig. 5). Then, algorithm ( 13) is applied, with finite element meshes T H on Ω and T h on Λ. It leads to a multi-scale solution (u H , u h ) of problem (5-6-7), given Fig. 6. So a new residual based error estimate can be computed using (20-21). For example, Fig. 7 and Fig. 8 provide illustrations of uniformly refined patches for several mesh size ratio on the coupling interface Γ. The patch is uniformly refined four times (each time a triangle a split into four smaller homothetic ones). Let U be the exact solution interpolated at the mesh nodes, and R(u n H ) be the residual at iteration n, see ( 14). Then, the relative residual at iteration n is defined as

n = ||R(u n H )|| 2 ||KU|| 2 .
Fig. 9 gives those residual evolution over iterations of the algorithm. Depending on the patch coarseness, about 30 iterations are sufficient to reach computer precision (in fact, only a few iterations are required to reach a reasonable precision, as shown on the figure and developed in the next section). Moreover, the more the patch is coarse, the less iterations are needed: Such phenomenon is due to the stiffness gap between the global model and the patch, see [START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF]].

We previously presented examples of uniformly refined patches for simplicity, but nothing prevents the local patch from being non-uniformly refined, as illustrated on Fig. 10, for which a classical mesh refinement procedure is used but only in the patch. By classical mesh refinement procedure, we mean that some triangles are divided according to a certain criterion and the mesh is made conformal after some mesh refinement iterations, the procedure being the one which exists in Code Aster [ Électricité de France 2014], the software we used in our numerical tests. • A good correspondence between the relative non-dimensional errors (estimated and reference). In fact, as we aim at mesh refinement, we are more interested in the local error map than in the global absolute error.

• A moderate dispersion of the normalized error ratio across the finite element mesh. Note that the factor ten which is observed on the scale is mainly due to corner elements. As the displacement is fully enforced by Dirichlet condition, e.g. at the top right corner, the error overestimation is not relevant.

• A simple, effective and cheap tool to drive local mesh refinement in a non-intrusive way.

Multi-patch and h-p refinement

Mesh refinement based on such algorithm allows a wide range of non-conforming couplings between the global model and the patch. In the following example (Fig. 11), a mortar method is used for coupling different element sizes and degrees: the global model uses linear quadrilaterals whereas the local patch is made of quadratic elements of smaller size. In this framework, non-conforming refinements using mortar method can easily be applied in order to locally refine quadrangular meshes, which is much more difficult when one wants to keep conforming meshes.

Finally, it would also be possible to mix different element types, for instance a patch with triangular mesh within a quadrangular global one, without any difficulty. Let us also note that remeshing techniques with patches may also be set up in a more general context than classical finite element method, e.g. in the frame of isogeometric analysis [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Chemin | Isogeometric local h-refinement strategy based on multigrids[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 the previous sections, we focused on adaptive mesh refinement. Each example has been computed using the non-intrusive iterative algorithm, and the presented results arose out of the converged solution of problem ( 10)-( 12), 20 namely when the numerical computer precision was reached. Of course, the quality of the computed solution may become acceptable much earlier. In fact, it is worthless to carry on additional iterations once the iterative process has reached a point such that the coupling error is smaller than the finite element one. The aim of this section is to investigate this point.

Convergence error versus finite element error

As the numerical solution to ( 8) is obtained through the iterative algorithm ( 10)-( 12), at iteration number n, we have an approximation u n Hh ∈ V Hh , such that u n Hh = u n H outside Λ and u n Hh = u n h inside Λ, u n h being such that

Λ σ(u n h ) : γ(v h ) = Λ f • v h + Γ λ n h • v h , ∀v h ∈ V h .
Now, let us introduce space M H as the trace space of V H on Γ. Then, we define λ n H ∈ M H such as

Ω\Λ σ(u n H ) : γ(v H ) - Ω\Λ f • v H = Γ λ n H • v H , ∀v H ∈ V H .
Here Γ λ n H • v H represents the generalized reaction stress of domain Ω\Λ on interface Γ. Then, if we set λ n Hh = λ n H + λ n h , adding the two previous equations yields to

a(u n Hh , v Hh ) = Ω f • v Hh + Γ λ n Hh • v Hh , ∀v Hh ∈ V Hh (24)
And we have the following result.

Proposition 8. Let u be the solution to (1) and u n Hh be constructed on iteration number n of the algorithm ( 10)-( 12). Under the hypotheses of Proposition 5, we have

||u -u n Hh || E Θ(u n Hh ) + Θ n α ( 25
)
where Θ(u n Hh ) is given by ( 19) but calculated with u n Hh instead of u Hh and

Θ n = L Ω,Γ         x∈N H (Γ) 1 H x Γ λ n Hh φ x 2         1 2 (26)
where N H (Γ) is the set of nodes of the coarse mesh located on interface Γ. Moreover, for any x ∈ N H (Γ), φ x is the shape function in V H at node x and H x is the mesh size at those node x, i.e. the length of an edge of T H adjacent to x. Finally, L Ω,Γ is a constant depending on Ω and Γ.

Proof. Set w = uu n Hh . Then, with (24), we have for any w Hh ∈ V Hh

||u -u n Hh || 2 E = a(u -u n Hh , w) = a(u -u n Hh , w -w Hh ) - Γ λ n Hh • w Hh = Ω σ(u -u n Hh ) : γ(w -w Hh ) - Γ λ n Hh • w Hh
Now, we take w Hh = I Hh w. Following the proof of Proposition 5, it is easy to bound the first term by Θ(u n Hh ) |w| 1,Ω , Θ(u n Hh ) being given by ( 19).

For the second term, we first observe that:

w Hh = w H = x∈N H (Γ)
w H (x) φ x on Γ. Then, using the properties of Scott-Zhang interpolation operator, it is possible to set w H (x) = I H w(x), for any x ∈ N H (Γ), in such a way it is calculated as a properly weighted average of w over edges adjacent to x and lying on Γ, and we have by inverse inequality. Hence (27) shows L Ω,Γ should be homogeneous to length. It can be chosen in practice as some combination, having the right dimension, of characteristic lengths, such as the diameter of Ω and the length of Γ. In order to get more insight into a plausible value of this constant, let us consider the following simple geometry: Ω is a disk of radius R centered at the origin and Γ is a circle of radius ρ < R centered at the same point. Assume also ∂Ω D = ∂Ω.

Γ λ n Hh • w Hh = x∈N H (Γ) w H (x) Γ λ n Hh • φ x ≤         x∈N H (Γ) 1 H x

Inequality

Taking any function w ∈ H 1 (Ω) vanishing on ∂Ω and introducing the polar coordinates (r, θ), we can estimate

w 2 0,Γ = 2π 0 w(ρ, θ) 2 ρdθ = 2π 0 R ρ ∂w ∂r (r, θ)dr 2 ρdθ ≤ 2π 0 R ρ ∂w ∂r (r, θ) 2 rdr R ρ 1 r dr ρdθ = ρ ln R ρ 2π 0 R ρ ∂w ∂r (r, θ) 2 rdrdθ = ρ ln R ρ |w| 2 1,Ω
where we have used the fundamental theorem of calculus, having in mind that w(R, θ) = 0, and Cauchy-Schwarz inequality. This suggests that one should take L Ω,Γ = ρ ln R ρ in such a geometry.

Remark 9. The previous value is optimal in this case, as seen on the function w defined on Ω by w(r, θ) = ln R r for ρ ≤ r ≤ R and w(r, θ) = ln R ρ for 0 ≤ r ≤ ρ.

Numerical test

To illustrate this result, we consider a test case similar to the previous ones. However, it is changed in the following way. In addition to mesh refinement, the local patch also modifies the geometry of domain Λ, see Fig. 12. Moreover, within the global model, the hole is replaced by a rigid inclusion, which Young's modulus is a hundred time higher than in the rest of Ω. Those hole is only taken into account in the local patch. At the end, such a situation will not differ from the previous one (see Fig. 1), as the converged solution of algorithm ( 10)-( 12) satisfies (5-7). As equation ( 5) is defined on Ω\Λ, the global solution u H on Λ is only a fictitious prolongation, which has no physical meaning and which values upon convergence only depend on the algorithm initialization.

Filling the hole with a rigid inclusion leads to a greater number of iterations. Such a choice is done in order to get a more readable situation. Indeed, in the previous one, the convergence is too fast and the number of iterations is not sufficient to be able to analyse the convergence error against the finite element error.

Figure 12: Uniformly refined path with geometric changes As far as coefficient L Ω,Γ , which appears in estimator Θ n , see ( 26), is concerned, we proceed as follows. Let us recall that, in the particular case where Ω is a disk of radius R centered at the origin and Γ a circle of radius ρ < R centered at the same point, we obtained L Ω,Γ = ρ ln R ρ . Introducing L Γ the length of interface Γ and L Ω the diameter of Ω, in the above mentioned particular case, we have: L Γ = 2 π ρ and L Ω = 2 R, which leads to :

L Ω,Γ = L Γ 2 π ln π L Ω L Γ . (28)
This is the value we take in our numerical test.

We focus now on convergence properties. Iterative process ( 14) can be seen as a modified Newton's method (also called Chord method), see [START_REF] Kelley | Solving Nonlinear Equations with Newton's Method[END_REF]]. So the convergence error estimator Θ n should satisfy Θ n+1 ≤ c Θ n , with c ∈ [0, 1]. Thus, it provides linear convergence as one can expect from this iterative method (see Fig. 13, where we observe numerically that c ≈ 0.685). Consequently, a smart way to decide whether iterations may be stopped, is to compare the convergence error Θ n with the finite element error Θ(u n Hh ). Indeed, at a time, we shall have Θ n Θ(u n Hh ). For instance, one can decide that the computed solution u n Hh at iteration n is acceptable as soon as the convergence error represents a small percentage of the total error (see Fig. 14). For example, in our case, we see that after only 10 iterations, the convergence error Θ n is less than one percent of the total error Θ(u n Hh ) + Θ n .

Last but not least, Fig. 15 gives a comparison between the estimated error Θ(u n Hh ) + Θ n and the reference error Θ re f . It can be seen that relation ( 25) is satisfied and provides good accuracy over iterations. Moreover, with the choice of L Ω,Γ given by ( 28), the ratio between Θ(u n Hh ) + Θ n and Θ re f remains between 2 and 4, which is rather nice in view of the way it was obtained. In order to better appreciate the errors comparison, note that ||u|| E ≈ 28 kJ. Remark 10. The reader may wonder about the accuracy of L Ω,Γ in a general case. In fact, as stated previously, the convergence of global/local coupling algorithm is linear: that means a substantial error on this constant implies few additional iterations only. So having a coarse approximation of the order of magnitude of this constant is sufficient. Moreover, and it is more important, to be relevant, it has to be very cheap compared to the computational cost of one iteration. Besides, even with this coarse approximation of L Ω,Γ , Fig. 15 exhibits a good agreement between the error estimated using this approximation of the constant and the exact error.

Finally, we show the error map evolution for the 9 first iterations of the non-intrusive algorithm. At the beginning, the error is concentrated around the interface, due to the lack of equilibrium on Γ. But after only few iterations, error near the interface becomes insignificant compared to error on the rest of the domain (essentially at hole edge). Such an observation is in accordance with the comparison provided between Θ(u n Hh ) and Θ n . In fact, less than 10 iterations are needed to ensure that the coupling error is lower than the finite element discretization error. Practically, in such a case, it means that the non-intrusive algorithm should be stopped after about 10 iterations and that a new refinement iteration might be done, for example if a given convergence criteria on the expected error is not achieved.

It is clear that computing the a posteriori error estimation at each iteration of the patch algorithm could appear time consuming, compared to the cost of a single iteration. Nevertheless, for some applications, saving a few iterations can be wholesome, even if it requires the computation of Θ(u n Hh ) several times, for example in the case of nonlinear problems using a rather crude error estimator.

Conclusion

This paper presents the construction of an ad hoc residual based explicit error estimator, for the assessment of multiscale solutions, associated with a non-intrusive mesh refinement technique. The extension of such residual error estimators would be straightforward for higher degree polynomial functions and affine family of triangulations. Moreover, the examples show that the proposed method is efficient with the given error estimate. However, it should be interesting to investigate other approaches, with more refined weightings such as [START_REF] Carstensen | Constants in Clément-interpolation error and residual based a posteriori error estimates in finite element methods[END_REF], in order to try to improve even further the results. This method helps to overcome the inherent cost of classical mesh adaptation, and make it easier for some complex situations (local refinement of quadrilateral or hexahedral meshes, non-conforming meshes, local h-p refinement), as illustrated on various examples. Moreover, in order to reduce the global cost of the process, a pragmatic criterion is introduced, based on a posteriori estimate of the convergence error of the non-intrusive algorithm, which allows to reduce the number of iterations. This first study calls for further analyses including for example local or global nonlinearities [START_REF] Duval | Non-intrusive Coupling: Recent Advances and Scalable Nonlinear Domain Decomposition[END_REF], or shell/solid coupling [START_REF] Guguin | Non intrusive coupling between 3D and 2D laminated composite models based on finite element 3D recovery[END_REF], non-matching meshes [START_REF] Bouclier | Local enrichment of nurbs patches using a non-intrusive coupling strategy: Geometric details, local refinement, inclusion, fracture[END_REF], model error... It is why such an approach, which enables to use any error estimator and, in the case of complex and large structures, several patches, is particularly attractive. Actually, it would make possible an easy use of different strategies of mesh refinement associated with different error estimators (for example goal oriented error estimators) for each patch, depending on the nature of the problem.

Figure 1 :

 1 Figure 1: Patch location and meshes

Figure 3 :

 3 Figure 2: Test case depiction

Figure 4 :

 4 Figure 4: Relative residual error on the initial mesh

Figure 7 :Figure 8 :Figure 9 :

 789 Figure 6: Test case solution with a uniformly refined patch

Figure 10 :

 10 Figure 10: Relative residual errors on non-uniformly refined patches: residual based error θ T -reference error θ re f Tnormalized error ratio

Figure 11 :

 11 Figure 11: Multi-patch h-p refinement

Figure 13 :

 13 Figure 13: Evolution of the convergence error Θ n

Figure 14 :

 14 Figure 14: Evolution of the convergence error Θ n and the finite elements error Θ(u n Hh )

Figure 15 :

 15 Figure 15: Evolution of the estimated total error Θ(u n Hh ) + Θ n and the reference error Θ re f

Figure 16 :

 16 Figure 16: Scaled error map evolution: θ T (u n Hh )

|w H (x)| ||w|| 0,E √ H x

Ainsworth M. (1993) A unified approach to a posteriori error estimation using element residual methods. Numerische Mathematik 65(1):23-50 Ainsworth M. and Oden J.T. (2000) A Posteriori Error Estimation in Finite Element Analysis. John Wiley & Sons

Acknowledgements

The authors would like to acknowledge the financial support of the Agence Nationale de la Recherche under Grant ICARE ANR-12-MONU-0002.