Raju Basava
email: basava-raju.akula@mines-paristech.fr

Akula

Julien Vignollet
email: julien.vignollet@safrangroup.com

Vladislav A Yastrebov
email: vladislav.yastrebov@mines-paristech.fr

Stabilized MorteX method for mesh tying along embedded interfaces

Keywords: mesh tying, embedded interface, MorteX method, stabilization, mortar method, X-FEM

We present a unified framework to tie overlapping meshes in solid mechanics applications. This framework is a combination of the X-FEM method and the mortar method, which uses Lagrange multipliers to fulfill the tying constraints. As known, mixed formulations are prone to mesh locking which manifests itself by the emergence of spurious oscillations in the vicinity of the tying interface. To overcome this inherent difficulty, we suggest a new coarse-grained interpolation of Lagrange multipliers. This technique consists in selective assignment of Lagrange multipliers on nodes of the mortar side and in non-local interpolation of the associated traction field. The optimal choice of the coarse-graining spacing is guided solely by the mesh-density contrast between the mesh of the mortar side and the number of blending elements of the host mesh. The method is tested on two patch tests (compression and bending) for different interpolations and element types as well as for different material and mesh contrasts. The optimal mesh convergence and removal of spurious oscillations is also demonstrated on the Eshelby inclusion problem for high contrasts of inclusion/matrix materials. Few additional examples confirm the performance of the elaborated framework.

Introduction

