On the construction of approximation space to model discontinuities and cracks with linear and quadratic extended finite elements

This paper presents a robust enrichment strategy to model weak and strong discontinuities as well as cracks for industrial applications. First, numerical issues encountered with popular extended finite element approximation spaces are pointed out. Then, the paper gives indications on how to circumvent those issues. The very originality of the paper relies on questioning the theoretical approximation spaces with respect to numerical results and to modify accordingly their design. The relationship between the new design and the previous designs is clearly established, in order to highlight the very small implementation cost of the modifications exposed here. Hence with minimal additional computational cost, gains in accuracy can be significant as shown later in the paper.


Introduction
Strain localization is usually an issue for conventional finite element approaches due to numerical issues in the softening regime of the stress-strain relation, when the problem becomes ill-posed. Apart from taking care directly of the localization by an adaptation of the mesh to the discontinuity, different methods have been used in the literature to circumvent this difficulty. Smeared-cracked models were first proposed [1,2] with perturbations of the fields across the interface. Discontinuous embedded elements appeared almost at the same time with initial papers of Ortiz et al. [3] or Dvorkin et al. [4], with arbitrary orientations through an element but independent from an element to its neighbor. A little bit later, special interface elements [5,6] were proposed localized in between conventional elements, which require frequent re-meshing and refined meshes in order to allow for crack propagation in the correct direction. The eXtented Finite Element Method X-FEM [7] and the Generalized Finite Element Method GFEM [8] finally allow meshes not to respect the crack geometry while providing a continuous transfer of information from one element to the next one about the crack surface localization unlike the generalized class of embedded discontinuous finite element approaches [3,4]. These methods with nodally based enrichments have managed to combine performance and robustness, considering non-meshed cracks in a finite element framework. X-FEM and GFEM use the Partition of Unity [9] and enrich the classical basis of shape functions with discontinuous functions [10]. The discontinuity of the displacement field across the crack surface is then introduced by a generalized Heaviside function, and adding asymptotic fields at the front crack gives good precision in linear elastic fracture mechanics [10,11]. The main advantage of these methods in comparison to mesh-less methods is their easy implementation in a general finite element software, and their capabilities to be applied to various fields: large transformations [12] or plasticity [13] in the X-FEM context for example… One can say that X-FEM and GFEM extend the possibilities of the FEM, while keeping all its advantages. A useful amelioration has been proposed by Sukumar et al. [14] with the introduction of level set functions to represent discontinuities (cracks, voids,…). These approaches are extremely handy in 3D to treat crack propagation [15,16].
Moreover, these methods have been implemented to solve linear elasticity fracture mechanics problems [8,17] with better accuracy than with finite element methods. However, the extension of those methods to quadratic elements in three dimensions is a big challenge. Meanwhile, industrial numerical tools use extensively quadratic elements and simulations in three dimensions are almost the norm. In the wake of [18][19][20], this article attempts to close in the gap between the theoretical methods and issues related to implementation in industrial software.
Minor concerns in one dimension with linear elements turn out to be daunting challenges in 3D or with quadratic elements. Those concerns have been described in [17,[21][22][23] but have raised little interest since recent publications [24,25]. Firstly, they are related to the ill conditioning of strong or weak discontinuous approximations in the general case of non-conforming interfaces and secondly, to the numerical issues related to "geometrical enrichment" techniques, near the crack-tip.
The analysis will be developed more particularly in the case of strongly discontinuous approximations because a direct link with conditioning can be clearly established. For weak discontinuities involved in bi-materials for instance, conditioning issues are still present but they are coupled with the quality of the approximation space to represent continuous solutions with discontinuous gradients [14,[26][27][28]. A specific section will be dedicated to this analysis.
In "Strong discontinuity approximation conditioning" and "Singular enrichment space optimality and conditioning" sections, conditioning issues related to strong discontinuity and singular functions are investigated, as well as, strategies available in the literature to solve those issues. In "Approximation spaces in the literature" and "Numerical behavior of strongly discontinuous and singular approximations" sections, those strategies are benchmarked with linear and quadratic elements. Further analysis unfolds that quadratic elements emphasize conditioning and accuracy issues almost unseen with linear elements. The results exposed here are general enough, given the wide range of approximation spaces considered in the paper.

Strong discontinuity approximation conditioning
In mechanical engineering general purpose software with extended finite elements, strong discontinuities are positioned rather arbitrarily within the bulk of the mesh. A sensitivity study is clearly needed to check out whether strongly discontinuous enrichment strate-gies are dependent on the position of the interface or not (see Fig. 1). In the literature, the answer is far from obvious. Authors suggested that there is indeed a sensitivity of the solution with respect to the position of the interface; nonetheless, they considered this sensitivity as a numerical side-effect and complex numerical solutions have been elaborated to deal with this issue [21,24]. Although those techniques may work with linear elements, very little is said about their efficiency with quadratic elements. As a matter of fact, the asymptotic behavior of strongly discontinuous approximations is not well understood. Conditioning and accuracy of solutions with strong discontinuities deteriorate severely when the interface gets close to the vertices of the mesh. For curved interfaces or unstructured meshes, the interface has a strong probability of getting close to one vertex at least. Dealing with this configuration is a requirement to prevent unexpected results when switching meshes or when moving the interface during industrial numerical simulations.
Those numerical issues have been studied in the case of X-FEM formulations [24]. With X-FEM, the enrichment functions and the classical shape functions are almost collinear when the interface cuts the mesh near one of its nodes (Fig. 2). Whenever those configurations appear, conditioning deteriorates very quickly. In such cases, there are plenty of strategies to limit the condition number large increase. Those strategies are detailed below. Fig. 1 Positions of the interface investigated in the paper. The condition number is optimal in two cases: firstly, when the interface cuts through the middle of the elements and secondly, when the interface is on the boundary of the elements. Apart from those cases, the condition number soars. This behavior will be investigated later in the paper Fig. 2 Example of configurations leading to ill conditioning. The speed of increase of the condition number is related to the distance between the interface and the nodes of the mesh. The condition number soars when either the distance N 1 IP goes to zero or when N 2 IP goes to zero First of all, we can consider the elimination of degrees of freedom, based on straightforward criteria [18,24]. Those criteria are rough estimates of the "Heaviside information" on both sides of the discontinuity. Basically, if the Heaviside information is heavily unbalanced, we have: -either: -or: where, Supp(Φ i ) refers to the support of the shape function and + and − are the domains defined on Fig. 3. Practically, this means that node numbered i is either on the side − or on the side + and does not "see" the information from the opposite side. From now on, the information from the opposite side will be called complementary information.
The criterion used here replaces the measurement of volumes [18,24] by the measurement of distances between the nodes of the mesh and the interface on cut edges. When these criteria are associated to a relocation of the interface they are usually named fit-to- Level-sets (abbreviated "lsn") are a handy tool to evaluate those lengths, since they do not require computing the coordinates of the intersection points.
From the stand of node N 1 , we have if lsn refers to the normal level set: where N 1 and N 2 are the vertices of the linear element pictured (Fig. 2). IP is a given intersection point on the linear element. lsn(x) represents the level-set function at point x.
Then, the level-set values of both nodes are compared in order to enrich or not node N 1 . A similar approach consists in changing directly the value of the level-set of the N 2 node to zero. In any case, a threshold is needed to decide whether or not the enrichment needs to be modified. This threshold based on the ratio of characteristic dimensions is generally chosen between 10 −2 and 10 −3 [18,24].
A simple criterion expresses as: • Node N 1 is eliminated or the level set is moved to N 2 by resetting lsn (N 2 ) to zero if, • Node N 2 is eliminated or the level set is moved to N 1 by resetting lsn (N 1 ) to zero if, Apart from distance or volume weighting criteria associated or not to fit-to-vertex, a more evolved criterion may be used [24]. Here, the additional idea is about weighting the "balance of Heaviside information" inside the stiffness matrix, instead of weighting the geometrical information of physical domains. This "stiffness" criterion attempts to correlate the condition number of the "stiffness" matrix to the numerical threshold of elimination. On second hand, to control the condition number, one can consider an algebraic preconditioner [21], dedicated to orthogonalize shape functions and enrichment functions, with a Cholesky factorization.
The blocks of P C express as follows: where the inverse terms are computed from Cholesky's factorization of local blocks of K :

Singular enrichment space optimality and conditioning
In linear elasticity, the asymptotic displacement solution at the tip of the crack satisfies (in the local basis defined, Fig. 4) [10,11]: where κ = 3 − 4υ (under the plane strain hypothesis) is Kosolov's constant and υ is Poisson's ratio. The first eigenvalue λ 0 = 0.5 represents the less regular term of this expansion series. The associated functions (∝ √ r) are responsible for a significant loss of accuracy within regular FEM [17].
It is the focus point of the X-FEM enrichment to recover at least an order of convergence in energy norm of min (1.5, m), where m is the interpolation order of the elements (1.5 is related to the regularity of the next eigenvalue λ 1 ). For linear elements, an order of convergence close to m = 1 is expected in the energy norm. For quadratic elements, an order of convergence close to m = 2 is expected in the energy norm, when only the components related to the 0.5 eigenvalue are activated by the boundary conditions.
Laborde et al. [17] showed that a fixed enrichment zone is needed, around the cracktip, to ensure the optimal accuracy of singular approximations (see Fig. 5). This technique uncouples the size of the enrichment zone and the size of the mesh, so that the enrichment zone does not shrink to zero with mesh refinement.
However, this enrichment technique has given rise to many concerns almost from the very beginning: • The first concern is the definition of the enrichment area and the related enrichment strategy in blending elements. In the literature, two main strategies emerge: the "cutoff" strategy of Nicaise et al. [29] and the geometrical enrichment strategy [17]. Firstly, the "cutoff" strategy adds global singular d.o.f. and an additional function to soften the transition between the enriched zone and the non-enriched zone. Secondly, the geometrical enrichment introduces local singular d.o.f. which values are set to zero outside the enrichment zone. Nicaise and et al. [29] shows that both strategies are relevant and do not alter the convergence rates. • The second concern is the swift increase of the condition number. Latest works of Chevaugeon et al. [30] and Guptaa et al. [25] give strong research directions to improve the condition number. A vectorial enrichment reduces drastically the number of d.o.f. needed to describe the singular solution. Guptaa et al. [25] suggests that vectorial enrichment could be further improved to permanently remove conditioning issues. The idea of [25] is to subtract from the enrichment function its linear interpolation on the set of elements where the singular enrichment is defined.

Approximation spaces in the literature
In the literature, there are numerous approximation spaces to model strong discontinuities and cracks within the general framework of extended finite element methods and alike.
The earliest work of [7] introduced the standard X-FEM approximation space of Table 1 in 1999. Then in 2000, the work of [8] extrapolated GFEM [22] to crack modeling. Although X-FEM and GFEM both aimed to model cracks, GFEM was designed to deal with a wider scope of problems than X-FEM. Initially, X-FEM and GFEM were assumed as different approaches. X-FEM introduced scalar functions to model Irwin modes at the crack-tip. On the opposite, GFEM used directly Irwin modes as vector functions to model the singular behavior at the crack-tip. In "Approximation spaces for a cracked domain" section, this difference is investigated. We will find out that X-FEM and GFEM are closely related. In 2004, Laborde et al. [17] paper questioned the numerical accuracy of the X-FEM approximation, particularly with higher order elements. Laborde et al. [17] suggested that X-FEM with one layer of enriched elements at the crack-tip, is not very accurate. When the enrichment zone enlarges, accuracy improves, but conditioning deteriorates. Therefore, it suggested a new approximation space to take care of conditioning issues and with a better behavior for quadratic elements. As this new approximation space is restricted to 2D analysis, this approximation is not fitted for industrial numerical simulations. Finally, to address higher order modeling, X-FEM design evolved into vectorial enrichment [30], close to GFEM.
Very recently, in 2013, even the GFEM approximation space has evolved into a new space called SGFEM [25,31]. SGFEM addresses conditioning issues of GFEM with a new definition of the vectorial functions. Following [32] it is even defined by the fact that the condition number of the associated stiffness matrix has to be of the same order with respect to mesh refinement than the one of FEM. However, SGFEM as proposed in [25,31] has some drawbacks that will be discussed in "Approximation spaces for a cracked domain" section. Specifically, to model strong discontinuities, many approximations have been released in the literature. The earliest work of [7] introduced a Heaviside function to model the jump d.o.f. on nodes in the vicinity of a strongly discontinuous interface (Fig. 3). The next relevant approximation space was the one of Hansbo et al. [33] that introduces discontinuous polynomials within the partition of unity. Numerous spaces were derived from those previous approximation spaces and far too many to be studied thoroughly in this paper [19,20,34,35]. Hence, we will focus, in "Strong discontinuity representation" section, on the ones of Moës et al. [7], Hansbo et al. [33] and Belytchko et al. [36], which is somehow intermediary between those of [7] and [33].
From the mathematical point of view, we will find out that the three spaces are very similar.

Strong discontinuity representation
In this section, we model only a strong discontinuity, with a media fully split by an interface. Hence, we do not need the singular functions in that case. Singular functions will be investigated in the next section.
Thus, let us assume the following X-FEM approximation space V h , modeling a strong discontinuity, where Φ i are lagrangian polynomials up to order two, in the scope of this article. I denotes the nodes of the finite element mesh, I H the set of enriched nodes and H the Heaviside function. A node belongs to I H if, and only if, its support is split by the interface . We recall that H is defined as follows: Let's consider the basic case where the whole domain is divided into two parts as in Fig. 3: For which, • u h + represents the approximation of the solution on the dimain + extended by continuity to the domain − over the interface , • u h − represents the approximation of the solution on the dimain − extended by continuity to the domain + over the interface .
The representation of the kinematic jump is crucial for many interfacial laws, such as cohesive laws and contact laws. The jump at every point (commonly integration point) located on the interface, follows the convention: In order to simplify the notations in Table 1, we define for a given node j ∈ I H , j,1 = χ + j and j,2 = χ − j , where χ A denotes the characteristic function of the set A. Using these notations, we have: Linear shape functions on enriched nodes Displacement approximation Jump approximation X-FEM [7] approximation space is on the left, Hansbo et al. [33] is in the middle, and Fries [26] approximation space is on the right where x j is the location of j. The function H − H x j is used to define the formulation introduced in [26]. Let us compare side-by-side, the XFEM representation of the displacement jump introduced by [7] with the other well-known formulations of [33] and [26]. Even though the jump approximations have different expressions (Table 1), we will investigate next whether the approximation spaces are really different or not.
Areias et al. [37] showed that [7] and [33] involve different bases but that both enrichments still represent the same approximation space. Each enrichment can be expressed in terms of the other one with a suitable change of variables. With the same analysis as in [37], the different enrichments can be expressed in terms of the other ones for the remaining couple of formulations X-FEM/Shifted and Domain/Shifted. Thus, from the mathematical point of view, the approximation spaces involved in the three formulations are equivalent (see Table 2). The same can be said of other formulations encountered in the literature [19,34,35]. We would like to stress the fact that this equivalency only holds in the case of strong discontinuity. When singular functions are injected in the approximation, the relationship between the different approximation spaces may be different as studied in the next section.

Approximation spaces for a cracked domain
The model problem is a cracked domain (Fig. 6), under a linear elasticity assumption. The material is also assumed to be homogeneous and isotropic. Dirichlet boundary conditions are applied on the boundary D and Neumann boundary conditions are applied on N .
The space of admissible displacements is:  The weak form of the equilibrium problem is: where u is the displacement, σ is the Cauchy stress, ε is the strain, 1 the identity tensor, f is the body force applied on ω, g is the traction applied on N and λ and μ are the Lamé parameters. Now, the discrete approximation spaces from the literature will be investigated. In the following sections, we will prove that spaces available in the literature can be gathered into two classes: "straightforward" enrichment and "bubble" enrichment. An illustration of those two classes is provided on Fig. 7. "Straightforward" enrichment means that the singular function is introduced directly into the approximation space without modification. "Bubble" enrichment means that the singular function is reshaped into a new function before its introduction into the approximation space.
In the following sections, those classes are described: • "Straightforward" singular enrichment class" section shows that X-FEM and GFEM are "straightforward" enrichments. Moreover, GFEM is a subspace of X-FEM. This statement implies that X-FEM is somehow more accurate than GFEM, but both methods are nonetheless very close. • "Bubble" enrichment class" section shows that a "bubble" enrichment as the one of Gupta et al. [25] is not a "straightforward" enrichment. This statement implies that "bubble" enrichments have numerical properties not identical to "straightforward" enrichments that will be investigated in "Singular approximation at a crack tip" section.

"Straightforward" singular enrichment class
Let V h be the approximation space using vector asymptotic functions to describe the crack kinematics: where d is equal to two in the bidimensional case and three in the tridimensional case, E 1 , E 2 , E 3 is the usual Cartesian basis, K α are the vector asymptotic functions and I CT the set of nodes enriched by such functions. In this contribution, I CT is defined as the set of nodes placed at a distance lower than R from the crack-tip (cf. Fig. 5).
The vector asymptotic functions K α are defined as follows: with u 0 I , u 0 II u 0 III given by (15). Let V h CUTOFF (W 2,∞ ) be the approximation space defined with a cutoff function χ ∈ W 2,∞ [29], which value is set to one if its distance from the crack front is lower than r 0 , is set to zero if its distance from the crack front is greater than r 1 with r 1 > r 0 and varies linearly in between: and let V h XFEM be the standard X-FEM approximation space [7]: where F α are defined as: We remark that: Consequently, the vector functions K α and the scalar functions F α are related by the following linear relations: i.e. the three approximation spaces are related as follows, Proof The proof relies on the work of [29].
• Equation (27) shows that K α can be expressed as In case of a straight crack, the local basis is constant.
Thus, the local basis expresses in cartesian coordinates as, Then, • In case of a curved crack, the result of the lemma holds with a discretization of the local basis of [30]. Let us consider a discrete local basis at each node e 1 k , e 2 k , e 3 k , such as, The straight crack proof still holds here, with the additional node index "k".
Remark as a practical consequence of Lemma-1, the error bound of X-FEM should be lower than the ones of V h and of the cutoff [29]. With X-FEM, the optimization process performs on a larger space than with other formulations and reaches a closer infimum to the exact solution of the problem. Nevertheless, X-FEM introduces more d.o.f. than GFEM and cutoff enrichments which increases its condition number.

"Bubble" enrichment class
Now, let us consider the SGFEM space for linear elastic fracture mechanics similar to the one proposed in [25,31]: δI CT is defined here as the transition layer of nodes between crack-tip elements and other elements of the mesh for which K α unknowns are set to zero. This allows to connect the enriched layer with the remaining of the mesh and to avoid blending issues [26].
, which can be shown by contradiction.
• Let us assume • Singular functions cannot be represented by discontinuous polynomials basis, so we necessarily have b i,j,k = 0 • Singular functions cannot be represented by continuous polynomials as well, • Φ i K α introduces at least three higher order polynomials per node, Remark As a practical consequence of Lemma-2 the "bubble" class of "Bubble" enrichment class" section does not belong to the "straightforward" class of "Straightforward" singular enrichment class" section, in which X-FEM is found to be the largest space. Hence, the definition of two classes of enrichment makes perfectly sense.

Numerical behavior of strongly discontinuous and singular approximations
In this section the numerical behavior of the different approximation spaces introduced previously is investigated through a couple of benchmarks. First of all, we assess the asymptotic behavior of strongly discontinuous approximations, through a one dimensional case. Secondly, we consider the convergence study of a two-dimensional problem with a crack and singular approximation.

Numerical analysis of strongly discontinuous approximations with linear elements
Let us consider the basic case of the traction of a one-dimensional bar with a fictitious bimaterial interface. The arbitrary interface is positioned along the abscissa x (ε) (Fig. 8). The problem even though continuous is treated as a discontinuous one, with gluing interface boundary conditions imposed through a Lagrange multiplier. Actually, if the materials were different on each side of the interface, the resulting problem would be equivalent to enforcing a weak discontinuity with a strong discontinuity framework. This artefact is used so as to establish convergence results in energy norm. If it was not the case, the solution obtained would be that of two rigid bodies with prescribed displacements on one of their end. • The problem is the solution of the following differential equation: • With arbitrary Dirichlet boundary conditions: • And with continuity conditions at the interface expressed in terms of the continuity of the derivative of the displacement field (stress continuity at the interface): This continuity condition results from the equivalent Lagrangian form of (30) in which a Lagrange multiplier λ is used to impose a continuous displacement across the interface.
(u, λ) is the saddle point of the following functional: The expected continuous solution satisfying the boundary conditions above, is: The resulting weak form of the problem is: The X-FEM discrete linear space is: Then, the 6 × 6 matrix associated with the discretization of the weak form on the X-FEM space of linear functions (first column of Table 1) With the domain formulation of the second column of Table 1, the discrete matrix is: ⎡ With the shifted formulation of the third column of Table 1, the discrete matrix is: ⎡ We enforce boundary conditions directly on the discrete problem, as equality constraints through Lagrange multipliers. The position of the interface depends only on ε. Nonetheless, the solution does not depend on the ε parameter, so that the relative error does not depend on ε, which makes our analysis relevant. Let us consider the following arbitrary numerical values: α = 10/9 and β = 10.
Then, the solution is computed in standard 64-bit arithmetic, with the linear solver UMFPACK, to reproduce the behavior of direct solvers.
The numerical error has to be close to zero, for any given position of the interface. The error is evaluated here in terms of the H 1 -norm: Although all three formulations represent the same approximation space, they do not have exactly the same numerical behavior (Fig. 9), at least concerning the error in H 1 -norm. The X-FEM error increases steadily while the "shifted" formulation and the "domain" formulation levels of error stay close to machine precision (about 2.2 · 10 −16 ). Moreover, all three enrichment conditionings increase as the distance between the interface and node N 3 decreases (Fig. 10). This means that all three enrichment strategies are very sensitive to the position of the interface. A pre-conditioner is needed for the three enrichments. In the case of X-FEM, we studied Béchet et al. [21] pre-conditioner. For other formulations we used a diagonal pre-conditioner to scale the d.o.f.
The diagonal pre-conditioner allows a scaling of rows and columns through a multiplication with a diagonal matrix to the left and to the right: In linear elasticity, D c is usually defined as, Thus, on Fig. 11, condition numbers are reduced drastically. However, the error still increases in case of X-FEM. The error is clearly not only related to conditioning. Another explanation should be considered. When looking at the X-FEM discrete stiffness matrix, it is noticeable that the unknowns a 2 and b 2 obey to an almost similar equation (as H 2 ≈ − 2 ). The difference between those equations lies in the terms 1 − 2ε and 2 − 2ε. Those terms imply the sum of heterogeneous quantities that leads to a tremendous loss of accuracy on the difference information 2ε. For instance, let us consider the following string of calculations in "double precision" arithmetic (8-octet storage):  The use of dedicated pre-conditioners improves the condition number but has little effect on X-FEM accuracy. According to the curves above, conditioning and accuracy are clearly two separated issues: improving the condition number does not automatically solve accuracy issues. Conditioning and accuracy are both symptoms of the faulty behavior of partition of unity with asymptotically small domains The final result 6.348E-12, is quite different from the exact result 6.3487497084E-12: only four digits remain of the "difference information" between the equations. Hence ill conditioning indicates also a lesser accuracy within the assembled stiffness matrix due to round-off errors and truncated information. Even highly effective pre-conditioners [21] cannot recover the truncated information. It is not surprising that, besides preconditioning, authors have used triple precision arithmetic to recover accuracy [24].
Remark Because the error curves are shattered, the phenomenon might be also probabilistic and related to random round-off errors in interaction with the solver algorithm. A comprehensive study with a wide range of solvers is presented (Figs. 12 and 13). This study shows that the type of solver interferes with the level of error, which makes numerical analysis on errors quite sensitive: • For iterative methods (PCG), the error is generally higher because those solvers are very sensitive to conditioning, • For matrix factorization based solvers (UMF), the error is lower than with iterative methods, • For matrix inversion based solvers, good results can be obtained frequently.

Numerical analysis of strongly discontinuous approximations with quadratic elements
The same one-dimensional case is considered here. The same discretization of the weak form is used with quadratic shape functions. The results are plotted on Fig. 14.
For the three approximations, conditioning worsens very quickly, at about three times the rate of linear elements. More worrying is that the error increases. X-FEM enrichment error worsens more quickly than the others.  12 Overall error analysis of a wide range of solvers for the 1D case. Only X-FEM's enrichment without pre-conditioner is studied here. The whole trend suggests a global increase of the sensitivity of solvers (with respect to epsilon), when the "epsilon" parameter goes to zero Fig. 13 The use of dedicated pre-conditioner stabilizes the behavior of solvers. But the error still increases with the X-FEM method, as explained in "Numerical analysis of strongly discontinuous approximations with linear elements" section On Fig. 15, the same pre-conditioners than in "Numerical behavior of strongly discontinuous approximations" section are considered. The wrong behavior noticed above still applies.
This numerical behavior with quadratic elements is well explained by the coupling between enrichment functions, as pictured on Fig. 16 with domain enrichment [33].
When measure( + ∩ Supp(Φ S ) ) → 0 for a given vertex node, the middle node satisfies likewise measure( + ∩ Supp(Φ M ) ) → 0, on the vertex patch. For such a configuration, the enrichment shape functions χ + Φ S and χ + Φ M become almost homothetic (Fig. 16).  • In case of Lagrange polynomials, we have: Both functions can be related with the following expression, when ξ goes to 1, the difference information between those functions decreases accordingly to (1 − ξ ) 2 . This generates a vanishing sub-space (consisting of polynomials only defined on the small support + ∩ Supp(Φ S ) ), which implies a conditioning slope of at least 3 with quadratic shape functions (Fig. 15). In other words, the condition number κ(K ) is bounded by: Proof Let κ(K ) be the condition number of the stiffness matrix: where 2 is the Euclidian norm 2.
In order to estimate the lower bound of the condition number, we consider the polynomials e 1 and e 2 illustrated on Fig. 17.
Those polynomials belong to the discrete space of [33] (likewise [7] and [26]). They can be used to determine a lower bound of the condition number. Functions e 1 and e 2 are expressed in the discrete space of [33] as: In the discrete space of [33], we have: so that we have e 1 2 = e 2 2 = 1. As 1 2 e T 1 Ke 1 and 1 2 e T 2 Ke 2 represent the elastic energy of these solutions, we have: Thus, The asymptotic behavior of Φ M (ξ ) [at the first order of (1 − ξ )] is: Then, the related d.o.f. become redundant, which leads to ill conditioning. As enrichment strategies are equivalent, the redundant d.o.f. pollutes also the approximation spaces of the other formulations [7,26].
• In case of Bernstein polynomials, we have: when ξ goes to 1, there is no coupling between Bernstein polynomials as observed with Lagrange polynomials. However, the polynomial Φ S cancels out accordingly to (1 − ξ ) 2 , which also leads to the same conditioning slope of three ( Fig. 18) (as shown in the previous section). It is noticeable that [33] has better accuracy with Bernstein polynomials than with Lagrange polynomials (Figs. 18, 19).

Numerical analysis of crack approximations with linear elements: crack opening in mode 1
Given the capabilities of the software we used, only two approximations were tested here: • X-FEM/GFEM vectorial enrichment [30], • X-FEM scalar enrichment [7].
Other approximations are out of scope because, -Gupta et al. SGFEM kind of enrichment [25,31] (let us recall that accordingly to [32] a genuine SGFEM enrichment is defined by the fact that the condition number of the associated stiffness matrix is of the same order-2 [38]-with respect to mesh refinement than the one of FEM) needs far too many additional Heaviside d.o.f. which are difficult to implement and may lead to conditioning issues, The analytical solution in displacement corresponds to the exact mode I limited to the √ r term in Williams series expansion [11] and is given by the first term of (15). Then, Dirichlet boundary conditions are applied on three sides. On the left side, Neumann conditions are preferred. As the left side is cut by the crack, Dirichlet conditions are more difficult to apply given the discontinuity of the displacement field (Fig. 20).
The condition number is estimated by MUMPS solver, so that the results shown here have to be taken only as rough estimates.
On Fig. 21, the condition number of X-FEM almost skyrockets as noticed in [21]. We did not consider the pre-conditioner of [21] here to show how both enrichment strategies, scalar and vectorial, behave without expensive treatment. Even the conditioning slope of the vectorial enrichment is far from optimal. The expected optimal conditioning slope is two [25].
Then, the relative error in energy norm is computed accordingly to the following formula:

Fig. 19
Asymptotic behavior with Bernstein polynomials with dedicated pre-conditioners (same pre-conditioners as with linear elements). As with Lagrange polynomials, dedicated pre-conditioners do not solve accuracy issues  Fig. 22, the rate of convergence in energy norm of the X-FEM scalar enrichment is under-optimal (around 0.91) as noticed in [17]. Optimality is recovered with the vectorial enrichment, the convergence rate being then 0.989. Nonetheless, X-FEM is more accurate than the vectorial enrichment as predicted in Lemma 1.
Example 2 (Inclined crack opening in mode I with linear elements). We introduce a more realistic case with a chosen inclination of the crack at 44.9 • , in order to test the conditioning with both Heaviside and singular enrichments. The analytical solution is still the one of a mode I problem limited to the √ r term in Williams series expansion given by the first term of (15) (Fig. 23). The inclined crack analytical solution is expressed through a rotation of 44.9 • of the horizontal crack solution. The strain and stress tensors are rotated as well of the same angle. Then the same mix of Neumann and Diriclet boundary conditions  are applied (Fig. 24), similarly to the case of the horizontal crack. We consider the same regular meshes as above which are not oriented accordingly to the crack surface here. The pre-conditioner of [21] is not applied here as in the previous case.
On Fig. 25, the conditioning of the X-FEM vectorial enrichment becomes very high and reaches a steady value around 10 12 . This behavior may be explained by the ill conditioning of the Heaviside enrichment. As the "Heaviside information" becomes very unbalanced, its ill conditioning overcomes the regular increase of conditioning expected with the vectorial enrichment (see Fig. 21) as explained in "Strong discontinuity approximation conditioning" section. Nonetheless, both enrichments have optimal convergence accuracy (see Fig. 26).

Numerical analysis of crack approximation with quadratic elements: crack opening in mode 1
Only the linear vectorial enrichment approximation [30] is working with a reasonable condition number. We recall that other approximations are not tested because, • X-FEM conditioning skyrockets with quadratic elements, even with a small enrichment zone, • Gupta et al. SGFEM kind of enrichment [25,31] needs far too many additional Heaviside d.o.f., so that it is difficult to implement and may lead to ill conditioning, • Cut-off with global d.o.f. assembling requires a "macro" element, which is not straightforward to implement in finite element software. Moreover the method does not extend properly in 3D [29].
Even if the linear vectorial enrichment converges at an optimal rate (in energy norm) with quadratic elements, the condition number is still very high (see Fig. 27). We did not use the pre-conditioner of [21], in order to test the conditioning of the vectorial enrichment. In "Singular approximation at a crack tip" section, a new enrichment strategy with better intrinsic numerical behavior, will be exposed.
Remark a dedicated integration scheme is needed at the tip of the crack, to get an optimal rate of convergence in energy norm with quadratic elements [17]. Here, we used a Gauss-Radau integration rule of order 20 [39]. As the aim of the paper is not about singular integration, we did not try to optimize the procedure. Optimization of integration schemes has been well studied in [30,40].

Improving the design of quadratic approximation spaces to deal with previous numerical issues
In this section, basic improvements are suggested to deal with numerical issues stated in previous sections: -For strongly discontinuous approximations, we suggest a correction to the approximation spaces for quadratic elements. Note that as denoted by [17] optimal con- vergence rates can only be achieved if the discontinuous approximation space is of the same order than the one of the continuous space, that is to say quadratic. The correction exposed here on the quadratic form is slightly different from the ones discussed in "Strong discontinuity approximation conditioning" section (fit-to-vertex and elimination of d.o.f.). -For singular enrichment, a reshape of the approximation space is considered. The new approximation sums up benefits of work on strong approximation spaces with the significant improvement associated to the "bubble" space as discussed in the next section.

Partition of unity alteration
As discussed in "Numerical analysis of strongly discontinuous approximations with quadratic elements" section, the use of vertex node and middle node shape functions (resp. χ + Φ S and χ + Φ M in Fig. 16) leads to an incorrect asymptotic behavior of strongly discontinuous formulations, both for Bernstein and Lagrange polynomials. When getting rid of the middle node d.o.f. χ + Φ M for ε around 10 −3 , the condition number decreases sharply (Fig. 28). This threshold value is close to estimates of [18,24], but we stress on the fact that only the middle d.o.f. is removed and that the position of the interface is not shifted. This alteration process allows the analysis to proceed beyond ε = 10 −12 as with linear elements.
Let us precise the differences between the removal of the middle node and straightforward criteria from "Strong discontinuity approximation conditioning" section: • The fit-to-vertex and volume criterion remove both the vertex node and the middle node d.o.f. (χ + Φ S and χ + Φ M ). Then, with the fit-to-vertex, the interface position is switched to the closest node.
• With the new strategy, only the minimal information is removed from the approximation (a single d.o.f., the nearest to the interface) with a tuned threshold. The interface is not switched.
However, the alteration of the partition of unity has a noticeable impact on enrichment [33] with Lagrange polynomials (Fig. 28): the H1-error increases sharply after the alteration, instead of decreasing as observed with other formulations. This behavior might be explained by the partition of unity and the nature of the solution represented here: • The exact solution is continuous, • In [7,26] discontinuous functions are added to continuous functions, so that the removal of the discontinuous middle node d.o.f. does not alter the partition of unity of continuous space functions [9] and consequently the representation of continuous functions.
This is in contrast with [33], as the removal of the enriched middle node d.o.f. prevents the representation of continuous functions on some parts of the domain where the partition of unity is broken.
A multi-axial loading is applied. Given the elastic behavior, the stress is homogeneous within the block at a value of 10 MPa, for any given position of the interface. The interface equation is parameterized with δ as: The error is analyzed with the energy norm: As the Hansbo et al. [33] enrichment performs as well as the "shifted" enrichment, only the "shifted" enrichment is tested here. The strategies discussed in 1D are extrapolated in 3D (an example is given Fig. 29): • Preconditioning strategies are unchanged (Béchet et al. pre-conditioner [21] with X-FEM and diagonal scaling with shifted enrichment), • The elimination threshold for the middle nodes, is fixed at a ratio of relative distance to the vertex of 10 −3 along the edge (see on the 1D example illustrated by Fig. 16), or, could be expressed with a volume ratio criterion [24]: On (Fig. 30), the shifted formulation is more accurate than the X-FEM one, as observed in the 1D case. The difference in accuracy reaches a 10 10 factor: accuracy issues observed in 1D with X-FEM are heightened in 3D. Furthermore, the volume ratio criterion on the middle node is very satisfactory. The results improve as soon as the elimination process is activated around a node to interface distance (here δ/ √ 3) of 10 −3 m.

Synthesis
Although, all three strongly discontinuous formulations considered in the paper represent the same approximation space, they do not have the same numerical performance. From the results above, domain and shifted enrichments are less sensitive to the position of the interface than X-FEM. With quadratic elements, the three formulations need a special care regarding the asymptotic behavior of the shape functions as shown in "Numerical

Fig. 30
Behavior of strongly discontinuous enrichments in 3D with the corrections studied previously analysis of strongly discontinuous approximations with quadratic elements" section. The X-FEM jump formulation is far less accurate in 3D. From a practical point of view, the "shifted" formulation seems to offer the best compromise, because it provides: • A simpler pre-conditioner than X-FEM, which decreases the computational cost, • Less discontinuous d.o.f. than [33], which eases the implementation and improves the accuracy of the formulation in case of elimination of d.o.f. with Lagrange polynomials, • A good accuracy in 3D with quadratic elements.

New "bubble" approximation space
We introduce a simpler space than the one of Gupta et al. [25,31], to reduce its too many additional Heaviside d.o.f. As a matter of fact Gupta et al. approach requires up to 12 Heaviside d.o.f. per node instead of 3 so that it is difficult to implement and may lead to ill conditioning.
In the elements where K α is discontinuous, a technique similar to the ghost node interpolation used in [41] is preferred. The principle consists in using two continuous interpolations by prolongation of the one on + and the one on − to the whole domain (an example is given Fig. 31). Hence, two continuous extensions of K α (r, θ) are introduced: In the elements where K α has a singular point, the interpolation of K α is smoothened, so that subsequent interpolation errors of the singular function K α do not pollute the accuracy of the displacement field as observed in [25].

Fig. 31
Comparison between double node interpolation and ghost node interpolation for an arbitrary function. The ghost node method uses the same element to interpolate both extrapolated branches. In contrast, the double node refines the mesh around the discontinuity. Thus, the double node is twice more accurate than the ghost node, but it requires a meshing of the interface (a strategy not relevant in the framework of X-FEM/GFEM) Let n be the order of the interpolation used for the crack tip degrees of freedom. I CT is the corresponding set of nodes enriched by the vectorial asymptotic functions i.e. I CT,1 corresponds to vertex nodes and I CT,2 contains both vertex nodes and middle nodes. If we consider a quadratic interpolation for the crack tip degrees of freedom, we have n = 2 and I CT,2 = I CT . Let I T be the associated subset of I CT with singular points i.e. the set of nodes such that a crack-tip belongs to their support. The nodes of I T are not taken into account in the interpolation of K α , so that the interpolation of K α is smooth when the polar coordinate r goes to zero. Let m be the order of the interpolation of K α , which can be different from n. The new interpolation operator reads: where k,m are the shape functions corresponding to order m interpolation and I CT,m is the subset of I CT corresponding to order m interpolation, i.e. I CT,1 is the restriction of I CT to vertex nodes and I CT,2 contains both vertex nodes and middle nodes of I CT .
Since we consider quadratic approximation in this section, we have I CT,2 = I CT . The set I CT,m ∪ δI CT,m contains all the nodes lying in the support of a node of I CT . This interpolation should be interpreted, for an element which is split by the interface (i.e. an element that does not contain a crack-tip), as follows • Assuming the evaluation point (gauss point) is located on + , • Assuming the evaluation point (gauss point) is located on − , The new "bubble" approximation space, based on the new interpolation operator, is: where Φ i are the shape functions corresponding to the order of the elements, i.e. p = 2 in this section and n refers to the order of the shape functions i,n used for the vectorial enrichment discretization.
Remark • Let us stress on the fact that the ghost node interpolation does not introduce additional nodes. The two branches of discontinuous functions are interpolated over the element through extrapolation of branches (Fig. 31); • The difference between the original SGFEM and the "new" bubble space V h bub,m,n lies in the behavior of the interpolation in the bandwidth of elements where the enrichment function becomes unsmooth or discontinuous. The original SGFEM uses the same interpolation operator whether the enrichment function is smooth or discontinuous. The new "bubble" space replaces the enrichment function with smooth continuous extensions in the bandwidth of elements where the enrichment function becomes discontinuous. Thanks to the ghost node interpolation, both smooth extensions are interpolated without additional node (Fig. 31).

Comparison with previous approximations
The mode I problem of "Numerical behavior of singular approximations at the crack tip" section is considered here, with a similar enrichment radius of r = 0.1. In this section, we focus on the convergence analysis carried out with quadratic elements.
On Fig. 32, the "bubble" enrichment has a better condition number than the linear vectorial enrichment. On Fig. 33, all enrichment strategies have optimal convergence rates in energy norm. The bubble transformation of the singular function preserves the optimality of the vectorial formulation.
Nonetheless, bubble spaces V h bub,m = 2,n = 1 and V h bub, m = 2, n = 2 are about five times more accurate than the bubble space V h bub, m = 1, n = 1 . The bubble space V h bub, m = 1, n = 1 has the same accuracy than X-FEM vectorial enrichment. Clearly, the bubble spaces with quadratic interpolation of the singular function outperform X-FEM and the linear bubble space. As a matter of fact, among the bubble spaces with quadratic interpolation of singular function denoted V h bub, m = 2, n = 1 and V h bub, m = 2,n = 2 , we recommend to choose n = 1 (linear enrichment for the singular part) but m = 2 (quadratic interpolation of the singular function) to have the best accuracy and lowest condition number with an order with respect to mesh refinement around 3.25 for an optimal expected value of 2 [31,32,38] (see Fig. 32).

Remark • V h
bub, m = 1, n = 2 is not very interesting because it has a poor numerical behavior. The shift with a lower order of interpolation has not been considered even in the original method [25]. In addition, an accuracy study of stress intensity factors is also shown with convergence orders with respect to mesh refinement close to the theoretical value of 2 (Fig. 36).  are estimated with a domain integral method (G-theta method [42,43]). The "bubble" formulation with quadratic elements reaches the limit of numerical accuracy of the estimator used so that a convergence analysis cannot be performed here. For formulations with linear elements, the convergence rates are close to the optimal rate of convergence with mesh refinement, around 2

Extension to weak discontinuities
Similarly to the analyses performed for strong discontinuities and singular enrichment, this section will describe several enrichments for weak discontinuities found in the literature. Those are used in case of bi-materials. We will show that getting a result with correct accuracy does not depend only on conditioning issues but also on the quality of the approximation space chosen to represent continuous fields with discontinuous gradients. These two separate issues are often mixed in the literature as it will be shown later in this section.

Enrichment functions for weak discontinuities
The general enrichment in case of weak discontinuity takes the following form: with the same notations than in "Strong discontinuity representation" section. The first sum in the right hand part of the previous equation corresponds to the standard finite element approximation while the second one stands for the enriched term. Different choices are present in the literature for the function F (x) and the ones in [14,27,28] are possible ( Table 3). Some of them are related: as a matter of fact, the choice by [28] is a shifted version [26] of [14] and [27] is a bubble transformation of [14].
We analyze in the following the influence of these choices on the overall behavior of the solution with respect to the satisfaction of simple patch tests convergence and accuracy properties.
Remark Sukumar's and Belytschko's enrichment functions yield almost the same approximation space. It can be shown straightforwardly through a change of variables as in "Strong discontinuity approximation conditioning" section. Let a j and b j be the classical and enriched degrees of freedom for Belytchko's enrichment and a j and b j refer to the classical and enriched degrees of freedom for Sukumar's enrichment. We have for j ∈ I H Thus, in the next sections both approximations may display similar numerical properties. Note also that the zero of the level set term i∈I H lsn i Φ i (x) which appears in [14,27,28] represents a discretization of the interface. In the following we took the same orders of interpolation for the functions F (x) and the level sets, which seemed to be quite natural.

Numerical analysis of weakly discontinuous approximation
In order to see if conditioning issues are also relevant for weak discontinuities we performed a numerical analysis similar to the one of "Numerical behavior of strongly discontinuous approximations" section, relying on a one dimensional bi-material problem.
We will show that even if these conditioning issues can be prevented by a change of formulation, convergence issues can still be observed, due to the fact that simple patch tests cannot be represented correctly for linear or quadratic elements.

Analytical solution for the one dimensional bi-material problem
To capture piecewise linear displacements as in "Numerical analysis of strongly discontinuous approximations with linear elements" section, we simulate a one-dimensional bi-material problem of length L in a two dimensional domain since our numerical software is 2D and 3D and as was done in [14]. To ease calculations, we reduce it to the one-dimensional bi-material problem.
• For the one dimensional bi-material bar, the solution verifies the following differential equation , where x 0 is the location of the interface, • With arbitrary Dirichlet boundary conditions we choose to impose on the central element nodes for the sake of simplicity, • Material properties are given by: , where x 0 is the location of the interface, • Assuming small strains, the expected solution satisfying the boundary conditions above is:

Numerical analysis of weakly discontinuous approximations with linear elements
In the numerical analysis of "Numerical analysis of weakly discontinuous approximations with linear elements" section, to ease the calculations, we consider = [− 1, 2] with the following Dirichlet boundary conditions u (-1) = 0, u (2) =ū so that L = 3 m. The interface is = {x 0 }.
The weak form of the problem becomes: Find: The X-FEM discrete linear space is: Then, the 6 × 6 matrix associated with the discretization of the weak form on the X-FEM space of linear function with Sukumar et al. formulation [14] is: With the formulation of Belytschko et al. [28] it becomes: Finally, with Moës et al. formulation [27], the discrete matrix is: where: For all cases, we have taken ε = 1 − x 0 , where x 0 denotes the interface position. On considering the 4 × 4 sub-matrix that is free from boundary conditions, the diagonal pre-conditioner can easily scale the formulation of Moës et al. [27] as entire rows (symmetrically columns) vanishes when ε = 0 or ε = 1 (see Fig. 37). For the enrichments of Sukumar et al. [14] and Belytschko et al. [28], no conditioning issues can be observed, in case of linear approximation. However, we will see in the next paragraph that these latest two formulations are far from optimal with respect to convergence issues, being unable to represent simple patch tests.

Ability of the enrichment to catch piecewise linear and quadratic solution
First of all, we would like to test that the proposed enrichments (see Table 4) are able to capture piecewise linear and quadratic solutions, so that the enrichment function for the material interface problem preserves the equivalence between the FEM and X-FEM discretized spaces. For that purpose, we will consider the previous problem of the onedimensional bi-material bar for the linear case [14]. Then to capture piecewise quadratic displacements we adapt the bi-material volumic load test of [27] using a constant volumic load, instead of a quadratic one.
To simplify the calculations, the bar is subdivided into three elements (see Fig. 8) of equal length L, in order to have cut and uncut elements. The central element which carries the interface is such that 0 ≤ x ≤ L. Dirichlet boundary conditions are applied on this central element such that u(0) = 0, u (L) = u. The linear shape functions we use in the central element are of the form 1 = 1 − x L , 2 = x L , while the quadratic ones are taken as: Fig. 37 Conditioning of X-FEM formulation with quadratic elements for weak discontinuities representation. With X-FEM's ridge enrichment functions no conditioning issues are observed through the use of a diagonal pre-conditioner Unknowns that we search to recover the linear and quadratic patch test solutions are the degrees of freedom a 1 , a 2 , b 1 , b 2 and a 1 , a 2 , a 3 , b 1 , b 2 , b 3 , in (58) respectively. For Sukumar and Belytschko formulations, the X-FEM blending elements which are adjacent to those completely cut by the interface have a quadratic interpolation [even with linear elements, due to the interpolation (58) and the expression of F (x) in Table 4] or cubic interpolation [with quadratic elements and the expressions (58) and F (x)] with just one active b i degree of freedom (see Fig. 8, for the left and right blending elements). Since the patch test is linear, these degrees of freedom must be equal to zero which cancels out all discontinuity. Hence linear patch tests cannot be satisfied and are just approximated by these formulations, which will affect the energy norm as noticed by [14]. Note that Moës formulation with the same order of interpolation for F (x) and the shape functions Φ(x) is able to catch the linear patch test solution, and that Moës formulation with the linear interpolation of F (x) is also able to catch the solution for quadratic or linear shape functions. Table 5 gives a summary of the principal results obtained for the linear patch test. Table 6 represents the relative error in energy norm for this bi-material bar under traction for different positions of the interface. One can see that the enrichments of Sukumar and Belytschko, even when quadratic (see Table 7), do not yield proper results due to the transition of F (x) in blending elements [26], while the problem is avoided with Moës enrichment.
Moreover, as we cannot recover the exact linear solution for this patch test, an order of convergence of 0.5 in energy norm with respect to mesh refinement is obtained for Sukumar and Belytschko enrichments (Fig. 38), which is consistent with the previous work of Moës et al. [44] extended to X-FEM with an enrichment that does not allow to catch the discontinuity. Analytical solution for the volumic load problem Consider a plate of side L with two materials E 1 and E 2 with ν 1 = ν 2 = 0, the interface being localized at x = x 0 in the plate, as in Fig. 39. Plain strain conditions are assumed. The plate is embedded on the left side and a variable volumic force is imposed on the plate along the x direction, given by f = f . The case is computed in two dimensions since our numerical software works only in two and three dimensions, using a similar approach to the one in [14]. Table 3 for the one dimensional bi-material bar problem Sukumar [14] M o e s [ 27] Belytschko [28] Linear Table 5 Coefficients used to capture the solution with the finite element approximation for the linear patch test Sukumar [14] M o ë s [ 27] Belytchko [28] Linear φ(x) and

Table 4 Expressions of F(x) for weakly discontinuous enrichments with linear and quadratic approximations resulting from
The analytical solution cannot be represented. Incompatible with blending elements.
The analytical solution cannot be represented. Incompatible with blending elements.
Quadratic φ(x) and The analytical solution cannot be represented. Incompatible with blending elements.
• if The analytical solution cannot be represented. Incompatible with blending elements.
• if The analytical solution cannot be represented. Incompatible with blending elements.
The analytical solution cannot be represented. Incompatible with blending elements.  The exact displacements in the two materials are given by: No linear element with linear F (x) is able to capture this quadratic solution directly. For Sukumar and Belytschko formulations with quadratic elements, the xfem blending elements which are adjacent to those completely cut by the interface give a cubic interpolation (linear F (x)) or order four interpolation (quadratic F (x)) with just one b i degree of freedom (see Fig. 8, for the left and right elements). Since the patch test is quadratic, these degrees of freedom must be equal to zero which cancels out all discontinuity. Hence quadratic patch tests cannot be satisfied and are just approximated by these formulations, which will affect the energy norm as noticed by [14]. Note that Moës formulation with quadratic shape functions Φ(x) and linear F (x) is able to catch the quadratic patch test solution, and that Moës formulation with quadratic shape functions Φ(x) and quadratic F (x) is not able to catch it. Table 7 gives a summary of the principal results obtained for the quadratic patch test. Figure 40 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the quadratic patch test, using linear shape functions Φ(x) and linear F (x) and quadratic shape functions Φ(x) and quadratic F (x). Using linear elements we obtain the optimal convergence rates of 1. Using quadratic elements we obtain 1.5. In this later case, a convergence rate of 2 would be expected since the analytical solution does not lie in the finite element space (cf. Table 8). We will not prove formally that the obtained convergence rates are in fact the expected ones, nevertheless we can give some arguments that those rates are indeed the expected ones.
In the quadratic patch test case, the first and second derivatives of the analytical solution are discontinuous. The finite element space with enrichment we use with quadratic F (x) and quadratic shape functions Φ(x) does not allow to capture this solution. However it allows to capture the linear patch test solution with a discontinuous first derivative. We can deduce from this analysis that its rate of convergence on energy will be greater than 1 but lower than 2. Then the loss of 0.5 in convergence rate observed numerically for quadratic elements, with a value of 1.5 observed numerically instead of an expected 2, is coherent with the analysis for linear elements (0.5 instead of 1 when elements do not catch the discontinuity on the first derivative) performed by Moës et al. [44].
From this analysis, we recommend to use Moës formulation with quadratic elements and linear F (x).

Convergence analysis for weak discontinuities
In this paragraph we will investigate the convergence of our X-FEM formulations dedicated to weak discontinuities when the solution is not in the approximation space. The analytical solution cannot be represented. Incompatible with blending elements. Fig. 40 Convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the quadratic patch test, using linear shape functions and a linear enrichment function (linear, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic)

Volumic load problem with straight interface
The test we propose in the following was already investigated by Moës et al. [18]. Let us consider again the squared plate of Fig. 39 with side L = 2 m separated in its middle by two materials E 1 = 1 MPa and E 2 = 10 MPa and with ν 1 = ν 2 = 0. This time the variable volumic force is f = f x L 2 imposed on the plate along the x direction.
The exact displacements in the two materials are: We considered an unstructured mesh with 5, 10, 20 and 40 triangular elements on each side. Figure 41 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the volumic load problem, using linear shape functions Φ(x) and linear F (x), quadratic shape functions Φ(x) and linear F (x) and quadratic shape functions Φ(x) and quadratic F (x). Using linear elements or quadratic elements and linear enrichment function we obtain the optimal convergence rates. Using quadratic elements and quadratic enrichment function we obtain 1.5. Note that the second derivative of the analytical displacement is discontinuous in this case too.

Circular inclusion with imposed displacement
To investigate the behavior of weak discontinuity enrichments in presence of curved interfaces, Moës et al. [18] suggested the following test case. A circular inclusion 1 of radius a = 0.4 m is placed in the center of a material 2 . A linear displacement is imposed on the boundary 2 (r = b = 2 m). The same material parameters are chosen: E 1 = 1 MPa, ν 1 = 0.25 and E 2 = 10 MPa, ν 2 = 0.3. This problem has an analytical solution provided in [18]. The analytical displacement reads: Fig. 41 Convergence rates on the relative error in terms of the energy for the volumic load problem. The enrichment strategies are: linear shape functions and linear enrichment function (linear, linear), quadratic shape functions and a linear enrichment function (quadratic, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic)

Fig. 42
Convergence rates on the relative errors with respect to mesh refinement for the circular inclusion problem. The enrichment strategies are: linear shape functions and linear enrichment function (weak, linear, linear), quadratic shape functions and linear enrichment (weak, quadratic, linear) and quadratic shape functions and quadratic enrichment function (weak, quadratic, quadratic). Two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface have also been investigated: one with quadratic shape functions and a linear approximation of the level set (strong, quadratic, linear) and one with quadratic shape functions and a quadratic approximation of the level set (strong, quadratic, quadratic) with: The domain of computation is a square 2 of side L = 2 m. We considered an unstructured mesh with 5, 10, 20, 40, 80 and 160 triangular elements by side. Figure 42 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for circular inclusion, using linear shape functions Φ(x) and linear F (x), quadratic shape functions Φ(x) and linear F (x) and quadratic shape functions Φ(x) and quadratic F (x). Using linear elements we obtain the optimal convergence rate. Using quadratic elements we obtain 1.5. In the case of quadratic F (x), this behavior arises from the discontinuity of the second derivative of the analytical displacement. The case of linear F (x) is slightly different, due to its link with the geometrical approximation of the level sets and more particularly of the interface established at the end of "Enrichment functions for weak discontinuities" section. Ferté et al. [45,46] showed for strong discontinuities that the approximation of the level set impacts the expected theoretical convergence rates, in the case of quadratic elements: 2 is expected if a quadratic approximation of the level set is used (equivalently quadratic F (x)), while 1.5 is expected if a linear representation of the level set is used (equivalently linear F (x)). To illustrate this behavior, we also considered two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface: one with quadratic shape functions and a linear approximation of the level set and one with quadratic shape functions and a quadratic approximation of the level set. The convergence rates obtained for these two models are also shown on Fig. 42.
Finally for weak discontinuities, optimal rates of convergence in energy norm are difficult to obtain for quadratic elements. If linear approximations of the level sets (equivalently linear F (x)) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the geometry. If quadratic approximations of the level sets (equivalently quadratic F (x)) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the approximation spaces. It appears that both convergence rates end up being the same, with a loss of 0.5 in energy norm with respect to optimal orders of convergence of 2. Optimal orders of convergence are recovered when using a strong discontinuity approximation space with lagrangian multipliers to represent weak discontinuities.

Conclusion
This paper outlines numerical issues around popular approximation spaces to model weak discontinuities or strong discontinuities and cracks, particularly with quadratic elements. We analyzed those issues on a few benchmarks. Those benchmarks revealed that popular enrichment strategies do not behave well with asymptotic configurations (when the discontinuity gets close to a node of the mesh, typically for a distance below 10% of the length of the cut edge) for strongly discontinuous enrichments and that the singular enrichment could be improved with respect to conditioning and accuracy.
Instead of working on complex and expensive strategies to solve those issues for strongly discontinuous enrichments, we preferred to apply a slight reshape of the approximation spaces, as it is in our sense, the most effective approach to deal with those issues for industrial applications. In fact, the new approximation spaces get through the proposed benchmarks, with significant improvement of the numerical behavior and accuracy of the formulations studied.
Nonetheless, more theoretical work is needed to understand why very unsmooth "bubble" functions tend to give better results than "straightforward" singular functions. At least, our results confirmed the numerical properties of "bubble" spaces denoted in the literature [25,31,32]. Moreover, the integration procedure needs to be improved, particularly in 3D with quadratic elements, to further the analysis on complex crack surface and crack front geometries.
At last, in case of weakly discontinuous enrichments, it was shown that most approximation spaces of the literature did not exhibit conditioning issues. However, for quadratic elements, orders of convergence in energy norm were 0.5 lower than those obtained with an equivalent strongly discontinuous enrichment, so that we recommend using the latest approach even in the case of weak discontinuities, associated to the strategy exposed in the present work to avoid conditioning issues.