The finite element method (FEM) is used to solve a wide range of physical and engineering problems. Based on a variational formulation and a discretized representation of the geometry, this method is extremely flexible in handling complex geometries, non-linear and heterogeneous constitutive equations and multi-physical/multi-field problems. A classification of finite element models can be proposed based on the strategy to represent the boundary of the computational mesh. Classical FE meshes fall into the category of "boundary fitted" (BF) methods, where the boundaries of the physical and computational domains coincide [Fig. 1(a)]. Alternatively, for "embedded/immersed boundary" (EB) methods, the computational domain is a mesh or a Cartesian grid hosting another physical domain [Fig. 1(b)]. Note that material properties or even the governing equations of the host medium and the embedded one can be different. Within the EB method, the geometry contour can be embedded either fully or partially. The BF methods [see Fig. 1(a)] can be used to solve boundary value problems where the boundary conditions are prescribed on surfaces explicitly represented by a mesh. The EB methods [see Fig. 1(b)] can handle a broader class of problems. First, it can be used to solve the same boundary value problems as BF methods but with the boundary represented by a rather general level-set function without explicit discretization of the physical domain. Second, the embedded boundary can serve as a material interface to include features such as inclusions, voids or even cracks for linear/quadratic background meshes [START_REF] Fries | Higher-order conformal decomposition FEM (CDFEM)[END_REF]. In addition, within the so-called CutFEM [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF][START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF], debonding of the interface can be taken into account as well as contact between faces. This class of methods can be easily used to create complex geometries [START_REF] Belytschko | Structured extended finite element methods for solids defined by implicit surfaces[END_REF], however the prescription of boundary conditions on the embedded surface is not straightforward [START_REF] Duboeuf | Embedded solids of any dimension in the X-FEM. Part I -Building a dedicated P1 function space[END_REF][START_REF] Duboeuf | Embedded solids of any dimension in the X-FEM. Part II -Imposing Dirichlet boundary conditions[END_REF]. A particular combination of BF and EB methods shown in Fig. 1(c) deals simultaneously with two or several superposed meshes, which can represent different physics or physical properties. This framework is suited for applications involving fluid-structure interactions (FSI) [START_REF] Baaijens | A fictitious domain/mortar element method for fluid-structure interaction[END_REF][START_REF] Fournié | A fictitious domain approach for Fluid-Structure Interactions based on the eXtended Finite Element Method[END_REF][START_REF] Puso | An embedded mesh method using piecewise constant multipliers with stabilization: mathematical and numerical aspects[END_REF], where the background mesh represents the fluid and the embedded mesh represents the solid. In this work, we place ourselves in the context of continuum solid mechanics and the finite element method to tie two overlapping meshes: one of these is a boundary-fitted mesh and is referred to as the patch mesh, which is fully or partly embedded in a non-boundary fitted host mesh (the physical domain and the mesh do not coincide). Note that there is no restriction on the number of patches that can be embedded into a host domain. The primary issue addressed here is the continuity of fields across the embedded boundaries. The standard methods to impose stiff continuity or Dirichlet boundary conditions along embedded surfaces include the penalty methods, the method of Lagrange multipliers, Nitsche methods and their variants [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Sanders | On methods for stabilizing constraints over enriched interfaces in elasticity[END_REF][START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF][START_REF] Ramos | A new formulation for imposing Dirichlet boundary conditions on non-matching meshes[END_REF]. Here we consider a framework combining the features of the mortar domain decomposition method and the extended finite element methods (X-FEM) to impose continuity constraints. We refer to this framework as the MorteX method.

The X-FEM is an enrichment method based on the partition of unity (PUM) [START_REF] Babuška | The partition of unity method[END_REF]. In this framework, embedded surfaces, cracks or material interfaces can be modeled without explicit mesh conformity. In X-FEM, it is achieved without compromise on the optimal convergence by means of enrichment functions which are added to the finite element approximation using the framework of PUM. The X-FEM methods are extensively used in applications such as fracture mechanics, voids and inclusions modeling, modeling discontinuities in the porous medium, crack propagation using cohesive elements, simulation of shock wave front and oxidation front propagation, and other applications involving discontinuities both strong and weak [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF][START_REF] Diez | A stable extended FEM formulation for multi-phase problems enforcing the accuracy of the fluxes through Lagrange multipliers[END_REF][START_REF] Gross | An extended pressure finite element space for two-phase incompressible flows with surface tension[END_REF][START_REF] Ji | A hybrid extended finite element/level set method for modeling phase transformations[END_REF][START_REF] Savvas | Homogenization of random heterogeneous media with inclusions of arbitrary shape modeled by XFEM[END_REF][START_REF] Faivre | 2d coupled HM-XFEM modeling with cohesive zone model and applications to fluid-driven fracture network[END_REF][START_REF] Ferté | 3d crack propagation with cohesive elements in the extended finite element method[END_REF][START_REF] Sanchez-Rivadeneira | A stable generalized/eXtended FEM with discontinuous interpolants for fracture mechanics[END_REF].

The mortar methods provide us with a comprehensive framework for mesh tying. They are a subclass of domain decomposition methods (DDM) [START_REF] Wohlmuth | Discretization Methods and Iterative Solvers Based on Domain Decomposition[END_REF][START_REF] Gosselet | Non-overlapping domain decomposition methods in structural mechanics[END_REF][START_REF] Keyes | Domain decomposition methods in science and engineering XVI[END_REF][START_REF] Mathew | Domain decomposition methods for the numerical solution of partial differential equations[END_REF], that are tailored for non-conformal spatial interface discretizations [START_REF] Bernardi | A new nonconforming approach to domain decomposition: the mortar element method, Nonliner Partial Differential Equations and Their Applications[END_REF], and were originally introduced as DDM for spectral elements [START_REF] Belgacem | A spectral element methodology tuned to parallel implementations[END_REF][START_REF] Bernardi | Coupling finite element and spectral methods: First results[END_REF]. The coupling and tying of different physical models, discretization schemes, and/or non-matching discretizations along interfaces can be ensured by mortar methods. The mathematical optimality and applicability of the mortar methods in spectral and finite element frameworks were studied extensively for elliptic problems in [START_REF] Bernardi | Coupling finite element and spectral methods: First results[END_REF][START_REF] Belgacem | A spectral element methodology tuned to parallel implementations[END_REF][START_REF] Wohlmuth | Discretization Methods and Iterative Solvers Based on Domain Decomposition[END_REF]. The mortar methods have been successfully adapted to solve contact problems [START_REF] Belgacem | The mortar finite element method for contact problems[END_REF][START_REF] Mcdevitt | A mortar-finite element formulation for frictional contact problems[END_REF][START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Gitterle | Finite deformation frictional mortar contact using a semi-smooth Newton method with consistent linearization[END_REF][START_REF] Farah | A mortar finite element approach for point, line, and surface contact[END_REF].

To ensure displacement continuity across the tied sub-domains, the mortar methods employ Lagrange multiplier fields; as such it is mixed finite element formulation. It is known that the choice of Lagrange multipliers interpolations strongly affects the mesh convergence rate and may lead to loss of accuracy in the interfacial tractions reported for mixed variational formulations as a result of non-satisfaction of Ladyzhenskaya-Babuška-Brezzi (LBB) also called inf-sup condition [START_REF] Babuška | The finite element method with Lagrangian multipliers[END_REF][START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. Issues resulting from the imposition of Dirichlet boundary conditions using Lagrange multipliers methods has been a topic of interest in various domains, such as the classical FEM [START_REF] Barbosa | The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition[END_REF], Interface-enriched Generalized Finite Element Method (IGFEM) [START_REF] Ramos | A new formulation for imposing Dirichlet boundary conditions on non-matching meshes[END_REF], the fictitious domain methods [START_REF] Burman | Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method[END_REF], the mesh free methods [START_REF] Fernández-Méndez | Imposing essential boundary conditions in mesh-free methods[END_REF], etc. This problem has also been dealt extensively within the context of the X-FEM. In [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF][START_REF] Hautefeuille | Robust imposition of Dirichlet boundary conditions on embedded surfaces[END_REF], the authors propose a strategy to construct an optimal Lagrange-multiplier space for the embedded interfaces which permits to apply Dirichlet boundary conditions. As opposed to the strategy of modifying the Lagrange multiplier spaces, the authors in [START_REF] Sanders | On methods for stabilizing constraints over enriched interfaces in elasticity[END_REF] propose a stabilization method to mitigate the oscillatory behaviour of the standard spaces. In the current work, we extend the strategy of modified Lagrange multiplier spaces, which hereinafter will be referred to as coarsening of Lagrange multiplier spaces. This strategy allows us to address specific problems of mesh-locking, which are inherent to mortar methods for overlapping domains, in particular in presence of a strong contrast of material properties in the vicinity of the interface.

Few contributions harnessing the advantages of the mortar method and the X-FEM (but in a different way from what is presented here) are listed below. In [START_REF] Chahine | Etude mathématique et numérique de méthodes d'éléments finis étendues pour le calcul en domaines fissurés[END_REF][START_REF] Chahine | A non-conformal eXtended Finite Element approach: Integral matching Xfem[END_REF] the authors used the mortar methods to ensure weak continuity conditions across the interface between a coarser mesh domain and non-intersecting finer mesh surrounding the crack, which in turn is represented by the X-FEM formulation. The tying in this case is limited to the interface with matching geometries but non-conformal discretizations, which is a classical application of the mortar method. A dual mortar contact formulation integrated into X-FEM fluid-structure interaction approach is introduced in [START_REF] Mayer | 3d fluid-structure-contact interaction based on a combined XFEM FSI and dual mortar contact approach[END_REF]. There, the combined X-FEM fluid-structure-contact interaction method (FSCI) allows to compute contact of arbitrarily moving and deforming structures embedded in a fluid.

The proposed method of coupling mortar and X-FEM competes with the volumetric coupling of the Arlequin method [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF] and the Polytope FEM for embedded interfaces [START_REF] Zamani | Embedded interfaces by polytope FEM[END_REF]. The Arlequin method involves superpositioning of the mechanical states in transition zone, and energy redistribution between these states using weight functions. The Polytope FEM involves the decomposition of elements cut by the embedded interface into new polytope elements. This is achieved by the creation of new degrees of freedom (DoFs) along the interface. In [START_REF] Sanders | A Nitsche embedded mesh method[END_REF], the authors used the Nitsche method to impose tying constraints for the overlapping domains circumventing mesh-locking; this method however requires appropriate stabilization. With the emphasis laid on the interface discretizations handled by the combination of the X-FEM and the mortar methods, many applications could be cited: sub-structuring, inclusion of arbitrary geometrical features into the existing mesh, meshing complex micro-structures, localized mesh refinement near crack tips, general static and dynamic mesh refinement (Fig. 2). However, applications are not limited to mesh tying. In analogy to mortar methods, which were extended to contact problems [START_REF] Belgacem | The mortar finite element method for contact problems[END_REF][START_REF] Fischer | Frictionless 2d contact formulations for finite deformations based on the mortar method[END_REF][START_REF] Fischer | Mortar based frictional contact formulation for higher order interpolations using the moving friction cone[END_REF][START_REF] Popp | Mortar methods for computational contact mechanics and general interface problems[END_REF], the method elaborated here can be used to solve contact problems between a virtual surface (represented by the X-FEM) and an explicitly represent surface of the homologue solid [START_REF] Yastrebov | Numerical methods in contact mechanics[END_REF]. In a separate paper [START_REF] Akula | MorteX method for contact along real and embedded surfaces: coupling X-FEM with the Mortar method[END_REF] we present the extension of this framework to contact problems, which also allows to treat efficiently wear problems.

The paper is organized as follows. In Section 2 we present the formulation of the interfacial mesh tying problem of overlapping domains; its weak form using the method of Lagrange multipliers is derived in Section 3. Section 4 presents the core methodologies of the X-FEM in terms of inclusion/void modeling and of the mortar discretization for mesh tying problems. The effect of the underlying host mesh interpolations on the tying problem is presented as a case study in Section 5. Two main remedies for mesh locking are presented, namely the triangulation of the blending elements and coarse-graining of Lagrange multipliers, which is detailed in Section 6. In Section 7 we consider two patch tests (tension/compression and bending) to emphasize the inherent problems of embedded boundaries in mixed finite element method. In Section 8 we study the mesh convergence properties of the method using the Eshelby inclusion problem. In Section 9 we illustrate the method's performance on various applications. In Section 10 we give the concluding remarks and state the prospective works.

Mesh tying problem

We consider two open domains Ω 1 and Ω 2 with an overlap region Ω 1 ∩ Ω 2 [Fig. 3(a)]. Solid Ω 1 has only outer surfaces which are split into Dirichlet, Neumann and tying boundaries 1 . We refer to the domain Ω 2 as the "host" domain as it hosts the partially embedded domain Ω 1 . In addition to outer boundaries (the Dirichlet boundary Γ 2 u and Neumann boundary Γ 2 t ), the host domain contains an embedded boundary Γ2 g = Γ 1 g ∩ Ω 2 (the subscript "g" refers to the "gluing"). We assume that the solids are glued together along the interface formed by the boundaries Γ2 g and Γ 1 g ; the physics in the overlap zone is determined by solid Ω 1 . Therefore, in the reference configuration Γ2 g = Γ 1 g and must remain so in any configuration. The mesh tying problem is concerned with the enforcement of displacement continuity along the interface Γ 1 g and Γ 1

Γ 1 u , Γ 1 t , Γ 1 g , respectively, such that Γ 1 u ∪ Γ 1 t ∪ Γ 1 g = ∂Ω
g . The standard boundary value problem (BVP) involves ensuring the balance of linear momentum along with the imposed boundary conditions for the two bodies (i = 1, 2):

∇ • σ 1 + f 1 v = 0 in Ω 1 , ∇ • σ 2 + f 2 v = 0 in Ω2 (1)
σ i • n i = ti at Γ i t , (2)
u i = ûi at Γ i u , (3) (4)
where σ is the Cauchy stress tensor, f i v represent the density of body forces, n i is the unit outward normal to Ω i and ti , ûi are the prescribed tractions and displacements. Ω2 = Ω 2 \ Ω1 is the effective non-overlapping region of the host domain volume [Fig. 3(b)], where the bar-notation denotes the open domain united with its closure, i.e. Ω1 = Ω 1 ∪ ∂Ω 1 . The solution is searched in the appropriate functional space

U 0 = {u i ∈ H 2 (Ω i ) | u i = ûi on Γ i u }
, where H 2 is second order Sobolev space. Displacement continuity between the two domains is enforced as:

u 1 = u 2 along
the interface formed by the boundaries Γ2 g and Γ 1 g .

An equivalent continuum problem is depicted in Fig. 3(b).

Weak form

Eqs. (1)-(3) represent the strong form of the standard solid mechanics BVP, and Eq. ( 5) represents the tying constraints. The weak form of the BVP is accordingly divided into the standard structural part and the contribution from the tying. In the weak form, the requirements on the solution functional space U i is slightly relaxed compared to the strong form, and an additional field of test functions is selected from another functional space V i given respectively by :

U i = {u i ∈ H 1 (Ω i ) | u i = ûi at Γ i u } (6)
V i = {δu i ∈ H 1 (Ω i ) | δu i = 0 at Γ i u }, (7)
where H1 (Ω i ) denotes the standard first order Sobolev space. In infinitesimal strain formulation, the virtual work δW s corresponding to the structural part is given by:

δW s = Ω 1 σ 1 : δε 1 dV + Ω2 σ 2 : δε 2 dV δW int - Ω 1 ∪ Ω2 f i v • δu i dV - Γ i t ti • δu i dS δW ext , ( 8
)
where ε is the infinitesimal strain tensor, and δW int , δW ext are the change in internal energy and the virtual work of forces, respectively. Note that the integration of the virtual internal energy of the host domain is carried out only over the effective volume Ω2 . The equality constraints (5) are enforced using the method of Lagrange multipliers, whose functions as well as their variations are selected from the space H -1/2 . The multipliers are defined over the surface Γ 1 g , which was selected from practical consideration of storing Lagrange multipliers on the explicitly represented surface :

λ, δλ ∈ H -1/2 (Γ 1 g ), (9)
where H -1/2 is the fractional Sobolev space defined over the tying boundary. The Lagrangian of the constrained problem leads to a mixed variational saddle point problem which includes both structural δW s and tying δW g contributions. The variation of the Lagrangian is given by:

δL(u, λ) = δW s + δW g = δW s + Γ 1 g λ • (δu 1 -δu 2 ) + δλ • (u 1 -u 2 ) dΓ = 0, (10)
where the integration is carried out over the boundary Γ 1 g . The first term in square brackets in Eq. 10 represents the virtual work contribution from the tying interface and second term represents the weak form of the equality constraints.

Methodology

The evaluation of the internal virtual work restricted to the effective volume of the host solid Ω2 is accomplished with the X-FEM method. The mortar method is extended to enforce the displacement equality constraint over the interface between the overlapping domains, i.e. between the boundary of the embedded domain (patch) Γ 1 g and the corresponding virtual surface Γ2 g of the host solid. The main features of the proposed method are illustrated on an example shown in Fig. 4. It represents the discretized finite-element setting for the overlapping domains. A discretized 1 square patch with a circular hole Ω 0 ∪ Ω 1 with surfaces Γ 0 = ∂Ω 0 and Γ 1 g = ∂ Ω 1 ∪ Ω0 is embedded into a host mesh Ω 2 . As before the bar notation is used to denote the open domain united with its closure, here Ω0 = Ω 0 ∪ ∂Ω 0 . Note that the necessity to consider explicitly the hole as an extra domain comes from the particularities of the problem and was introduced for the sake of avoidance of misinterpretation: the physics inside the contour Γ 1 g is fully determined by the patch with a circular hole. The intersection of the patch's boundary Γ 1 g and the host domain represents the virtual surface Γ2 g = Γ 1 g ∩ Ω 2 . The X-FEM is used to account for the virtual work only in the effective domain volume Ω2 = Ω 2 \ { Ω0 ∪ Ω1 }. The mortar method is is brought into play to tie together the two domains Ω 1 and Ω2 along the interface made of Γ 1 g and Γ2 g . Note that in the presented example, the tying boundary is fully embedded.

Extended finite element method

The virtual surface Γ2 g of the host domain is treated as an internal discontinuity. This is modeled within the X-FEM framework, thereby nullifying the presence of the overlap region Ω1 ∪ Ω0 ∩ Ω2 in the domain Ω 2 [see Fig. 4(c)]. The X-FEM relies on enhancement of the FEM shape functions used to interpolate the displacement fields. Here the enrichment functions describing the field behavior are incorporated locally into the finite element approximation. This feature allows the resulting displacement to capture discontinuities. The subdivision of the host mesh is defined by indicator function φ(X) : R dim → {0, 1} (where X is the spatial position vector in the reference configuration in domain Ω 2 ) [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics[END_REF]. The indicator function is non-zero only in the non-overlapping part of domain Ω 2 :

φ(X) =        1, if X ∈ Ω2 ; 0, elsewhere.

The discontinuity surface Γ2

g can be seen as a level-set defined as follows:

Γ2 g = X ∈ Ω 2 : ∇φ(X) 0
As a result, the indicator function φ(X) partitions the elements of the host domain Ω 2 into three distinct categories [Fig. 5(a)], namely standard elements, blending elements and discarded elements.

In practice, the enrichment of shape functions in case of void/inclusion problem can be simply replaced by a selective integration scheme [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF]. For the standard elements, there is no change in volume of integration and the discarded elements are simply excluded from the volume integration procedure. In order to obtain the effective volume of integration for each blending element, we perform the clipping of the blending elements by the discretized surface Γ2

g [Fig. 5(b)] . The clipping of a single element could result in one or several various polygons 2 both convex, and non-convex , which represent the effective volumes of integration.

To selectively integrate the internal virtual work in the effective volume only, the resulting polygons are virtually remeshed into standard convex elements (for example, triangles). Note that this remeshing is merely performed to use a Gauss quadrature for integration [Fig. 5(c)], and does not imply the creation of additional degrees of freedom, and as such does not change the topological connectivity of nodes. The displacement field is evaluated using the standard shape functions and the original DoFs; only the integration is changed. To carry out this remeshing, we applied the ear clipping triangulation algorithm [START_REF] Mei | Ear-clipping based algorithms of generating high-quality polygon triangulation[END_REF] to the polygons. The DoFs associated with the elements outside the integration domain Ω2 [marked with red crosses in Fig. 5(a)] are removed from the global system of equations.

MorteX discretization

Within the mortar discretization framework, the tied domains are classified into mortar and non-mortar sides. The superscript "1" refers to the mortar side of the interface and "2" to the non-mortar side; the former stores the Lagrange multipliers (dual DoFs) in addition to displacement degrees of freedom (primal DoFs). If the host is selected as a mortar side, the context of the problem becomes similar to the one considered in [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF][START_REF] Hautefeuille | Robust imposition of Dirichlet boundary conditions on embedded surfaces[END_REF], where it was shown that strong restrictions apply on the choice of Lagrange multiplier spaces in order to fulfill the inf-sup condition. The algorithm for construction of such spaces is not straightforward. Therefore, to avoid these difficulties, we select the patch side as the mortar surface, which provides us with a more flexible setting. This choice was already reflected in the fact that the host boundary was chosen as the integration side for tying conditions [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF]. However, under specific problem settings, the choice of employing standard interpolation functions for the Lagrange multipliers on the embedded interface still leads to spurious oscillations of interfacial tractions. Remedies for this problem will be discussed in Section 7.

Displacements on the mortar side Γ 1 g are given by classical one-dimensional shape functions with the interpolation order equal to that of the underlying mesh:

u 1 (ξ 1 , t) = N 1 m (ξ 1 )u m (t), m ∈ [1, M] ( 11
)
where M is the number of nodes per mortar element's edge and ξ 1 ∈ [-1; 1] is the parametric coordinate of the mortar side. The displacements along the virtual surface Γ2 g running through the host mesh elements are characterized by the two dimensional host mesh interpolations, and can be expressed as follows:

u 2 (ζ, t) = N 2 i µ 2 (ζ), η 2 (ζ) u i (t), i ∈ [1, N] ( 12
)
where ζ is the one-dimensional parametric coordinate of integration segment of non-mortar side, µ 2 , η 2 ∈ [-1; 1] are the classical two-dimensional parametric coordinates of the host element, and N is the number of nodes of this element. The use of two-dimensional interpolation on the non-mortar side marks the difference between the classical mortar and the presented MorteX frameworks. The Lagrange multipliers (defined on the mortar side) are interpolated using shape functions Φ:

λ(ξ 1 , t) = Φ l (ξ 1 ) λ l (t), l ∈ [1, L] (13)
where L can be less than or equal to M. It enables to select shape functions for dual variables independently of the primal shape functions. Few remarks could be made here. We need first to introduce the notion of segment: it is a straight line whose vertices can either be a clip intersection 3 or a kink point, the later being a node of the mortar side lying inside the host element. Note that segments are defined both on mortar and non-mortar sides as they coincide. Since a host element can be intersected by several patch segments S i [Fig. 6 (a)], the functions µ 2 (ζ), η 2 (ζ) can be piece-wise smooth, which implies that the underlying displacement (and coordinate) can also be piece-wise smooth. Second remark: if the host elements are quadrilateral, then each segment interpolation is given by p 1 × (p 2 + 1), where p 1 is the mortar interpolation order, and p 2 is the non-mortar interpolation order; the later is augmented by one since the virtual interface passes inside the element, where the interpolation order in quadrilateral 2D elements is one order higher than along the edges. For a triangular host mesh, the interpolation order is simply the product of interpolation orders of the host and patch meshes p 1 × p 2 . To give an example, let us assume that the patch mesh is linear and the host mesh is linear quadrilateral; then let us imagine that a host element is cut into two parts by a (straight) patch segment. Then the displacement in the host element along this straight segment will be second order polynomial function of the parameter ζ [see Fig. 6(b)]. However, the displacement along such a cut would remain linear for triangular linear host elements [see Fig. 6(c)].

Mortar interface element

A mortar element is formed with a single mortar segment (on the patch side) and a single non-mortar element (on the host side). Each tying element consists of (M + N) nodes, M from the mortar segment, and N from non-mortar element, and stores (L × dim) Lagrange multipliers, where (dim) is the spatial dimension. The choice of L is guided by the inf-sup condition requirement of the discrete Lagrange multiplier spaces (usually L ≤ M).

Substituting the interpolations ( 11), ( 12), ( 13) into the weak form [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF] and extracting only the terms related to the mesh tying of a single mortar element, we obtain ( 14)

δW el g = λ l • δu m + δλ l • u m D lm -λ l • δu i + δλ l • u i M li , l ∈ [1, L], i ∈ [1, N], m ∈ [1, M] (15)
where D lm and M ln are the mortar integrals evaluated over the mortar-side segments S el ⊂ Γ 1 g ∩ Ωe i , where Ωe i is the current host element forming the mortar element.

D lm = S el Φ 1 l (ξ 1 )N 1 m (ξ 1 ) dΓ, (16)
M li = S el Φ 1 l (ξ 1 )N 2 i µ 2 (ζ 2 ), η 2 (ζ 2 ) dΓ. ( 17
)
The nodal blocks of the mortar matrices denoted as D (L × M) and M (L × N) can be expressed as:

D(l, m) = D lm I (18) M(l, i) = M li I (19)
where, I is the identity tensor of the spatial dimension of the problem. Using matrix notations, Eq. ( 15) reads

δW el g =           D • L -M • L D • U 1 -M • U 2           •           δU 1 δU 2 δL           (20)
where arrays U 1 , U 2 , L store current values of associated nodal primal (on mortar and non-mortar sides) and dual (mortar) DoFs:

U 1 = u 1 , . . . , u N , U 2 = u 1 , . . . , u M , L = λ 1 , . . . , λ L
whereas their variations are denoted δU 1 , δU 2 , δL. The tangent operator for the mortar interface element is obtained by taking the derivatives of ( 20) with respect to its DoFs:

K =           0 0 D 0 0 -M D -M 0           . ( 21
)
This tangent operator has zero blocks for primal DoFs, non-zero blocks to link primal with dual DoFs, and zero blocks on the main diagonal, which is the typical structure of a saddle-point system.

Evaluation of mortar integrals

The evaluation of mortar integrals is performed over segments S el (see Figs. 6,7). The evaluation of the integrals D lm [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF] is straightforward, as it involves the product of shape functions from the mortar side only. In contrast, the integral M li (17) combines shape functions from both mortar and non-mortar sides. The evaluation of this integral requires a mapping between the mortar and non-mortar sides. The integration domain is parametrized by γ ∈ [-1, 1] (Fig. 7), which needs to be linked with the parametrization of the mortar side, which is given by:

ξ 1 (γ) = 1 2 (1 -γ)ξ 1 a + 1 2 (1 + γ)ξ 1 b , (22)
where ξ 1 a and ξ 1 b define the limits of the integration on the mortar side as shown in Fig. 7. To evaluate the integrals using the Gauss quadrature, we need to find the location of Gauss points γ G in terms of mortar (ξ 1 G ) and non-mortar (µ 2 G , ν 2 G ) parametrization. While the former is straightforward using [START_REF] Ferté | 3d crack propagation with cohesive elements in the extended finite element method[END_REF], the latter can be done by solving the following equation:

N 1 m (ξ 1 G )X 1 m = N 2 i (µ G , η G )X 2 i , (23)
where the physical location of the Gauss point is given by

X G = N 1 m (ξ 1 G )X 1 m .
With these notations, the mortar integrals can now be evaluated using the Gauss quadrature rule as:

D lm = S el Φ l (ξ 1 )N 1 m (ξ 1 ) dΓ = N G G=1 w G Φ l (ξ 1 G )N 1 m (ξ 1 G )J seg (ξ 1 G ) (24)
M li = S el Φ l (ξ 1 )N 2 i (µ 2 , η 2 ) dΓ = N G G=1 w G Φ l (ξ 1 G )N 2 i µ 2 G , η 2 G J seg (ξ 1 G ) ( 25
)
where as previously l ∈

[1, L], m ∈ [1, M], n ∈ [1, N]
and N G is the number of Gauss integration points, w G are the Gauss weights, J seg is the Jacobian of the mapping from the parent space ξ 1 to the real space including the adjustment of the integral limits:

J seg (ξ 1 ) = ∂N m ∂ξ 1 ∂ξ 1 ∂γ X m . (26)

Intra-element interpolation of displacements in the host mesh

For the overlapping domains the coupling is made between the patch mesh boundary and virtual surface running through the host mesh elements. This is reflected in the mortar matrix M (25) that contains the integral of a product of volumetric (in the host mesh) and surface (patch mesh) shape functions. To demonstrate the effect of the interpolation choice, we use the set-up shown in Fig. 8 (a). The following dimensions are used: h 1 = 1 mm, h 2 = 1.25 mm, h * = 0.25 mm and l = 1.5 mm. A uniform pressure of σ 0 = 1 MPa is applied on the top surface of Ω 1 . The meshes Ω 1 and Ω 2 are tied along the interface Γ 1 g . The domain Ω 2 is discretized with a triangular (T1) and two quadrilateral (Q1, Q2) elements, all elements use first order interpolation, and the patch domain is discretized into a rectangular elements [Fig. 8(b)]. The two domains are made of the same linearly elastic material (E = 1 GPa, ν = 0); the reference solution for the selected boundary conditions is a uniform stress field (σ xx = σ xy = 0 and σ yy = σ 0 ) whether it be under plane strain or plane stress formulation. Of course, the displacement along the tying line is uniform. However, as shown in Fig. 9 the selected discretization does not allow to obtain the reference solution. The solid line in Fig. 9 shows the resulting vertical displacement along the tying line Γ 1 g ; it consists of a combination of linear and non-linear portions. We hypothesize that the inability to reproduce the reference solution is somehow related to the interpolation order of displacement in host elements. As seen from the figure, the order of the solution is matching the maximum available interpolation order of displacements p max . This maximum interpolation order is one (p max = 1) along any straight line inside a linear triangle (T1). However, this order raises to two (p max = 2) along straight lines inside quadrilateral elements, as long as both parent coordinates µ and η change along these lines (Q2). This is not the case in Q1 (p max = 1), where one of parent coordinates remains constant along the virtual interface, i.e.

µ 2 (ζ) = const or η 2 (ζ) = const.
This observation motivates us to test triangulation of blending element to limit the maximum interpolation order in host element to one (p max = 1). This operation does not change the number of DoFs and is easy to handle in practice. This procedure enables us to obtain the reference solution in the considered case as demonstrated in Fig. 9. The general applicability of this method will be tested in the following sections.

Coarse Grained Interpolation (CGI) of Lagrange multipliers

The stability of the proposed mixed formulation is guided by the requirement to satisfy the inf-sup condition [START_REF] Babuška | The finite element method with Lagrangian multipliers[END_REF], which is not a trivial task. For example, the use of Lagrange multipliers with standard interpolation leads to nonphysical oscillations along the interface when used to enforce Dirichlet boundary conditions [START_REF] Barbosa | The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition[END_REF]. Having a much x, mm

Q 1 Q 2 T 1
Figure 9: Displacement u y profile along nodes of Γ 1 g,h for the set-up shown in Fig. 8: direct tying using the MorteX method results in non-linear displacement field (solid line); triangulation of blending elements results in a perfect match (dash-dotted line) with the reference solution u y = -1 µm.

stiffer patch than the host material represents a case approximately similar to the imposition of Dirichlet boundary conditions, therefore in most presented examples stiffer patch will be used. The technique presented in [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] for the X-FEM framework, involves coarsening of the Lagrange multipliers to avoid spurious oscillations. This is achieved by algorithmically selecting nodes referred to as "winner nodes" along the Dirichlet boundary. The algorithm favors nodes close to the boundary, or nodes from which many edges intersecting the boundary emanate. The Lagrange multiplier space is built using the winner nodes. The size of the multiplier space is the size of the winner nodes set. Inspired from this technique, we suggest to coarse grain the interpolation of Lagrange multipliers (dual DoFs) to avoid spurious oscillations. In Fig. 10, for example, the number of mortar nodes (each of which carries Lagrange multipliers in the standard approach) per host element considerably outnumbers the associated DoFs to render the associated constraints independent of physical deformations the host elements allow [START_REF] Sanders | A Nitsche embedded mesh method[END_REF]. Therefore the system shown in this figure is overconstrained and requires an appropriate stabilization. Coarse graining of Lagrange multiplier interpolation functions enables to reduce the number of constraints and thus improves the problem stability. In this approach, not every mortar node is equipped with a Lagrange multiplier. Therefore, the interpolation functions become non-local, i.e. they span more than one patch segment. For this purpose, we choose a 1D parametric space ξ CG ∈ [-1, 1], spanning multiple mortar-side segments. Such parametrization can be chosen such that length L i of the corresponding super-segment in the physical space is comparable to the size of host elements. As shown in Fig. 10(a), the mortar-surface is segmented into three super-segments of lengths L 1 , L 2 and L 3 . The end nodes of these segments are termed the "master" nodes (they carry the dual DoFs λ), other mortar-nodes are termed "slave" nodes. We introduce the local coarse-graining parameter κ that determines the number of segments contained in a supersegment, and thus (κ -1) determines the number of slave nodes per super-segment. In Fig. 10(a), the coarse-graining parameter takes the values κ = 4, 9, 5, for the super-segments of lengths L 1 , L 2 and L 3 , respectively.

In theory, the coarse graining is achieved through defining dual DoFs only on master nodes. In practice, this can also be done by keeping the Lagrange multipliers at all mortar nodes and using a multi-point constraint (MPC) to enforce a linear interpolation between the master nodes. Hence, for a given slave node, the dual DoFs are given by:

λ slave (ξ CG ) = Φ CG l (ξ CG )λ l master , l = 1, 2, (27)
The parametrization of the super-segment can be chosen such that ∀i: (ξ CG i+1ξ CG i )/L i+1,i = const, where the numerator represents the spacing between two mortar (patch) nodes in the coarse-grained parametric space and L i+1,i corresponds to the physical length of the corresponding segment. An average ratio of number of mortar segments per number of blending elements, which is termed "mesh contrast" m c , can be used to guide the selection of the coarse-graining parameter κ: for example, m c = 6 (18 mortar segments per 3 host elements) in the example shown in Fig. 10. The coarse-graining parameter κ, for an open mortar surface, can take values in the range κ ∈ [1, N m ], where N m is the number mortar segments; for a closed mortar surface the upper limit is one less.

The optimal choice of coarse-graining parameter κ is studied on particular problem settings in Sections 7 and 8 for open and closed mortar surfaces. The limit case κ = 1 corresponds to the standard Lagrange interpolation (SLI). For the case of approximately regular discretization on both sides, a global coarse-graining parameter can be chosen, and its value is set to be approximately equal to the mesh contrast parameter κ ≈ m c as shown in Fig. 10(b). In case of non-regular mesh discretizations on mortar or/and on host sides, the coarse-graining parameter should be selected element-by-element according to the local mesh contrast as shown in Fig. 10(a). However, in all the examples considered below, we use regular discretizations, and thus the global coarse-graining parameter will be used.

Patch tests

In this section, the algorithms introduced in the previous sections, implemented in the finite element suite Zset [START_REF] Besson | Large scale object-oriented finite element code design[END_REF], are tested on simple problems of tied overlapped domains of different discretizations and different material contrasts subject to bending or tensile/compressive boundary conditions: the patch mesh is Ω 1 and the host mesh is Ω 2 . Linear elastic material properties are used for both the patch (E 1 , ν 1 ) and the host (E 2 , ν 2 ). The geometric set-up of the patch and host domains are illustrated in [Fig. 11 (a)]. The following two extreme cases will be considered: Case 1. a finer and stiffer patch mesh is superposed onto the host mesh, and Case 2. a coarser and stiffer patch mesh is superposed onto the host mesh.

These two particular cases are chosen for the validation since they are prone to severe manifestations of the mesh locking [START_REF] Sanders | A Nitsche embedded mesh method[END_REF] as in the case of enforcing Dirichlet boundary conditions using Lagrange multipliers along embedded surfaces [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF][START_REF] Hautefeuille | Robust imposition of Dirichlet boundary conditions on embedded surfaces[END_REF]. Additionally, the host domain is meshed with "distorted" quadrilaterals which is classical in patch test studies to exacerbate potential anomalies. Moreover, as was shown by a simple example Section 5 the tying along distorted quadrilateral elements is prone to considerable errors if the mortar-type tying is used directly. The material contrast is introduced by choosing E 1 /E 2 = 1000, and E 1 = 1 GPa, both domains have the same Poisson's ratio ν 1 = ν 2 = 0.3. The discretizations for the two cases are shown in Fig. 11 (b,c). The mesh contrast m c ≈ 11 and m c ≈ 0.1 and the number of mortar segments N m = 191 and N m = 35 are used for Case 1 and 2, respectively. Note that all the stress fields σ xx , σ yy along the tying interface, which are presented on the plots below, are found by the data extrapolated from Gauss points and averaged at mortar nodes.

Master mortar nodes with λ Host mesh nodes

Slave mortar nodes with interpolated λ g cutting through the blending elements (shaded in yellow) of the coarser host mesh. The coarse graining of Lagrange multipliers can be implemented with respect to the local (a) or global (b) contrast in mesh densities.

L 1 L 2 L 3 ≈ L 3 ≈ L 2 L 1 (a) (b)

Tension/compression patch test

A uniform pressure σ 0 is applied on the top surface, the bottom surface is fixed in all directions u = 0 [Fig. 12(a)]. This is a classical patch test in contact mechanics, which is used here to test the tying of different materials. This material contrast requires additional lateral conditions (lateral sides are fixed in normal direction u x = 0) to avoid singularities at extremities of the interface. The reference solution for σ yy is a uniform field σ yy = σ 0 . As expected, in case of stiffer and finer patch mesh (Case 1), spurious oscillations are observed. They have large amplitude that reaches 300 % of the reference solution, moreover, they are not confined to the interface but propagate into the bulk [Figs. 13(a), 14(a)]. In case of stiffer and coarser patch mesh (Case 2), the spurious oscillations are of considerably lower amplitude (under 1 %), they are rather localized in the host mesh in close vicinity of the interface and do not extend in the patch mesh [Figs. In order to quantify the improvement achieved with the suggested coarse grained interpolation (CGI) and with triangulation technique, we introduce the L 2 -norm of the error in the σ yy stress component:

E r (σ yy ) = σ ref yy -σ yy L 2 (Γ 1 g ) σ ref yy L 2 (Γ 1 g ) , (28)
where the norm means f (x) -g(x) L 2 (Γ 1 g ) = i f (x i ) -g(x i ) 2 , where x i ∈ [0, L] are the x-coordinate of mortar nodes in the reference configuration, and L is the length of the surface Γ 1 g . In Fig. 15 we demonstrate the performance of the CGI technique. As seen from the figure, the error in stress greatly reduces compared to the standard interpolation (SLI), when coarse-graining parameter κ increases. However, the error saturates at ≈ 10 -3 and the convergence to the reference solution is missing. On the contrary, the triangulation technique (Fig. 16) enables to achieve a superior precision as shown in Figs. 15(a

Bending patch test

A linear distribution of pressure

σ yy = 2σ 0 (x/l -1/2) ( 29
)
with σ 0 = 1 MPa is applied on the bottom surface, while keeping the top surface fixed in vertical direction, only the corner point is fixed in horizontal direction, the lateral sides remain free σ xy = σ xx = 0 [Fig. 12 (b)]. The same linear distribution of the vertical stress component [START_REF] Belgacem | A spectral element methodology tuned to parallel implementations[END_REF] through the two solids should take place. This case study (Case 1, 2) was inspired from the work [START_REF] Sanders | A Nitsche embedded mesh method[END_REF], where the authors also used the combination of the mortar method and the selective integration. It was shown that Case 1, in particular, results in high-amplitude spurious oscillations in the interface contrary to Case 2 that has a smoother stress profile along the interface. Under the standard Lagrangian interpolation set-up we could reproduce similar results, see Fig. 18, 19. These oscillations could be removed with Nitsche method provided some adjustment of the stabilization penalty parameters on each side of the interface [START_REF] Sanders | A Nitsche embedded mesh method[END_REF]. We demonstrate below that using the coarse-grained interpolation for Lagrange multipliers, as suggested in Section 6, also permits avoiding these oscillations in the MorteX framework. As discussed in Section 6, the choice of the optimal coarsegraining parameter is governed by local or global mesh contrast, which can be easily determined either for every segment or for the whole interface. It renders the choice of the coarse-graining parameter κ fully automatic. Contrary to the stabilized Nitsche method, no knowledge about local material contrast is needed. In Fig. 20, the vertical stress component is shown for the Case 1 when the coarse-graining parameter for Lagrange multipliers is set to (a) κ = 6 and (b) κ = 12. It is clearly seen in the figures that κ = 6 does not sufficiently relaxes the over-constraining of the Lagrange multipliers space (we recall that m r = 11) to obtain a smooth reference solution, even though the amplitude of oscillations is slightly reduced compared to standard Lagrange multipliers, which can be seen in Fig. 21(a) where the standard solution (obtained with the standard interpolation of Lagrange multipliers SLI) is compared with the coarse-grained interpolation (CGI). With the coarse-graining parameter κ = 12 we obtain a much improved result comparable to the reference solution.

The relative L 2 error E r (σ yy ) Eq. ( 28) is shown in Fig. 21(b) for different values of κ. The error becomes acceptable only for κ ≥ m c , however, the accuracy of the solution constantly improves with a further increase of κ. The fast drop of the error for κ ∈ [1, m c ] is associated with the graduate removal of spurious oscillations in the stress distribution, whereas for κ > m c improves further the error by better approximation of stresses at extremities of the interface. Since the reference stress distribution is linear, only two Lagrange multipliers are sufficient to capture it, leading to an error reduction up to κ = N m . In general, as will be shown later, too coarse a representation of Lagrange multipliers leads to deterioration of the solution (see Section 8). The triangulation of blending element is also tested in the bending patch test, see Figs. [START_REF] Ferté | 3d crack propagation with cohesive elements in the extended finite element method[END_REF]. In contrast to the compression patch test, triangulation here does not help with the removal of spurious oscillations.

Summary of patch tests

Here, we present the ensemble of patch-test results for various combinations of host/patch meshes and material contrasts. As before, we consider two cases: Case 1 corresponds to a fine patch mesh [mesh-density contrast m c = 10, Fig. 23(a)] which is tied with a coarser host mesh made of triangular, aligned or distorted quadrilateral elements Fig. 23(b,c,d), respectively. Results of Case 1 are presented in Table 1. In Case 2, the patch mesh [mesh-density contrast m c = 0.1, Fig. 24(a)] is coarser than the host mesh, which again can be made of triangles, aligned or distorted quadrilateral elements, see Fig. 24(b,c,d), respectively. Results of Case 2 are presented in Table 2. Softer (E 1 /E 2 = 10 -3 ) and stiffer (E 1 /E 2 = 1000) patch materials are compared to the host material were considered. We also tested different interpolation order (p0 and p1) for standard Lagrange interpolation (SLI), and for coarsegrained interpolation (CGI) in which the parameter takes its maximum value κ = N m . The table clearly demonstrates that the tying performance is strongly dependent on the type of patch test. Cases which show a small error in bending test can demonstrate a slightly higher error in compression test as in case m c = 10, E 1 /E 2 = 1000 for distorted quads with CGI scheme. However, with a high fidelity it could be stated that if the tying method passes the bending patch test then it passes the compression patch test. The inverse is, in general, false. Interestingly, the triangulation of quadrilateral elements can considerably increase the error in case of SLI scheme, it does not happens with the CGI scheme. Clearly, from these tables it can be concluded that the CGI scheme outperforms the standard SLI scheme (both p0 and p1) in all studied combinations of mesh, element types and patch-test type (48 tests in total).

Circular inclusion in infinite plane: convergence study

Having demonstrated a general good performance of the coarse-graining interpolation, here we carry out a meshconvergence study. We focus on the worse case scenario (see Section 7) when the patch mesh is finer and stiffer than the host mesh. We consider a circular inclusion embedded in an infinite softer matrix in plane strain formulation, and subject to a uniform traction applied at infinity [START_REF] Sharma | Circular inclusion in an infinite elastic strip[END_REF][START_REF] Kachanov | Handbook of elasticity solutions[END_REF][START_REF] Hervé | Elastic behaviour of multiply coated fibre-reinforced composites[END_REF]. This particular problem represents a sub-case of a general Eshelby problem of an ellipsoidal inclusion in a matrix [START_REF] Eshelby | The elastic field outside an ellipsoidal inclusion[END_REF][START_REF] Muskhelishvili | Some basic problems of the mathematical theory of elasticity[END_REF]. Fig. 25(a) shows the used computational set-up: a circular inclusion Ω 1 (patch) with radius R = 0.1 mm, centered at origin, is superposed on a matrix Ω 2 (host) represented by a square of side L = 10 mm (L R). Linear elastic material properties are applied to both the inclusion (E 1 , ν 1 ) and the matrix (E 2 , ν 2 ). The inclusion is made more rigid than the matrix by choosing E 1 /E 2 = 1000, E 1 = 1 GPa, the same Poisson's ratio is used for both ν 1 = ν 2 = 0.3. A uniform pressure σ 0 = 0.1 MPa is applied on the right side as shown in Fig. 25(a), displacements on the left side are fixed in horizontal direction u x = 0 and the lower left corner is fixed. The inclusion patch is tied to the host matrix along the boundary of the inclusion Γ 1 g . The analytical solution for the stress state inside and outside the inclusion is given below in polar coordinates (r, θ) [START_REF] Kachanov | Handbook of elasticity solutions[END_REF]. Stress components inside the inclusion (r < R) are given by:

σ - rr = σ 0 2 β 1 + δ 1 cos 2θ (30)
σ - θθ = σ 0 2 β 1 -δ 1 cos 2θ (31)
σ - rθ = - σ 0 2 δ 1 sin 2θ (32)
where,

β 1 = µ 1 (k 2 + 1) 2µ 1 + µ 2 (k 1 -1) , δ 1 = µ 1 (k 2 + 1) µ 2 + µ 1 k 2 . ( 33
Γ * ⊂ Γ 1 g (θ ∈ [-π/2, π/2]).
Outside the inclusion (r > R), the stress components are given by:

σ + rr = σ 0 2 1 -γ 2 R 2 r 2 + 1 -2β 2 R 2 r 2 -3δ 2 R 4 r 4 cos 2θ (34)
σ + θθ = σ 0 2 1 + γ 2 R 2 r 2 + 1 -3δ 2 R 4 r 4 cos 2θ (35)
σ + rθ = - σ 0 2 1 + β 2 R 2 r 2 + 3δ 2 R 4 r 4 sin 2θ (36) (37)
where

β 2 = - 2(µ 1 -µ 2 ) µ 2 + µ 1 k 2 , δ 2 = µ 1 -µ 2 µ 2 + µ 1 k 2 , 2 = µ 2 (k 1 -1) -µ 1 2 - 2µ 1 + µ (k 1 1) . ( 38
)
For the considered plane strain formulation the constants µ 1,2 and k 1,2 are given by

µ 1,2 = E 1,2 2(1 + ν 1,2 ) , k 1,2 = 3 -4ν 1,2 . (39)

Mesh convergence

We are particularly interested in the stress state along the tying boundary = R) where possible spurious oscillations take place. The distribution of the radial stress component σ rr (θ) at r = R can be obtained either from 4Eqs. [START_REF] Bernardi | Coupling finite element and spectral methods: First results[END_REF] or [START_REF] Gitterle | Finite deformation frictional mortar contact using a semi-smooth Newton method with consistent linearization[END_REF]. This analytical solution is compared with the numerical one obtained using MorteX method along the interface Γ * ⊂ Γ 1 g , for which θ ∈ π/2] [see Fig. 25(b)]. For this purpose we use the L 2 error norm as in (28) defined along Γ * :

E r (σ rr ) = ||σ rr -λ rr || L 2 (Γ * ) ||σ rr || L 2 (Γ * ) ( 40
)
where σ rr is the analytical solution, and λ rr is the radial component of Lagrange multiplier vector obtained by projecting it on the radial basis vector of polar coordinates λ rr = λ • e r . The mesh refinement is carried out maintaining a constant mesh contrast between the patch and the host meshes, i.e. the ratio of the mortar segments to the number of blending elements is fixed to be m c ≈ 3 (Fig. 26); the number of mortar segments was varied N m ∈ {256, 512, 1024, 2048} We recall that to represent the interfacial tractions in compression and bending patch tests one or two Lagrange multipliers, respectively, were enough. In contrast to these patch tests, here the stress distribution is no longer affine along the tying interface, therefore it is expected to obtain a more practical result for the selection of the coarse graining parameter κ, which ensures optimal convergence. Results obtained for mesh contrast m c ≈ 6 and for N m = 1024 and different coarse-graining are shown in Fig. 27. The oscillations are clearly seen near the inclusion/matrix interface, especially inside the inclusion. However, for high enough κ these oscillations are completely removed and a uniform stress field is recovered inside the inclusion. A quantitative convergence study is presented in Fig. 28(a). It clearly demonstrates that for standard interpolation (SLI, i.e. κ = 1) or coarse-grained interpolation (CGI) used with small values of κ = {2, 4}, the presence of spurious oscillations induces very high errors in interfacial tractions [see Fig. 29(a)]. When κ = 16, 32, the error reaches its minimum. This is due to the fact that for the given discretization, this level of coarse graining offers an appropriate balance between on the one hand the relaxation of the over-constraining of the Lagrange multipliers, and on the other hand the ability to accurately describe the complex traction field at the interface. For higher values of κ = {64, 128} the error increases again because of too coarse representation of interfacial tractions. It is thus expected that, in general case, there exists a range for κ which ensures oscillation free and accurate enough solution. It is also expected that optimal κ is determined by the global mesh density contrast m c . However, as demonstrated in Fig. 26, the local mesh density contrast can be more pronounced than the average one, therefore it is expected that the optimal value of coarse graining parameter κ lies in the range κ > m c ; for the considered case the error is minimized for κ/m c = {2.667, 5.333}, probably, the optimal value lies in between. The effect of optimal κ is clearly demonstrated in Fig. 29(b) where interfacial tractions λ rr for different κ are plotted.

For the fixed mesh contrast m c = 6 and optimal κ = 16 and sub-optimal κ = 8, the mesh convergence was carried out with meshes of different densities N m ∈ {128, 256, 512, 1024, 2048}. In Fig. 28(b) we plot the error decay with decreasing mesh size, for which we select the length of mortar edge normalized by the total length of the interface h/2πR = 1/N m . For the selected error-measure along the inclusion/matrix interface, the standard interpolation for Lagrange multiplier (SLI) results in optimal convergence (E r ∼ h). However, even though the convergence is optimal, the error remains very high due to the spurious oscillations, implying that an excessively fine mesh would be required to achieve an acceptable error. For example, to reach E r = 0.1 % in case of SLI, 51 200 elements on the mortar side would be needed. Moreover, for quadrilateral host mesh in absence of triangulation of blending elements, the convergence is lost for very fine meshes. In contrast, the coarse graining technique (CGI) used with optimal κ = 16 results in the error below 0.1 %, even with the coarsest mesh used in our study N m = 128. At the same time, the optimal convergence is preserved. As expected, the triangulation of the blending elements slightly deteriorates the quality of the solution, but preserves the optimality of the convergence. The effect of the parameter m c on the amplitude of spurious oscillations is demonstrated in Fig. 30. For a fixed host-mesh discretization we increase the number of mortar edges N m = [64, 128, 256], which correspond to m c ≈ {1, 3, 6}. Fig. 30(a) demonstrates the increase in the amplitude of oscillations with increasing mesh contrast m c for SLI interpolation. The removal of spurious oscillations with the CGI scheme is shown in Fig. 30(b) for reasonable choice of coarse graining parameter κ = {2, 4, 8} for m c ≈ {1, 3, 6}, respectively.

Numerical examples

In this section we illustrate the method of mesh tying along embedded interfaces in light of potential applications. In all the presented examples, we use a linear elastic material model under plane strain assumption. All the problems are solved in the in-house finite element suite Z-set [START_REF] Besson | Large scale object-oriented finite element code design[END_REF]. All triangular and quadrilateral elements used in simulations possess three and four Gauss points for integration, respectively. The triangulated blending elements also use three Gauss points. The MorteX interface uses three integration points to evaluate the MorteX integrals.

Plate with a hole

As a first example, we solve the problem of a square plate with an embedded square patch containing a circular hole, which was used to illustrate the method in Section 4 (Fig. 4). This example demonstrates the ease with which arbitrary geometrical features can be included into the host mesh. Classically, in the X-FEM method a void can be easily included in the host mesh, however, in the vicinity of the void a stronger stress gradients take place, therefore the mesh around the void should be properly refined. It can be easily achieved by surrounding the void with a finer patch mesh, as done here, and by embedding this refined geometry in a coarse host mesh. The geometric dimensions used in the problem are the following: the plate's side is L H = 12 mm, the hole's radius is R = 0.75 mm and the patch's side is L P = 4.5 mm [Fig. 31(a)]. The patch and the host are made of the same material with Young's modulus E = 1000 MPa and Poisson's ratio ν = 0.3. The left edge of the host domain is fixed in x (u x = 0), the lower left corner is fixed, and a uniform traction σ 0 = 1 MPa is applied on the right edge, the upper and lower boundaries remain free σ xy = σ yy = 0. For comparison purposes, a reference solution is obtained with a classical monolithic mesh [Fig. 31(c)]. Fig. 32(a) shows a rather smooth contour of stress component σ xx , which was obtained using MorteX tying with SLI scheme. However, the seemingly smooth stress field near the interface, exhibits oscillations near the interface as can be seen in Fig. 32(b), where σ xx was plotted over a part of the interface Γ * . These oscillations have a smaller amplitude than in cases with high material contrast, and as previously, they can be efficiently removed when coarsegrained interpolation CGI is used, what is shown in the same figure. Coarse-graining parameter κ = 3 appears to be sufficient to remove them. Note that the slight difference between the MorteX tying and the monolithic mesh comes from inherently different extrapolation/interpolation of stresses to the interface nodes.

Crack inclusion in a complex mesh

In many engineering applications, the solids are subjected to cyclic loads and therefore modeling of structures with fatigue cracks appears essential for computational lifespan prediction. The structural finite element analysis can indicate potential locations of the onset of fatigue cracks, however, insertion of cracks is not always trivial [START_REF] Proudhon | 3d simulation of short fatigue crack propagation by finite element crystal plasticity and remeshing[END_REF][START_REF] Feld-Payet | A new marching ridges algorithm for crack path tracking in regularized media[END_REF], especially, in the common case where the original CAD model is not available . Moreover, the position of the onset of the crack is subjected to statistical perturbations, therefore it is often of interest to probe various scenarios in which the crack starts at different locations. Within the proposed framework, studying various fracture scenarios (crack in this case) merely implies placing the patch at a different location on a host mesh, avoiding potential creation of conformal geometries. Here we demonstrate an example of incorporating a crack in a model blade-disk fir-tree connection subject to a vertical tensile load. The frictionless contact is handled using the augmented Lagrangian method in the framework of the standard mortar method. The following dimensions [see Fig. 33(a)] are used for the blade disk assembly: L 1 = 35 mm, L 2 = 12 mm, L 3 = 14 mm and L 4 = 10 mm. The Young's modulus is E = 1000 MPa and ν = 0.3 for the blade, disk, and the patch containing the crack of length a = 0.3 mm. A vertical displacement u y = 0.2 mm is applied on the top surface of the blade. In Fig. 33(a) we present the resulting stress field for the case of intact structure and for the case of a structure with embedded crack, respectively. As seen in the later case, the stress fields are very smooth across the tying interface (shown as white dashed boundary) ensured by MorteX with SLI only. The coarse graining is not needed here as the mesh densities are comparable and the same material is used for the patch and for the host. Similarly to the presented case of crack insertion, the method can be used in general for introducing various other geometric features into the existing mesh. Using the MorteX method, the location/orientation of these features can be adjusted with ease and without remeshing to perform a sensitivity analysis.

Multi-level submodeling: patch in a patch

In this example we demonstrate the ability of the MorteX method to handle multi-level overlapping domains, i.e. when an embedded patch mesh hosts other domains. In Fig. 34(a) we present such a scenario where a patch with a notch (Ω 2 ) is embedded into a host domain Ω 1 , both made of the same material. At the same time, the patch Ω 2 itself hosts 2 circular inclusions (Ω 3,4 ), which are stiffer than the surrounding material. The following dimensions are used: , Ω 3 , Ω 3 , and Ω 4 respectively). A Poisson's ratio of ν = 0.3 is used for all the domains. A vertical displacement u y = 0.1 mm is applied on top surface of the Ω 1 , while the left surface is fixed in the x direction and bottom is fixed in all directions. The contour plots of σ yy for the cases of SLI and CGI are shown in Fig. 35(a) and (b), respectively. The oscillations in the stress are distinctly seen in case of SLI, but they are removed by applying CGI with κ = 4. Fig. 36 compares σ yy along Γ 3 * , Γ 4 * , which form π/2 portions of matrix/inclusion interfaces. This example illustrates the case where a host mesh with embedded domains of different material properties can be dealt within CGI MorteX scheme. Note that in contrast to the Nitsche method where the stabilization parameter needed to avoid mesh locking is dependent on the local material contrasts [START_REF] Sanders | A Nitsche embedded mesh method[END_REF], the CGI stabilization does not require the knowledge of this contrast. In the CGI scheme, knowing a local or a global contrast of mesh densities across the tying interface m c is enough to automatically select the coarse graining parameter κ, which efficiently stabilizes the mixed formulation and removes spurious oscillations present in the standard mortar scheme.

Conclusion

We presented a unified framework for mesh tying between overlapping domains. This framework was entitled MorteX as it combines features of the mortar and X-FEM methods. As known, the resulting mixed finite element problem may be prone to mesh locking phenomena especially for high material or mesh-density contrasts between the host and the patch meshes. Manifestation of the emerging spurious oscillations for different element types and various material as well as mesh contrasts was illustrated on two patch tests (bending and compression) and on selected examples. These oscillations strongly deteriorate solution in the vicinity of interfaces resulting in poor mesh convergence. Even though triangular elements help to avoid oscillations in compression patch tests, they do not perform well in bending patch test, nor in more complicated examples. These oscillations comes from the overconstraining of the interface in case of mesh-density contrast, when few mortar-side nodes located on the patch mesh are tied to displacement field of a single host element. To get rid of the resulting mesh locking, we suggested to coarse-grain interpolation (CGI) of Lagrange multipliers by interpolating the associated field along few mortar edges. It implies that only every (κ + 1) node along the mortar side stores a Lagrange multipliers and a linear interpolation is used in between. The value of coarse-graining spacing parameter κ controls the performance of the scheme. If κ is too small compared to mesh-density contrast, the spurious oscillations persist as in standard interpolation of Lagrange multipliers (SLI). If the value of κ is too high, the spatial variation of resulting interface tractions cannot be captured properly. Therefore, in general problem, there exists an optimal choice for the spacing parameter κ which can be automatically determined either by local mesh-density contrast between the patch and the host mesh or by the global mesh-density contrast. The performance of the MorteX method with coarse-grained interpolation was demonstrated on Eshelby problem for a stiff inclusion in a softer matrix (elastic contrast of 1000). Few other examples, demonstrating the ease with which the method can be used for: submodeling, local mesh refinement and inclusion of arbitrary geometrical features in the existing mesh, without remeshing. Among these examples, a multilevel/hierarchical overlapping is shown, where a patch is inserted into a host mesh, which in turn is inserted into another host mesh. The MorteX method equipped with CGI demonstrates a very good performance, removes the mesh locking oscillations, and ensures optimal convergence. The important feature of the method is that its stabilization requires knowledge of local mesh densities only, thus it presents a good alternative to the Nitsche method, which requires stabilization constructed on a priori knowledge of material stiffness in the interface. In analogy with the classical mortar method, the MorteX method can be extended to handle contact problems along virtual interfaces embedded in a mesh; this extension is presented in a separate paper [START_REF] Akula | MorteX method for contact along real and embedded surfaces: coupling X-FEM with the Mortar method[END_REF].

Figure 1 :

 1 Figure 1: Illustration of meshes with (a) fitted boundary, (b) embedded boundary (for example, level set), (c) embedded mesh.

Figure 2 :

 2 Figure 2: Few applications: (a) substructuring and inclusion of arbitrary geometrical features, (b) microstructure modeling, (c) localized mesh refinement, for example, around an inserted crack.

Figure 3 :

 3 Figure 3: (a) Continuum setting of the two overlapping domains with applied boundary conditions; (b) an equivalent continuum problem without overlap.

Figure 4 :

 4 Figure 4: (a) Two overlapping meshes: the host Ω 2 and the patch Ω 1 with a circular hole Ω 0 are tied along interface Γ 1 g ; (b) zoom on the interface between the host and patch meshes; (c) effective volume of the host mesh Ω2 = Ω 2 \ { Ω0 ∪ Ω1 }.

Figure 5 :

 5 Figure 5: (a) Element classification in X-FEM framework; (b) clipping of blending elements by Γ2 g , the volume colored in blue in (b.1-3) is the effective volume of integration ( Ωe ): (b.1) a convex polygon, (b.2) a non-convex polygon, (b.3) disjoint polygons; (c) selective integration is carried out over re-triangulated blending elements with reinitialized Gauss integration points (shown in red).

Figure 6 :

 6 Figure

Figure 7 :

 7 Figure 7: Example of mortar integration domain defined by segment S el connecting a clip intersection ξ 1 a and the kink point ξ 2 b , Gauss points are shown as triangles.

Figure 8 :

 8 Figure 8: Illustration of the effect of underlying mesh interpolations: (a) problem setting: tying of two rectangular overlapping domains, (b) discretized patch (upper solid) and host (lower solid) domains, (c) triangulated effective volume of the host domain Ω2 ; note that in (b,c) the two meshes are shown separated only for convenience.

Figure 10 :

 10 Figure 10: An example illustrating the embedded surface Γ 1 g cutting through the blending elements (shaded in yellow) of the coarser host mesh. The coarse graining of Lagrange multipliers can be implemented with respect to the local (a) or global (b) contrast in mesh densities.

  13(b), 14(b)].

  ) and 17.

Figure 11 :Figure 12 :

 1112 Figure 11: Validation tests set-up: (a) problem setting h 1 = 1.0 mm, h 2 = 1.25 mm, h * = 0.25 mm, l = 5.0 mm, the elastic contrast between the patch and the host is given by E 1 /E 2 = 1000; finite-element discretizations of the patch and host solids are shown in (b) for Case 1 ( the patch mesh is finer than that of the host, m c ≈ 11, N m = 191), and in (c) for Case 2 (the host mesh is finer than that of the patch, m c ≈ 0.1, N m = 35).

Figure 13 :Figure 14 :

 1314 Figure 13: Compression patch test: contour plots of stress component σ yy in (a) Case 1, and (b) Case 2.

Figure 15 :

 15 Figure 15: Compression patch test: (a) decay and saturation of the relative error E r (σ yy ) for κ = {1, 6, 12, 24, 48, 96, 192} for CGI in comparison with the error obtained with triangulation of blending elements; (b,c) contour plots of stress component σ yy for (b) κ = 12 and (c) κ = N m = 192.

gFigure

  Figure Patch test stress σ yy with Triangulation of blending elements: (a) Case 1, (b) Case 2.

Figure 18 :Figure 19 :

 1819 Figure 18: Bending stresses σ yy : (a) Case 1, (b) Case 2.

Figure 20 :

 20 Figure 20: Coarse grained Lagrange multiplier space for Case 1: (a) κ = 6, (b) κ = 12.

Figure 21 :

 21 Figure 21: (a) Case 1: comparison of bending stresses (σ yy ) for SLI and CGI with the reference solution; (b) decay of the relative error E r (σ yy ) for κ = {1, 6, 12, 24, 48, 96, 192} for CGI in comparison with the error obtained with triangulation of blending elements.

Figure 22 :

 22 Figure 22: Bending stresses σ yy in Case 1 for triangulated blending elements: (a) stress distribution along the tying interface; (b) contour plot of the stress component.

Figure 23 :

 23 Figure 23: Discretized setting of coarse host mesh and finer patch mesh (m c ≈ 10): (a) patch mesh (Ω 1 ), (b) host mesh with linear triangular elements (Ω 2 ), (c) host mesh with bilinear quadrilateral elements (Ω 2 ), (d) host mesh with bilinear distorted quadrilateral elements (Ω 2 ).

Figure 24 :

 24 Figure 24: Discretized setting of fine host mesh and coarse patch mesh (m c ≈ 0.1): (a) patch mesh (Ω 1 ), (b) host mesh with linear triangular elements (Ω 2 ), (c) host mesh with bilinear quadrilateral elements (Ω 2 ), (d) host mesh with bilinear distorted quadrilateral elements (Ω 2 ).

Figure 25 :

 25 Figure 25: Square matrix with circular inclusion: (a) problem setting (not to scale) with inclusion domain Ω 1 (radius R = 0.1 mm) superposed over the matrix domain Ω 2 (L = 10 mm), (b) inclusion and matrix domain discretizations, where Γ 1 g is the tying boundary and Γ * ⊂ Γ 1 g (θ ∈ [-π/2, π/2]).

(Figure 26 :

 26 Figure 26: Illustration of the local mesh density contrast. The ratio of mortar segment per number of blending elements is m c ≈ 3, but locally the number of mortar segments intersecting different blending can vary considerably m loc c = 3, 6, 6 in zooms shown on the right.

Figure 27 :

 27 Figure 27: Stress component σ xx computed using MorteX method: (a) standard interpolation (SLI), (b) coarse grained interpolation (CGI) for κ = 4, (c) CGI for κ = 16.

Figure 28 :

 28 Figure 28: Convergence study results for the circular inclusion problem for mesh contrast m c ≈ 6: (a) E r (σ rr ) error change with coarse-graining parameter κ, (b) comparison of convergence of SLI and CGI with and without triangulation of blending elements, the mesh size along the interface is normalized by the circumference h/2πR = 1/N m which is equivalent to the inverse of the number of mortar segments.

Figure 29 :

 29 Figure 29: Comparison of λ rr with analytical solution for various values of κ along Γ * (for the mesh with N m = 1024): (a) the standard Lagrange multiplier spaces (κ = 1), (b) the coarse grained Lagrange multiplier solution for κ ∈ [8, 32, 64, 128, 256].

Figure 30 :

 30 Figure 30: The effect of m c on λ rr along Γ * : (a) with standard Lagrange multiplier spaces (κ = 1), (b) with coarse grained Lagrange multiplier spaces (κ≈m c ).

Figure 31 :

 31 Figure 31: Plate with a hole: (a) overlapping domain setting, (b) discretized overlapping domains, (c) monolithic discretization used to obtain the reference solution.

Figure 32 :

 32 Figure 32: Regular host mesh with an embedded patch containing a circular hole (m c ≈ 3): (a) contour plot of σ xx stress component obtained using MorteX (SLI scheme), (mesh is shown only in the second quadrant, the interface is marked with a dashed line; (b) comparison of σ xx distributions along Γ * between the solution obtained with a monolithic mesh [Fig. 31(b)] and MorteX solution obtained with overlapping meshes [Fig. 31(a)] using SLI and CGI schemes.

Figure 33 :

 33 Figure 33: Problem setting for the model blade-disk assembly with an embedded patch mesh with a crack: (a) shows stress component σ yy contour plot and a zoom near the tying interface; (b) shows used host and patch mesh.

L 1 =

 1 5 mm, L 2 = 3 mm, L 3 = 1 mm, L 4 = 3 mm, R 3 = 0.2 mm, and R 4 = 0.4 mm. The material properties used are: E 1 = 1.0 MPa, E 2 = 1.0 MPa, E 3 = 100.0 MPa, E 4 = 1000.0 MPa (the upper indices correspond to the domains Ω 1

4 Figure 34 :

 434 Figure 34: (a) Multi-level overlapping domain set-up; (b) finite element meshes to be coupled (shown not in proportion), Γ 3 * , Γ 4 * denote interfaces over which CGI and SLI are compared in Fig. 36.

Figure 35 :

 35 Figure 35: Contour plots of σ yy for multi-level overlapping domains: (a) SLI, (b) CGI used with κ = 4.

Figure 36 :

 36 Figure 36: Comparison of σ yy stress for standard (SLI) and coarse grained (CGI) interpolations along : (a) Γ 3 * , (b) Γ 4 * .

Note that for brevity and simplicity hereinafter we preserve the same notations for discretized entities as were introduced in the continuous problem statement.

Hereinafter, we assume that all elements use first order interpolation, therefore all edges of elements are straight. It enables us to assume that an intersection or difference of elements can be always represented as one or several polygons.

The clip intersection points are located on the edges of blending elements intersected by the virtual surface Γ2 g .

Note however that σ θθ is not continuous across the interface

Acknowledgment

The authors acknowledge financial support of the ANRT (grant CIFRE no 2015/0799). We are also grateful to Nikolay Osipov and Stéphane Quilici from the Z-set team for their help with the implementation of the method.