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

M. Ndeffo, P. Massin, Nicolas Moes, A. Martin, S. Gopalakrishnan

► To cite this version:

M. Ndeffo, P. Massin, Nicolas Moes, A. Martin, S. Gopalakrishnan. On the construction of approximation space to model discontinuities and cracks with linear and quadratic extended finite elements. Advanced Modeling and Simulation in Engineering Sciences, 2017, 4 (1), pp.6. 10.1186/s40323-017-0090-3. [https://hal.science/hal-03838822](https://hal.science/hal-03838822)

HAL Id: hal-03838822

https://hal.science/hal-03838822

Submitted on 30 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_2.png" /></figure>

CrossMark

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

M. Ndeffo<fnref n="1" />, P. Massin<fnref n="1" />, N. Moës², A. Martin<fnref n="1" />* and S. Gopalakrishnan³

## Abstract

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.

**Keywords:** Fracture, Quadratic elements, Condition number, X-FEM, GFEM, SGFEM

## 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,

<figure><img src="image_3.png" /></figure>

<footnote>*Correspondence:
alexandre-externe.martin@edf.fr</footnote>

<footnote n="1">¹IMSIA, UMR
EDF/CNRS/CEA/ENSTA ParisTech
9219, Université Paris Saclay, 828
Boulevard des Maréchaux, 91762
Palaiseau Cedex, France
Full list of author information is
available at the end of the article</footnote>
<figure><img src="image_2.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_3.png" /></figure>

CrossMark

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

M. Ndeffo<fnref n="1" />, P. Massin<fnref n="1" />, N. Moës², A. Martin<fnref n="1" />* and S. Gopalakrishnan³

## Abstract

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.

**Keywords:** Fracture, Quadratic elements, Condition number, X-FEM, GFEM, SGFEM

## 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,

<figure><img src="image_4.png" /></figure>

<footnote>*Correspondence:
alexandre-externe.martin@edf.fr</footnote>

<footnote n="1">¹IMSIA, UMR
EDF/CNRS/CEA/ENSTA ParisTech
9219, Université Paris Saclay, 828
Boulevard des Maréchaux, 91762
Palaiseau Cedex, France
Full list of author information is
available at the end of the article</footnote>
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–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–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–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-
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–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–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–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.

<figure><img src="image_4.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**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_1IP$ goes to zero or when $N_2IP$ goes to zero</figcaption></figure>
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.

<figure><img src="image_5.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**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_1IP$ goes to zero or when $N_2IP$ goes to zero</figcaption></figure>
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:

$$ measure (\{\Omega_+ \cap \text{Supp}(\Phi_i)\}) \ll measure (\{\Omega_- \cap \text{Supp}(\Phi_i)\}) \quad (1) $$

– or:

$$ measure (\{\Omega_- \cap \text{Supp}(\Phi_i)\}) \ll measure (\{\Omega_+ \cap \text{Supp}(\Phi_i)\}) \quad (2) $$

where, $\text{Supp}(\Phi_i)$ refers to the support of the shape function and $\Omega_+$ and $\Omega_-$ are the domains defined on Fig. 3.

Practically, this means that node numbered $i$ is either on the side $\Omega_-$ or on the side $\Omega_+$ 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-

<figure><img src="image_6.png" /></figure>
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:

$$ measure (\{\Omega_+ \cap \text{Supp}(\Phi_i)\}) \ll measure (\{\Omega_- \cap \text{Supp}(\Phi_i)\}) \quad (1) $$

– or:

$$ measure (\{\Omega_- \cap \text{Supp}(\Phi_i)\}) \ll measure (\{\Omega_+ \cap \text{Supp}(\Phi_i)\}) \quad (2) $$

where, $\text{Supp}(\Phi_i)$ refers to the support of the shape function and $\Omega_+$ and $\Omega_-$ are the domains defined on Fig. 3.

Practically, this means that node numbered $i$ is either on the side $\Omega_-$ or on the side $\Omega_+$ 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-

<figure><img src="image_7.png" /></figure>
vertex. Along each edge (see example on Fig. 2), the balance of “Heaviside information” is weighted through the comparison of lengths [N₁IP] and [N₂IP] with respect to [N₁N₂]. 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₁, we have if lsn refers to the normal level set:

$$
measure(\{\Omega_+ \cap \text{Supp}(\Phi_1)\}) \propto \text{length}([N_1 \text{IP}]) \approx |\text{lsn}(N_1)| \quad (3)
$$

$$
measure(\{\Omega_{-} \cap \text{Supp}(\Phi_{1})\}) \propto \text{length}([N_{2}\text{IP}]) \approx |\text{lsn}(N_{2})| \quad (4)
$$

where N₁ and N₂ 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₁. A similar approach consists in changing directly the value of the level-set of the N₂ 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⁻² and 10⁻³ [18,24].

A simple criterion expresses as:

• Node N₁ is eliminated or the level set is moved to N₂ by resetting lsn (N₂) to zero if,

$$
\frac{|\mathrm{lsn}(N_2)|}{|\mathrm{lsn}(N_1)| + |\mathrm{lsn}(N_2)|} < 10^{-2} \qquad (5)
$$

• Node N₂ is eliminated or the level set is moved to N₁ by resetting lsn (N₁) to zero if,

$$
\frac{|\mathrm{lsn}(N_1)|}{|\mathrm{lsn}(N_1)| + |\mathrm{lsn}(N_2)|} < 10^{-2} \qquad (6)
$$

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 pre-
conditioner [21], dedicated to orthogonalize shape functions and enrichment functions,
with a Cholesky factorization.

Let's assume K is the stiffness matrix. Kᵢ,ⱼ is the term of the matrix associated with shape function Φᵢ and the enriched d.o.f. corresponding to the enrichment function FₓΦᵢ. As K is bilinear, symmetric and positive definite there is an inner product ⟨·, ·⟩ₖ such as:

$$
K_{i,x} = \langle \Phi_i, F_x \Phi_i \rangle_K \tag{7}
$$

Béchet et al. [21] pre-conditioner orthogonalizes Φᵢ and FₓΦᵢ (in the sense of ⟨·, ·⟩ₖ),
according to the following procedure:

$$
\begin{align}
K &\rightarrow \tilde{K} = P_C^T K P_C \nonumber \\
K u &= f \rightarrow \begin{cases} \tilde{K} \tilde{u} = P_C^T f \\ u = P_C \tilde{u} \end{cases} \tag{8}
\end{align}
$$

and,

$$
\tilde{K}_{i,x} = \langle \Phi_i, F_x \Phi_i \rangle_{\tilde{K}} = \langle P_C \Phi_i, P_C F_x \Phi_i \rangle_K = 0 \quad (9)
$$
vertex. Along each edge (see example on Fig. 2), the balance of “Heaviside information” is weighted through the comparison of lengths [N₁IP] and [N₂IP] with respect to [N₁N₂]. 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₁, we have if lsn refers to the normal level set:

$$
measure(\{\Omega_{+} \cap \text{Supp}(\Phi_{1})\}) \propto \text{length}([N_{1}IP]) \approx |\text{lsn}(N_{1})| \quad (3)
$$

$$
measure(\{\Omega_{-} \cap \text{Supp}(\Phi_{1})\}) \propto \text{length}([N_{2}IP]) \approx |\text{lsn}(N_{2})| \quad (4)
$$

where N₁ and N₂ 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₁. A similar approach consists in changing directly the value of the level-set of the N₂ 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⁻² and 10⁻³ [18,24].

A simple criterion expresses as:

• Node N₁ is eliminated or the level set is moved to N₂ by resetting lsn (N₂) to zero if,

$$
\frac{|\operatorname{lsn}(N_2)|}{|\operatorname{lsn}(N_1)| + |\operatorname{lsn}(N_2)|} < 10^{-2} \qquad (5)
$$

• Node N₂ is eliminated or the level set is moved to N₁ by resetting lsn (N₁) to zero if,

$$
\frac{|\mathrm{lsn}(N_1)|}{|\mathrm{lsn}(N_1)| + |\mathrm{lsn}(N_2)|} < 10^{-2} \qquad (6)
$$

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 pre-
conditioner [21], dedicated to orthogonalize shape functions and enrichment functions,
with a Cholesky factorization.

Let's assume K is the stiffness matrix. Kᵢₓ is the term of the matrix associated with shape function Φᵢ and the enriched d.o.f. corresponding to the enrichment function FₓΦᵢ. As K is bilinear, symmetric and positive definite there is an inner product ⟨·, ·⟩ₖ such as:

$$
K_{i,x} = \langle \Phi_i, F_x \Phi_i \rangle_K \tag{7}
$$

Béchet et al. [21] pre-conditioner orthogonalizes Φᵢ and FₓΦᵢ (in the sense of ⟨·, ·⟩ₖ),
according to the following procedure:

$$
\begin{align}
K &\rightarrow \tilde{K} = P_C^T K P_C \nonumber \\
K u &= f \rightarrow \begin{cases} \tilde{K} \tilde{u} = P_C^T f \\ u = P_C \tilde{u} \end{cases} \tag{8}
\end{align}
$$

and,

$$
\tilde{K}_{i,x} = \langle \Phi_i, F_x \Phi_i \rangle_{\tilde{K}} = \langle P_C \Phi_i, P_C F_x \Phi_i \rangle_K = 0 \quad (9)
$$
If the degrees of freedom for each node are adjacent (ranked contiguously), the shape of $P_C$ is:

$$P_C = \begin{bmatrix} P_C^1 & 0 & 0 & \cdots & \cdots \\ 0 & P_C^2 & 0 & \cdots & 0 \\ 0 & 0 & P_C^3 & \cdots & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \cdots & 0 & \cdots & \cdots & P_C^i \\ \cdots & \cdots & \cdots & \cdots & \ddots \end{bmatrix} \quad (10)$$

The blocks of $P_C$ express as follows:

$$P_C^i = \begin{cases} S_i^{-1} & \text{if the node is enriched} \\ I_d & \text{if not} \end{cases} \quad (11)$$

where the inverse terms are computed from Cholesky's factorization of local blocks of K:

$$K_{loc}^i = S_i^T S_i \quad (12)$$

$$K_{loc}^i = \begin{bmatrix} \langle \Phi_i, \Phi_i \rangle_K & \langle \Phi_i, H\Phi_i \rangle_K & \langle \Phi_i, F_x\Phi_i \rangle_K \\ \langle \Phi_i, H\Phi_i \rangle_K & \langle H\Phi_i, H\Phi_i \rangle_K & \langle H\Phi_i, F_x\Phi_i \rangle_K \\ \langle \Phi_i, F_x\Phi_i \rangle_K & \langle H\Phi_i, F_x\Phi_i \rangle_K & \langle F_x\Phi_i, F_x\Phi_i \rangle_K \end{bmatrix} \quad (13)$$

## 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]:

$$u(r, \theta) = \sum_{j=0}^{\infty} r^{\lambda_j} (k_1^j u_I^j + k_2^j u_{II}^j + k_3^j u_{III}^j) \quad \lambda_j = \frac{1}{2} + j \quad (14)$$

$$\left\{ \begin{array}{l} \underline{u}_I^0 = \frac{1}{2\mu} \sqrt{\frac{r}{2\pi}} (\cos \frac{\theta}{2} (\kappa - \cos \theta) \underline{e}^1 + \sin \frac{\theta}{2} (\kappa - \cos \theta) \underline{e}^2) \\ \underline{u}_{II}^0 = \frac{1}{2\mu} \sqrt{\frac{r}{2\pi}} (\sin \frac{\theta}{2} (2 + \kappa + \cos \theta) \underline{e}^1 + \cos \frac{\theta}{2} (2 - \kappa - \cos \theta) \underline{e}^2) \\ \underline{u}_{III}^0 = \frac{2}{\mu} \sqrt{\frac{r}{2\pi}} \sin \frac{\theta}{2} \underline{e}^3 \end{array} \right. \quad (15)$$

where $\kappa = 3 - 4\nu$ (under the plane strain hypothesis) is Kosolov's constant and $\nu$ is Poisson's ratio.

The first eigenvalue $\lambda_0 = 0.5$ represents the less regular term of this expansion series. The associated functions ($\propto \sqrt{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 $\lambda_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 crack-tip, 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:
If the degrees of freedom for each node are adjacent (ranked contiguously), the shape of $P_C$ is:

$$P_C = \begin{bmatrix} P_C^1 & 0 & 0 & \cdots & \cdots \\ 0 & P_C^2 & 0 & \cdots & 0 \\ 0 & 0 & P_C^3 & \cdots & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ \cdots & 0 & \cdots & \cdots & P_C^i \\ \cdots & \cdots & \cdots & \cdots & \ddots \end{bmatrix} \quad (10)$$

The blocks of $P_C$ express as follows:

$$P_C^i = \begin{cases} S_i^{-1} & \text{if the node is enriched} \\ I_d & \text{if not} \end{cases} \quad (11)$$

where the inverse terms are computed from Cholesky's factorization of local blocks of K:

$$K_{loc}^i = S_i^T S_i \quad (12)$$

$$K_{loc}^i = \begin{bmatrix} \langle \Phi_i, \Phi_i \rangle_K & \langle \Phi_i, H\Phi_i \rangle_K & \langle \Phi_i, F_x\Phi_i \rangle_K \\ \langle \Phi_i, H\Phi_i \rangle_K & \langle H\Phi_i, H\Phi_i \rangle_K & \langle H\Phi_i, F_x\Phi_i \rangle_K \\ \langle \Phi_i, F_x\Phi_i \rangle_K & \langle H\Phi_i, F_x\Phi_i \rangle_K & \langle F_x\Phi_i, F_x\Phi_i \rangle_K \end{bmatrix} \quad (13)$$

## 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]:

$$u(r, \theta) = \sum_{j=0}^{\infty} r^{\lambda_j} (k_1^j u_I^j + k_2^j u_{II}^j + k_3^j u_{III}^j) \quad \lambda_j = \frac{1}{2} + j \quad (14)$$

$$\left\{ \begin{array}{l} \underline{u}_I^0 = \frac{1}{2\mu} \sqrt{\frac{r}{2\pi}} (\cos \frac{\theta}{2} (\kappa - \cos \theta) \underline{e}^1 + \sin \frac{\theta}{2} (\kappa - \cos \theta) \underline{e}^2) \\ \underline{u}_{II}^0 = \frac{1}{2\mu} \sqrt{\frac{r}{2\pi}} (\sin \frac{\theta}{2} (2 + \kappa + \cos \theta) \underline{e}^1 + \cos \frac{\theta}{2} (2 - \kappa - \cos \theta) \underline{e}^2) \\ \underline{u}_{III}^0 = \frac{2}{\mu} \sqrt{\frac{r}{2\pi}} \sin \frac{\theta}{2} \underline{e}^3 \end{array} \right. \quad (15)$$

where $\kappa = 3 - 4\nu$ (under the plane strain hypothesis) is Kosolov's constant and $\nu$ is Poisson's ratio.

The first eigenvalue $\lambda_0 = 0.5$ represents the less regular term of this expansion series. The associated functions ($\propto \sqrt{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 $\lambda_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 crack-tip, 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:
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Plane crack local basis at the crack-tip. Definition of cylindrical coordinates at the crack-tip for plane crack</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 5** Sets of enriched nodes around the crack. If the crack ends in the middle of the domain, additional enrichment functions are needed to model the singularity at the vicinity of the crack-tip. **a** illustration of the sets of nodes for a linear mesh (for which $l_{CT} = l_{CT,1}$), **b** illustration of the sets of nodes for a quadratic mesh (for which $l_{CT} = l_{CT,2}$) and **c** illustration of the sets $l_{CT,1}$ and $\delta_{l_{CT,1}}$ for a quadratic mesh</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Plane crack local basis at the crack-tip. Definition of cylindrical coordinates at the crack-tip for plane crack</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5** Sets of enriched nodes around the crack. If the crack ends in the middle of the domain, additional enrichment functions are needed to model the singularity at the vicinity of the crack-tip. **a** illustration of the sets of nodes for a linear mesh (for which $l_{CT} = l_{CT,1}$), **b** illustration of the sets of nodes for a quadratic mesh (for which $l_{CT} = l_{CT,2}$) and **c** illustration of the sets $l_{CT,1}$ and $l_{CT,2}$ for a quadratic mesh</figcaption></figure>
* 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
* 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,

$$V^h = \left\{ u^h = \sum_{i \in I} a_i \Phi_i + \sum_{j \in I_H} b_j H \Phi_j \right\} \qquad (16)$$

where $\Phi_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 $\Gamma$. We recall that $H$ is defined as follows:

$$\forall x \in \Omega, x \notin \Gamma, \quad H(x) = \begin{cases} -1, & \text{if } \ln s(x) < 0 \\ +1, & \text{if } \ln s(x) > 0 \end{cases} \qquad (17)$$

Let's consider the basic case where the whole domain is divided into two parts as in Fig. 3:

* $\Omega_{(H=+1)} = \Omega_+$

* $\Omega_{(H=-1)} = \Omega_-$.

For which,

* $u^h|_{\Omega_+}$ represents the approximation of the solution on the dimain $\Omega_+$ extended by continuity to the domain $\Omega_-$ over the interface $\Gamma$,

* $u^h|_{\Omega_-}$ represents the approximation of the solution on the dimain $\Omega_-$ extended by continuity to the domain $\Omega_+$ over the interface $\Gamma$.

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:

$$\forall x \in \Gamma, \quad \Delta u^h(x) = \left\{ u^h|_{\Omega_+} - u^h|_{\Omega_-} \right\} (x) \qquad (18)$$

In order to simplify the notations in Table 1, we define for a given node $j \in I_H$, $\Phi_{j,1} = \chi_{\Omega_+} \Phi_j$ and $\Phi_{j,2} = \chi_{\Omega_-} \Phi_j$, where $\chi_A$ denotes the characteristic function of the set A. Using these notations, we have:

$$H \Phi_j - H(x_j) \Phi_j = \begin{cases} +2\chi_{\Omega_+} \Phi_j = +2\Phi_{j,1}, & \text{if } j \in \Omega_- \\ -2\chi_{\Omega_-} \Phi_j = -2\Phi_{j,2}, & \text{if } j \in \Omega_+ \end{cases}$$
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,

$$V^h = \left\{ u^h = \sum_{i \in I} a_i \Phi_i + \sum_{j \in I_H} b_j H \Phi_j \right\} \quad (16)$$

where $\Phi_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 $\Gamma$. We recall that $H$ is defined as follows:

$$\forall x \in \Omega, x \notin \Gamma, \quad H(x) = \begin{cases} -1, & \text{if } \ln s(x) < 0 \\ +1, & \text{if } \ln s(x) > 0 \end{cases} \quad (17)$$

Let's consider the basic case where the whole domain is divided into two parts as in Fig. 3:

* $\Omega_{(H=+1)} = \Omega_+$

* $\Omega_{(H=-1)} = \Omega_-$.

For which,

* $u^h|_{\Omega_+}$ represents the approximation of the solution on the dimain $\Omega_+$ extended by continuity to the domain $\Omega_-$ over the interface $\Gamma$,

* $u^h|_{\Omega_-}$ represents the approximation of the solution on the dimain $\Omega_-$ extended by continuity to the domain $\Omega_+$ over the interface $\Gamma$.

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:

$$\forall x \in \Gamma, \quad \Delta u^h(x) = \left\{ u^h|_{\Omega_+} - u^h|_{\Omega_-} \right\} (x) \quad (18)$$

In order to simplify the notations in Table 1, we define for a given node $j \in I_H$, $\Phi_{j,1} = \chi_{\Omega_+} \Phi_j$ and $\Phi_{j,2} = \chi_{\Omega_-} \Phi_j$, where $\chi_A$ denotes the characteristic function of the set A. Using these notations, we have:

$$H \Phi_j - H(x_j) \Phi_j = \begin{cases} +2\chi_{\Omega_+} \Phi_j = +2\Phi_{j,1}, & \text{if } j \in \Omega_- \\ -2\chi_{\Omega_-} \Phi_j = -2\Phi_{j,2}, & \text{if } j \in \Omega_+ \end{cases}$$
<table><caption>Table 1 Strong discontinuity enrichments comparison</caption><thead><tr><th></th><th>X-FEM [7]</th><th>Domain enrich. [33]</th><th>Shifted enrich. [26]</th></tr></thead><tbody><tr><td>Linear shape functions on enriched nodes</td><td></td><td></td><td></td></tr><tr><td>Displacement approximation</td><td>u<sup>h</sup>|<sub>Ω+</sub> = &sum;<sub>j&isin;H</sub> a<sub>j</sub>&Phi;<sub>j</sub> + &sum;<sub>j&isin;H</sub> b<sub>j</sub>&Phi;<sub>j</sub><br>u<sup>h</sup>|<sub>Ω-</sub> = &sum;<sub>j&isin;H</sub> a<sub>j</sub>&Phi;<sub>j</sub> - &sum;<sub>j&isin;H</sub> b<sub>j</sub>&Phi;<sub>j</sub></td><td>u<sup>h</sup>|<sub>Ω+</sub> = &sum;<sub>i&isin;H</sub> a<sub>i</sub>&Phi;<sub>i</sub> + &sum;<sub>j&isin;H</sub> a<sub>j</sub>&Phi;<sub>j</sub><br>u<sup>h</sup>|<sub>Ω-</sub> = &sum;<sub>i&isin;H</sub> a<sub>i</sub>&Phi;<sub>i</sub> + &sum;<sub>j&isin;H</sub> a<sub>j</sub>&Phi;<sub>j</sub></td><td>u<sup>h</sup>|<sub>Ω+</sub> = &sum;<sub>i&isin;H</sub> c<sub>i</sub>&Phi;<sub>i</sub> - &sum;<sub>j&isin;H</sub> 2d<sub>j</sub>&Phi;<sub>j</sub><br>u<sup>h</sup>|<sub>Ω-</sub> = &sum;<sub>i&isin;H</sub> c<sub>i</sub>&Phi;<sub>i</sub> + &sum;<sub>j&isin;H</sub> 2d<sub>j</sub>&Phi;<sub>j</sub></td></tr><tr><td>Jump approximation</td><td>&Delta;u<sup>h</sup> = &sum;<sub>j&isin;H</sub> 2b<sub>j</sub>&Phi;<sub>j</sub></td><td>&Delta;u<sup>h</sup> = &sum;<sub>j&isin;H</sub> a<sub>j</sub>&Phi;<sub>j</sub><sub>1</sub> - a<sub>j</sub>&Phi;<sub>j</sub><sub>2</sub>&Phi;<sub>j</sub><sub>2</sub></td><td>&Delta;u<sup>h</sup> = - &sum;<sub>j&isin;H</sub> 2d<sub>j</sub>&Phi;<sub>j</sub><sub>2</sub> - &sum;<sub>j&isin;H</sub> 2d<sub>j</sub>&Phi;<sub>j</sub><sub>1</sub></td></tr></tbody></table>

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
<table><caption>Table 1 Strong discontinuity enrichments comparison</caption><thead><tr><th></th><th>X-FEM [7]</th><th>Domain enrich. [33]</th><th>Shifted enrich. [26]</th></tr></thead><tbody><tr><td>Linear shape functions on enriched nodes</td><td></td><td></td><td></td></tr><tr><td>Displacement approximation</td><td>u<sup>h</sup>|<sub>Ω+</sub> = &sum;<sub>j&isin;I<sub>H</sub></sub> a<sub>j</sub>&Phi;<sub>j</sub> + &sum;<sub>j&isin;I<sub>H</sub></sub> b<sub>j</sub>&Phi;<sub>j</sub><br>u<sup>h</sup>|<sub>Ω-</sub> = &sum;<sub>j&isin;I<sub>H</sub></sub> a<sub>j</sub>&Phi;<sub>j</sub> - &sum;<sub>j&isin;I<sub>H</sub></sub> b<sub>j</sub>&Phi;<sub>j</sub></td><td>u<sup>h</sup>|<sub>Ω+</sub> = &sum;<sub>i&isin;J<sub>H</sub></sub> a<sub>i</sub>&Phi;<sub>i</sub> + &sum;<sub>j&isin;I<sub>H</sub></sub> a<sub>j</sub>&Phi;<sub>j</sub><br>u<sup>h</sup>|<sub>Ω-</sub> = &sum;<sub>i&isin;J<sub>H</sub></sub> a<sub>i</sub>&Phi;<sub>i</sub> + &sum;<sub>j&isin;I<sub>H</sub></sub> a<sub>j</sub>&Phi;<sub>j</sub></td><td>u<sup>h</sup>|<sub>Ω+</sub> = &sum;<sub>i&isin;I<sub>H</sub></sub> c<sub>i</sub>&Phi;<sub>i</sub> - &sum;<sub>j&isin;I<sub>H</sub>&cap;Ω<sub>-</sub></sub> 2d<sub>j</sub>&Phi;<sub>j</sub><br>u<sup>h</sup>|<sub>Ω-</sub> = &sum;<sub>i&isin;I<sub>H</sub></sub> c<sub>i</sub>&Phi;<sub>i</sub> + &sum;<sub>j&isin;I<sub>H</sub>&cap;Ω<sub>+</sub></sub> 2d<sub>j</sub>&Phi;<sub>j</sub></td></tr><tr><td>Jump approximation</td><td>&Delta;u<sup>h</sup> = &sum;<sub>j&isin;I<sub>H</sub></sub> 2b<sub>j</sub>&Phi;<sub>j</sub></td><td>&Delta;u<sup>h</sup> = &sum;<sub>j&isin;I<sub>H</sub></sub> a<sub>j</sub>&Phi;<sub>j</sub><sub>1</sub> - a<sub>j</sub>&Phi;<sub>j</sub><sub>2</sub>&Phi;<sub>j</sub><sub>2</sub></td><td>&Delta;u<sup>h</sup> = - &sum;<sub>j&isin;I<sub>H</sub>&cap;Ω<sub>-</sub></sub> 2d<sub>j</sub>&Phi;<sub>j</sub><sub>2</sub> - &sum;<sub>j&isin;I<sub>H</sub>&cap;Ω<sub>+</sub></sub> 2d<sub>j</sub>&Phi;<sub>j</sub><sub>1</sub></td></tr></tbody></table>

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 $\Omega$ (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 $\Gamma_D$ and Neumann boundary conditions are applied on $\Gamma_N$.

The space of admissible displacements is:

$$V = \{v \in H^1(\Omega); v = 0 \text{ on } \Gamma_D\} \qquad (19)$$

<table><caption>**Table 2** Change of variables when switching enrichment strategies. $V_j \in I_H$, the following changes of variables hold</caption><thead><tr><th></th><th>X-FEM</th><th>Domain enrich.</th><th>Shifted enrich.</th></tr></thead><tbody><tr><td>X-FEM</td><td>$\begin{cases} a_j = a_j \\ b_j = b_j \end{cases}$</td><td>$\begin{cases} \alpha_{j,1} = a_j + b_j \\ \alpha_{j,2} = a_j - b_j \end{cases}$</td><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} c_j = a_j + b_j \\ d_j = b_j \end{cases} \end{cases}$<br>$\begin{cases} j \in \Omega_- \\ \begin{cases} c_j = a_j - b_j \\ d_j = b_j \end{cases} \end{cases}$</td></tr><tr><td>Domain enrich.</td><td>$\begin{cases} a_j = \frac{\alpha_{j,1} + \alpha_{j,2}}{2} \\ b_j = \frac{\alpha_{j,1} - \alpha_{j,2}}{2} \end{cases}$</td><td>$\begin{cases} \alpha_{j,1} = \alpha_{j,1} \\ \alpha_{j,2} = \alpha_{j,1} \end{cases}$</td><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} c_j = \alpha_{j,1} \\ d_j = \frac{\alpha_{j,1} - \alpha_{j,2}}{2} \end{cases} \end{cases}$<br>$\begin{cases} j \in \Omega_- \\ \begin{cases} c_j = \alpha_{j,2} \\ d_j = \frac{\alpha_{j,1} - \alpha_{j,2}}{2} \end{cases} \end{cases}$</td></tr><tr><td>Shifted enrich.</td><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} a_j = c_j - d_j \\ b_j = d_j \end{cases} \end{cases}$</td><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} \alpha_{j,1} = c_j \\ \alpha_{j,2} = c_j - 2d_j \end{cases} \end{cases}$</td><td>$\begin{cases} j \in \Omega_- \\ \begin{cases} a_j = c_j + 2d_j \\ \alpha_{j,2} = c_j \end{cases} \end{cases}$</td></tr></tbody></table>
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 $\Omega$ (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 $\Gamma_D$ and Neumann boundary conditions are applied on $\Gamma_N$.

The space of admissible displacements is:

$$V = \{v \in H^1(\Omega); v = 0 \text{ on } \Gamma_D\} \qquad (19)$$

<table><caption>**Table 2** Change of variables when switching enrichment strategies. $V_j \in I_H$, the following changes of variables hold</caption><thead><tr><th></th><th>X-FEM</th><th>Domain enrich.</th><th>Shifted enrich.</th></tr></thead><tbody><tr><th>X-FEM</th><td>$\begin{cases} a_j = a_j \\ b_j = b_j \end{cases}$</td><td>$\begin{cases} \alpha_{j,1} = a_j + b_j \\ \alpha_{j,2} = a_j - b_j \end{cases}$</td><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} c_j = a_j + b_j \\ d_j = b_j \end{cases} \end{cases}$<br>..................<br>$\begin{cases} j \in \Omega_- \\ \begin{cases} c_j = a_j - b_j \\ d_j = b_j \end{cases} \end{cases}$</td></tr><tr><th>Domain enrich.</th><td>$\begin{cases} a_j = \frac{\alpha_{j,1} + \alpha_{j,2}}{2} \\ b_j = \frac{\alpha_{j,1} - \alpha_{j,2}}{2} \end{cases}$</td><td>$\begin{cases} \alpha_{j,1} = \alpha_{j,1} \\ \alpha_{j,2} = \alpha_{j,1} \end{cases}$</td><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} c_j = \alpha_{j,1} \\ d_j = \frac{\alpha_{j,1} - \alpha_{j,2}}{2} \end{cases} \end{cases}$<br>..................<br>$\begin{cases} j \in \Omega_- \\ \begin{cases} c_j = \alpha_{j,2} \\ d_j = \frac{\alpha_{j,1} - \alpha_{j,2}}{2} \end{cases} \end{cases}$</td></tr><tr><th>Shifted enrich.</th><td>$\begin{cases} j \in \Omega_+ \\ \begin{cases} a_j = c_j - d_j \\ b_j = d_j \end{cases} \end{cases}$</td><td>$\begin{cases} j \in \Omega_- \\ \begin{cases} a_j = c_j + d_j \\ b_j = d_j \end{cases} \end{cases}$</td><td>$\begin{cases} j \in \Omega_- \\ \begin{cases} a_j = c_j + 2d_j \\ \alpha_{j,2} = c_j \end{cases} \end{cases}$</td></tr></tbody></table>
<figure><img src="image_9.png" /></figure>

The weak form of the equilibrium problem is:

$$P: \left\{ \begin{array}{l} \text{Find } u \in V \text{ such that } a(u, v) = l(v) \quad \forall v \in V, \\ a(u, v) = \int_{\Omega} \sigma(u) : \varepsilon(v) d\Omega, \\ l(v) = \int_{\Omega} f \cdot v \, d\Omega + \int_{\Gamma_N} g \cdot v \, d\Gamma, \\ \sigma(u) = \lambda \operatorname{tr}(\varepsilon(u)) + 2\mu\varepsilon(u), \end{array} \right. \qquad (20)$$

where $u$ is the displacement, $\sigma$ is the Cauchy stress, $\varepsilon$ is the strain, $f$ is the body force applied on $\omega$, $g$ is the traction applied on $\Gamma_N$ and $\lambda$ and $\mu$ 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.
<figure><img src="image_10.png" /><figcaption>Fig. 6 Example of a cracked domain. Definition of notations used to label the domain and its boundaries</figcaption></figure>

The weak form of the equilibrium problem is:

$$ P: \left\{ \begin{array}{l} \text{Find } u \in V \text{ such that } a(u, v) = l(v) \quad \forall v \in V, \\ a(u, v) = \int_{\Omega} \sigma(u) : \varepsilon(v) d\Omega, \\ l(v) = \int_{\Omega} f \cdot v \, d\Omega + \int_{\Gamma_N} g \cdot v \, d\Gamma, \\ \sigma(u) = \lambda \operatorname{tr}(\varepsilon(u)) + 2\mu\varepsilon(u), \end{array} \right. \qquad (20) $$

where $u$ is the displacement, $\sigma$ is the Cauchy stress, $\varepsilon$ is the strain, $f$ is the body force applied on $\omega$, $g$ is the traction applied on $\Gamma_N$ and $\lambda$ and $\mu$ 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.
<figure><img src="image_10.png" /></figure>

### “Straightforward” singular enrichment class

Let $V^h$ be the approximation space using vector asymptotic functions to describe the crack kinematics:

$$V^h = \left\{ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j \right\} + \sum_{i \in I_{CT}} \sum_{\alpha=1}^{d} c_{i,\alpha} \Phi_i K_{\alpha} \quad (21)$$

where $d$ is equal to two in the bidimensional case and three in the tridimensional case, ($\underline{E}_1, \underline{E}_2, \underline{E}_3$) is the usual Cartesian basis, $K_\alpha$ 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_\alpha$ are defined as follows:

$$\left\{ \begin{array}{l} \underline{K}_1 = \underline{u}_I^0 \\ \underline{K}_2 = \underline{u}_{II}^0 \\ \underline{K}_3 = \underline{u}_{III}^0 \end{array} \right. \qquad (22)$$

with $\underline{u}_I^0, \underline{u}_{II}^0, \underline{u}_{III}^0$ given by (15).

Let $V_{CUTOFF}^h(W^{2,\infty})$ be the approximation space defined with a cutoff function $\chi \in W^{2,\infty}$ [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:
<figure><img src="image_11.png" /></figure>

### “Straightforward” singular enrichment class

Let $V^h$ be the approximation space using vector asymptotic functions to describe the crack kinematics:

$$V^h = \left\{ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j \right\} + \sum_{i \in I_{CT}} \sum_{\alpha=1}^{d} c_{i,\alpha} \Phi_i K_{\alpha} \quad (21)$$

where $d$ is equal to two in the bidimensional case and three in the tridimensional case, $(\underline{E}_1, \underline{E}_2, \underline{E}_3)$ is the usual Cartesian basis, $K_\alpha$ 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_\alpha$ are defined as follows:

$$\left\{ \begin{array}{l} \underline{K}_1 = \underline{u}_I^0 \\ \underline{K}_2 = \underline{u}_{II}^0 \\ \underline{K}_3 = \underline{u}_{III}^0 \end{array} \right. \qquad (22)$$

with $\underline{u}_I^0, \underline{u}_{II}^0, \underline{u}_{III}^0$ given by (15).

Let $V_{CUTOFF}^h(W^{2,\infty})$ be the approximation space defined with a cutoff function $\chi \in W^{2,\infty}$ [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:
$$V_{CUTOFF}^{h}(W^{2,\infty}) = \left\{ v^{h} = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_{i} E_{j} + \sum_{i \in I_{H}} \sum_{j=1}^{d} b_{i,j} H \Phi_{i} E_{j} + \sum_{\alpha=1}^{d} c_{\alpha} \chi \underline{K}_{\alpha} \right\} \quad (23)$$

and let $V_{XFEM}^h$ be the standard X-FEM approximation space [7]:

$$V_{XFEM}^{h} = \left\{ v^{h} = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_{i} E_{j} + \sum_{i \in I_{H}} \sum_{j=1}^{d} b_{i,j} H \Phi_{i} E_{j} + \sum_{k \in I_{CT}} \sum_{\alpha=1,4} \sum_{j=1}^{d} c_{k,\alpha}^{j} F^{\alpha} \Phi_{k} E_{j} \right\} \quad (24)$$

where $F^\alpha$ are defined as:

$$\left\{ \begin{aligned} F^1 &= \sqrt{r} \sin \frac{\theta}{2} \\ F^2 &= \sqrt{r} \cos \frac{\theta}{2} \\ F^3 &= \sqrt{r} \sin \frac{\theta}{2} \sin \theta \\ F^4 &= \sqrt{r} \cos \frac{\theta}{2} \sin \theta \end{aligned} \right. \quad (25)$$

We remark that:

$$\left\{ \begin{aligned} \sqrt{r} \cos \frac{\theta}{2} \cos \theta &= F^2 - F^3 \\ \sqrt{r} \sin \frac{\theta}{2} \cos \theta &= F^4 - F^1 \end{aligned} \right. \quad (26)$$

Consequently, the vector functions $\underline{K}_\alpha$ and the scalar functions $F^\alpha$ are related by the following linear relations:

$$\left\{ \begin{aligned} \underline{K}_1 &= \frac{1}{2\mu\sqrt{2\pi}} [(\kappa - 1)F^2 + F^3] \underline{e}^1 + \frac{1}{2\mu\sqrt{2\pi}} [(\kappa + 1)F^1 - F^4] \underline{e}^2 \\ \underline{K}_2 &= \frac{1}{2\mu\sqrt{2\pi}} [(\kappa + 1)F^1 + F^4] \underline{e}^1 + \frac{1}{2\mu\sqrt{2\pi}} [(1 - \kappa)F^2 + F^3] \underline{e}^2 \\ \underline{K}_3 &= \frac{2}{\mu\sqrt{2\pi}} F^1 \underline{e}^3 \end{aligned} \right. \quad (27)$$

**Lemma 1** $V_{CUTOFF}^h(W_h^{2,\infty})$ is a subspace of $V^h$ and $V^h$ is a subspace of $V_{XFEM}^h$, i.e. the three approximation spaces are related as follows,

$$V_{CUTOFF}^{h}(W_{h}^{2,\infty}) \subset V^{h} \subset V_{XFEM}^{h}$$

where $W_h^{2,\infty}$ is the subspace of interpolated cutoff functions.

*Proof* The proof relies on the work of [29].

$$\bullet \ V_{CUTOFF}^{h}(W_{h}^{2,\infty}) \subset V^{h}:$$

For $v^h \in V_{CUTOFF}^h(W_h^{2,\infty})$,

$$v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{\alpha=1}^{d} c_{\alpha} \chi_h \underline{K}_{\alpha}$$

$\chi_h$ is the interpolation of a cutoff function $\chi$ which vanishes outside $I_{CT}$,

$$\chi_h = \sum_{i \in I_{CT}} \chi(x_i) \Phi_i$$

Then,

$$v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{i \in I_{CT}} \sum_{\alpha=1}^{d} c_{\alpha} \chi(x_i) \Phi_i \underline{K}_{\alpha}$$

which leads to,

$$v^h \in V^h$$
$$V_{\text{CUTOFF}}^{h}(W^{2,\infty}) = \left\{ v^{h} = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_{i} E_{j} + \sum_{i \in I_{H}} \sum_{j=1}^{d} b_{i,j} H \Phi_{i} E_{j} + \sum_{\alpha=1}^{d} c_{\alpha} \chi \underline{K}_{\alpha} \right\} \quad (23)$$

and let $V_{\text{XFEM}}^h$ be the standard X-FEM approximation space [7]:

$$V_{\text{XFEM}}^{h} = \left\{ \begin{array}{l} v^{h} = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_{i} E_{j} + \sum_{i \in I_{H}} \sum_{j=1}^{d} b_{i,j} H \Phi_{i} E_{j} + \\ \qquad \sum_{k \in I_{CT}} \sum_{\alpha=1,4} \sum_{j=1}^{d} c_{k,\alpha}^{j} F^{\alpha} \Phi_{k} E_{j} \end{array} \right\} \quad (24)$$

where $F^\alpha$ are defined as:

$$\left\{ \begin{aligned} F^1 &= \sqrt{r} \sin \frac{\theta}{2} \\ F^2 &= \sqrt{r} \cos \frac{\theta}{2} \\ F^3 &= \sqrt{r} \sin \frac{\theta}{2} \sin \theta \\ F^4 &= \sqrt{r} \cos \frac{\theta}{2} \sin \theta \end{aligned} \right. \quad (25)$$

We remark that:

$$\left\{ \begin{aligned} \sqrt{r} \cos \frac{\theta}{2} \cos \theta &= F^2 - F^3 \\ \sqrt{r} \sin \frac{\theta}{2} \cos \theta &= F^4 - F^1 \end{aligned} \right. \quad (26)$$

Consequently, the vector functions $\underline{K}_\alpha$ and the scalar functions $F^\alpha$ are related by the following linear relations:

$$\left\{ \begin{aligned} \underline{K}_1 &= \frac{1}{2\mu\sqrt{2\pi}} [(\kappa - 1)F^2 + F^3] \underline{e}^1 + \frac{1}{2\mu\sqrt{2\pi}} [(\kappa + 1)F^1 - F^4] \underline{e}^2 \\ \underline{K}_2 &= \frac{1}{2\mu\sqrt{2\pi}} [(\kappa + 1)F^1 + F^4] \underline{e}^1 + \frac{1}{2\mu\sqrt{2\pi}} [(1 - \kappa)F^2 + F^3] \underline{e}^2 \\ \underline{K}_3 &= \frac{2}{\mu\sqrt{2\pi}} F^1 \underline{e}^3 \end{aligned} \right. \quad (27)$$

**Lemma 1** $V_{\text{CUTOFF}}^h(W_h^{2,\infty})$ is a subspace of $V^h$ and $V^h$ is a subspace of $V_{\text{XFEM}}^h$, i.e. the three approximation spaces are related as follows,

$$V_{\text{CUTOFF}}^{h}(W_{h}^{2,\infty}) \subset V^{h} \subset V_{\text{XFEM}}^{h}$$

where $W_h^{2,\infty}$ is the subspace of interpolated cutoff functions.

*Proof* The proof relies on the work of [29].

$$\bullet \ V_{\text{CUTOFF}}^{h}(W_{h}^{2,\infty}) \subset V^{h}:$$

For $v^h \in V_{\text{CUTOFF}}^h(W_h^{2,\infty})$,

$$v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{\alpha=1}^{d} c_{\alpha} \chi_h \underline{K}_{\alpha}$$

$\chi_h$ is the interpolation of a cutoff function $\chi$ which vanishes outside $I_{CT}$,

$$\chi_h = \sum_{i \in I_{CT}} \chi(x_i) \Phi_i$$

Then,

$$v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{i \in I_{CT}} \sum_{\alpha=1}^{d} c_{\alpha} \chi(x_i) \Phi_i \underline{K}_{\alpha}$$

which leads to,

$$v^h \in V^h$$
$$ \bullet \quad V^h \subset V_{XFEM}^h : $$

For $v^h \in V^h$,

$$ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{k \in I_{CT}} \sum_{\alpha=1}^{d} c_{k,\alpha} \Phi_k \underline{K}_{\alpha} $$

* Equation (27) shows that $\underline{K}_{\alpha}$ can be expressed as $\underline{K}_{\alpha} = \sum_{l=1}^{4} \sum_{m=1}^{d} k_{l,m,\alpha} F^l \underline{e}^m$, where $k_{l,m,\alpha}$ is a constant tensor.

In case of a straight crack, the local basis is constant.

Thus, the local basis expresses in cartesian coordinates as, $\underline{e}^m = \sum_{j=1}^{d} \mu_{m,j} \underline{E}_j$, where $\mu_{m,j}$ is a constant tensor.

Then,

$$ \underline{K}_{\alpha} = \sum_{l=1}^{4} \sum_{m=1}^{d} \sum_{j=1}^{d} k_{l,m,\alpha} \mu_{m,j} F^l \underline{E}_j $$

$$ \sum_{\alpha=1}^{d} c_{k,\alpha} \Phi_k \underline{K}_{\alpha} = \sum_{\alpha=1}^{d} \sum_{l=1}^{4} \sum_{m=1}^{d} \sum_{j=1}^{d} c_{k,\alpha} k_{l,m,\alpha} \mu_{m,j} F^l \Phi_k \underline{E}_j $$

$$ = \sum_{l=1}^{4} \sum_{j=1}^{d} \left( \sum_{\alpha=1}^{d} \sum_{m=1}^{d} c_{k,\alpha} k_{l,m,\alpha} \mu_{m,j} \right) F^l \Phi_k \underline{E}_j $$

With a suitable change of variable,

$$ c_{k,l}^j = \sum_{\alpha=1}^{d} \sum_{m=1}^{d} c_{k,\alpha} k_{l,m,\alpha} \mu_{m,j} $$

we have finally,

$$ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i \underline{E}_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i \underline{E}_j + \sum_{k \in I_{CT}} \sum_{l=1}^{4} \sum_{j=1}^{d} c_{k,l}^j F^l \Phi_k \underline{E}_j $$

which means that:

$$ v^h \in V_{XFEM}^h $$

* 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 $(\underline{e}_k^1, \underline{e}_k^2, \underline{e}_k^3)$, such as,

$$ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i \underline{E}_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i \underline{E}_j + \sum_{k \in I_{CT}} \sum_{\alpha=1}^{d} c_{k,\alpha} \Phi_k \underline{K}_{\alpha,k} $$

$$ \underline{K}_{\alpha,k} = \sum_{l=1}^{4} \sum_{m=1}^{d} k_{l,m,\alpha} F^l \underline{e}_k^m $$

$$ \underline{e}_k^m = \sum_{j=1}^{d} \mu_{m,j,k} \underline{E}_j $$

The straight crack proof still holds here, with the additional node index “k”.

□
$$ \bullet \quad V^h \subset V_{XFEM}^h : $$

For $v^h \in V^h$,

$$ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{k \in I_{CT}} \sum_{\alpha=1}^{d} c_{k,\alpha} \Phi_k \underline{K}_{\alpha} $$

* Equation (27) shows that $\underline{K}_{\alpha}$ can be expressed as $\underline{K}_{\alpha} = \sum_{l=1}^{4} \sum_{m=1}^{d} k_{l,m,\alpha} F^l \underline{e}^m$, where $k_{l,m,\alpha}$ is a constant tensor.

In case of a straight crack, the local basis is constant.

Thus, the local basis expresses in cartesian coordinates as, $\underline{e}^m = \sum_{j=1}^{d} \mu_{m,j} E_j$, where $\mu_{m,j}$ is a constant tensor.

Then,

$$ \underline{K}_{\alpha} = \sum_{l=1}^{4} \sum_{m=1}^{d} \sum_{j=1}^{d} k_{l,m,\alpha} \mu_{m,j} F^l \underline{E}_j $$

$$ \sum_{\alpha=1}^{d} c_{k,\alpha} \Phi_k \underline{K}_{\alpha} = \sum_{\alpha=1}^{d} \sum_{l=1}^{4} \sum_{m=1}^{d} \sum_{j=1}^{d} c_{k,\alpha} k_{l,m,\alpha} \mu_{m,j} F^l \Phi_k E_j $$

$$ = \sum_{l=1}^{4} \sum_{j=1}^{d} \left( \sum_{\alpha=1}^{d} \sum_{m=1}^{d} c_{k,\alpha} k_{l,m,\alpha} \mu_{m,j} \right) F^l \Phi_k E_j $$

With a suitable change of variable,

$$ c_{k,l}^j = \sum_{\alpha=1}^{d} \sum_{m=1}^{d} c_{k,\alpha} k_{l,m,\alpha} \mu_{m,j} $$

we have finally,

$$ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{k \in I_{CT}} \sum_{l=1}^{4} \sum_{j=1}^{d} c_{k,l}^j F^l \Phi_k E_j $$

which means that:

$$ v^h \in V_{XFEM}^h $$

* 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 $(\underline{e}_k^1, \underline{e}_k^2, \underline{e}_k^3)$, such as,

$$ v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} H \Phi_i E_j + \sum_{k \in I_{CT}} \sum_{\alpha=1}^{d} c_{k,\alpha} \Phi_k \underline{K}_{\alpha,k} $$

$$ \underline{K}_{\alpha,k} = \sum_{l=1}^{4} \sum_{m=1}^{d} k_{l,m,\alpha} F^l \underline{e}_k^m $$

$$ \underline{e}_k^m = \sum_{j=1}^{d} \mu_{m,j,k} E_j $$

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]:

$$V_{SGFEM}^h = \left\{ v^h = \sum_{i \in I} \sum_{j=1}^d a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^d \sum_{k=1}^4 b_{i,j,k} H \psi_{i,k} \Phi_i E_j \right. \\ \left. + \sum_{i \in I_{CT}} \sum_{\alpha=1}^d c_{i,\alpha} \Phi_i (\underline{K}_\alpha - \Pi \underline{K}_\alpha) \right\} \quad (28)$$

where $\{\psi_{i,k}\}_{k=1..4} = \{1, x - x_i, y - y_i, z - z_i\}$, and $\Pi \underline{K}_\alpha = \sum_{k \in \{I_{CT} \cup \delta I_{CT}\}} \underline{K}_\alpha (r_k, \theta_k) \Phi_k$ is the usual interpolation.

$\delta I_{CT}$ is defined here as the transition layer of nodes between crack-tip elements and other elements of the mesh for which $\underline{K}_\alpha$ unknowns are set to zero. This allows to connect the enriched layer with the remaining of the mesh and to avoid blending issues [26].

**Lemma 2** $V_{XFEM}^h$ is not a subspace of $V_{SGFEM}^h$ and $V_{SGFEM}^h$ is not a subspace of $V_{XFEM}^h$, so that:

$$V_{XFEM}^h \not\subset V_{SGFEM}^h \text{ and } V_{SGFEM}^h \not\subset V_{XFEM}^h \quad (29)$$

*Proof* — $V_{SGFEM}^h \not\subset V_{XFEM}^h$

It is obvious that $V_{SGFEM}^h \not\subset V_{XFEM}^h$ because $V_{SGFEM}^h$ involves higher order discontinuous polynomials $H \psi_{i,k} \Phi_i$, which are out of $V_{XFEM}^h$.

$$V_{XFEM}^h \not\subset V_{SGFEM}^h$$

* $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k$ belongs to the test function space $V_{XFEM}^h$ with $a_{i,j} = 0, b_{i,j} = 0, c_{k,1}^1 = 1$ and $c_{k,\alpha}^j = 0$ elsewhere (for $\alpha \neq 1$).

* However, $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k$ does not belong to $V_{SGFEM}^h$, which can be shown by contradiction.

* Let us assume $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k$ belongs to $V_{SGFEM}^h$,

$$F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k = \sum_{i \in I} \sum_{j=1}^d a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^d \sum_{k=1}^4 b_{i,j,k} H \psi_{i,k} \Phi_i E_j \\ + \sum_{i \in I_{CT}} \sum_{\alpha=1}^d c_{i,\alpha} \Phi_i (\underline{K}_\alpha - \Pi \underline{K}_\alpha)$$

* 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, so, $\sum_{i \in I} \sum_{j=1}^d a_{i,j} \Phi_i E_j - \sum_{i \in I_{CT}} \sum_{\alpha=1}^d c_{i,\alpha} \Phi_i \Pi \underline{K}_\alpha = 0$ as $\Phi_i \Pi \underline{K}_\alpha$ is also a polynomial.

* $\Phi_i \Pi \underline{K}_\alpha$ introduces at least three higher order polynomials per node,
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]:

$$V_{SGFEM}^h = \left\{ v^h = \sum_{i \in I} \sum_{j=1}^d a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^d \sum_{k=1}^4 b_{i,j,k} H \psi_{i,k} \Phi_i E_j + \sum_{i \in I_{CT}} \sum_{\alpha=1}^d c_{i,\alpha} \Phi_i (\underline{K}_\alpha - \Pi \underline{K}_\alpha) \right\} \quad (28)$$

where $\{\psi_{i,k}\}_{k=1..4} = \{1, x - x_i, y - y_i, z - z_i\}$, and $\Pi \underline{K}_\alpha = \sum_{k \in \{I_{CT} \cup \delta I_{CT}\}} \underline{K}_\alpha(r_k, \theta_k) \Phi_k$ is the usual interpolation.

$\delta I_{CT}$ is defined here as the transition layer of nodes between crack-tip elements and other elements of the mesh for which $\underline{K}_\alpha$ unknowns are set to zero. This allows to connect the enriched layer with the remaining of the mesh and to avoid blending issues [26].

**Lemma 2** $V_{XFEM}^h$ is not a subspace of $V_{SGFEM}^h$ and $V_{SGFEM}^h$ is not a subspace of $V_{XFEM}^h$, so that:

$$V_{XFEM}^h \not\subset V_{SGFEM}^h \text{ and } V_{SGFEM}^h \not\subset V_{XFEM}^h \quad (29)$$

*Proof* — $V_{SGFEM}^h \not\subset V_{XFEM}^h$

It is obvious that $V_{SGFEM}^h \not\subset V_{XFEM}^h$ because $V_{SGFEM}^h$ involves higher order discontinuous polynomials $H \psi_{i,k} \Phi_i$, which are out of $V_{XFEM}^h$.

$$V_{XFEM}^h \not\subset V_{SGFEM}^h$$

* $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k$ belongs to the test function space $V_{XFEM}^h$ with $a_{i,j} = 0, b_{i,j} = 0, c_{k,1}^1 = 1$ and $c_{k,\alpha}^j = 0$ elsewhere (for $\alpha \neq 1$).

* However, $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k$ does not belong to $V_{SGFEM}^h$, which can be shown by contradiction.

* Let us assume $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k$ belongs to $V_{SGFEM}^h$,

$$F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k = \sum_{i \in I} \sum_{j=1}^d a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^d \sum_{k=1}^4 b_{i,j,k} H \psi_{i,k} \Phi_i E_j + \sum_{i \in I_{CT}} \sum_{\alpha=1}^d c_{i,\alpha} \Phi_i (\underline{K}_\alpha - \Pi \underline{K}_\alpha)$$

* 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, so, $\sum_{i \in I} \sum_{j=1}^d a_{i,j} \Phi_i E_j - \sum_{i \in I_{CT}} \sum_{\alpha=1}^d c_{i,\alpha} \Phi_i \Pi \underline{K}_\alpha = 0$ as $\Phi_i \Pi \underline{K}_\alpha$ is also a polynomial.

* $\Phi_i \Pi \underline{K}_\alpha$ introduces at least three higher order polynomials per node,
* Which leads to, $\sum_{\alpha=1}^{d} c_{i,\alpha} \Phi_i \Pi \underline{K}_{\alpha} = 0$,

* Then, $c_{i,\alpha} = 0$,

* And, $a_{i,j} = 0$ because $\sum_{j=1}^{d} a_{i,j} \Phi_i \underline{E}_j = 0$.

* Thus, $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k = 0$ which is contradictory.

<figure><img src="image_11.png" /></figure>

**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 behavior of strongly discontinuous approximations

#### 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 bi-material interface. The arbitrary interface is positioned along the abscissa $x$ ($\varepsilon$) (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.

<figure><img src="image_12.png" /><figcaption>**Fig. 8** 1D case with X-FEM interface and linear elements. Left, the three elements used in the model. Right, the analytical solution expected. The use of a discontinuous function is possible because discontinuous degrees of freedom bear Neumann conditions that enforce the continuity of the solution (as in a contact problem with Lagrange multipliers at the interface)</figcaption></figure>
* Which leads to, $\sum_{\alpha=1}^{d} c_{i,\alpha} \Phi_i \Pi \underline{K}_{\alpha} = 0$,

* Then, $c_{i,\alpha} = 0$,

* And, $a_{i,j} = 0$ because $\sum_{j=1}^{d} a_{i,j} \Phi_i \underline{E}_j = 0$.

* Thus, $F^1 \underline{E}_1 \sum_{k \in I_{CT}} \Phi_k = 0$ which is contradictory. □

**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 behavior of strongly discontinuous approximations

#### 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 bi-material interface. The arbitrary interface is positioned along the abscissa $x$ ($\varepsilon$) (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.

<figure><img src="image_12.png" /><figcaption>**Fig. 8** 1D case with X-FEM interface and linear elements. Left, the three elements used in the model. Right, the analytical solution expected. The use of a discontinuous function is possible because discontinuous degrees of freedom bear Neumann conditions that enforce the continuity of the solution (as in a contact problem with Lagrange multipliers at the interface)</figcaption></figure>
• The problem is the solution of the following differential equation:

$$u''(x) = 0 \quad \forall x \in \Omega = [0, 3] \qquad (30)$$

• With arbitrary Dirichlet boundary conditions:

$$u(0) = \beta \quad u(3) = 3\alpha + \beta \qquad (31)$$

• 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):

$$\left. \frac{du}{dx} \right|_{x(\varepsilon^+)} = \alpha^{-} \quad \left. \frac{du}{dx} \right|_{x(\varepsilon^-)} = \alpha^{+} \quad \alpha^{-} = \alpha^{+} = \alpha \qquad (32)$$

This continuity condition results from the equivalent Lagrangian form of (30) in which a Lagrange multiplier $\lambda$ is used to impose a continuous displacement across the interface. $(u, \lambda)$ is the saddle point of the following functional:

$$L(u, \lambda) = \frac{1}{2} \int_{0}^{3} (u')^{2} dx + \lambda [u(x(\varepsilon^{+})) - u(x(\varepsilon^{-}))] = 0$$

The expected continuous solution satisfying the boundary conditions above, is:

$$u(x) = \alpha x + \beta, \quad \forall x \in [0, 3] \qquad (33)$$

The resulting weak form of the problem is:

$$\left. \begin{array}{l} u \in \prod : \{ w \in H^1/w(0) = \beta, w(3) = 3\alpha + \beta \}, \lambda \in \Re \\ \forall \nu \in \prod_0 : \{ w \in H^1/w(0) = 0, w(3) = 0 \} \\ \forall \mu \in \mathbb{R} \end{array} \right\}$$

$$\mapsto \left\{ \begin{array}{l} \int_0^3 u'v' dx = 0 \\ (\lambda - \alpha^+) \nu(x(\varepsilon^+)) - (\lambda - \alpha^-) \nu(x(\varepsilon^-)) = 0 \\ \mu [u(x(\varepsilon^+)) - u(x(\varepsilon^-))] = 0 \end{array} \right. \qquad (34)$$

The X-FEM discrete linear space is:

$$\prod_h : \left\{ \begin{array}{l} w/w = a_1\Phi_1 + a_2\Phi_2 + b_2H\Phi_2 + a_3\Phi_3 + b_3H\Phi_3 + a_4\Phi_4 \\ a_1 = \beta \quad a_4 = 3\alpha + \beta \end{array} \right\} \qquad (35)$$

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) is:

$$\left[ \begin{array}{cccccc} 1 & -1 & 1 & 0 & 0 & 0 \\ -1 & 2 & -2 + 2\varepsilon & -1 & 1 - 2\varepsilon & 0 \\ 1 & -2 + 2\varepsilon & 2 & 1 - 2\varepsilon & -1 & 0 \\ 0 & -1 & 1 - 2\varepsilon & 2 & 2\varepsilon & -1 \\ 0 & 1 - 2\varepsilon & -1 & 2\varepsilon & 2 & -1 \\ 0 & 0 & 0 & -1 & -1 & 1 \end{array} \right] \left\{ \begin{array}{c} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{array} \right\} \qquad (36)$$

With the domain formulation of the second column of Table 1, the discrete matrix is:

$$\left[ \begin{array}{cccccc} 1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 2 - \varepsilon & 0 & -1 + \varepsilon & 0 & 0 \\ 0 & 0 & \varepsilon & 0 & -\varepsilon & 0 \\ 0 & -1 + \varepsilon & 0 & 1 - \varepsilon & 0 & 0 \\ 0 & 0 & -\varepsilon & 0 & 1 + \varepsilon & -1 \\ 0 & 0 & 0 & 0 & -1 & 1 \end{array} \right] \left\{ \begin{array}{c} \alpha_1 \\ \alpha_{2,1} \\ \alpha_{2,2} \\ \alpha_{3,1} \\ \alpha_{3,2} \\ \alpha_4 \end{array} \right\} \qquad (37)$$
* The problem is the solution of the following differential equation:

$$u''(x) = 0 \quad \forall x \in \Omega = [0, 3] \qquad (30)$$

* With arbitrary Dirichlet boundary conditions:

$$u(0) = \beta \quad u(3) = 3\alpha + \beta \qquad (31)$$

* 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):

$$\left. \frac{du}{dx} \right|_{x(\varepsilon^+)} = \alpha^{-} \quad \left. \frac{du}{dx} \right|_{x(\varepsilon^-)} = \alpha^{+} \quad \alpha^{-} = \alpha^{+} = \alpha \qquad (32)$$

This continuity condition results from the equivalent Lagrangian form of (30) in which a Lagrange multiplier $\lambda$ is used to impose a continuous displacement across the interface. $(u, \lambda)$ is the saddle point of the following functional:

$$L(u, \lambda) = \frac{1}{2} \int_{0}^{3} (u')^{2} dx + \lambda [u(x(\varepsilon^{+})) - u(x(\varepsilon^{-}))] = 0$$

The expected continuous solution satisfying the boundary conditions above, is:

$$u(x) = \alpha x + \beta, \quad \forall x \in [0, 3] \qquad (33)$$

The resulting weak form of the problem is:

$$\left. \begin{array}{l} u \in \prod : \{ w \in H^1/w(0) = \beta, w(3) = 3\alpha + \beta \}, \lambda \in \Re \\ \forall \nu \in \prod_0 : \{ w \in H^1/w(0) = 0, w(3) = 0 \} \\ \forall \mu \in \mathbb{R} \end{array} \right\}$$

$$\mapsto \left\{ \begin{array}{l} \int_0^3 u'v' dx = 0 \\ (\lambda - \alpha^+) \nu(x(\varepsilon^+)) - (\lambda - \alpha^-) \nu(x(\varepsilon^-)) = 0 \\ \mu [u(x(\varepsilon^+)) - u(x(\varepsilon^-))] = 0 \end{array} \right. \qquad (34)$$

The X-FEM discrete linear space is:

$$\prod_h : \left\{ \begin{array}{l} w/w = a_1 \Phi_1 + a_2 \Phi_2 + b_2 H \Phi_2 + a_3 \Phi_3 + b_3 H \Phi_3 + a_4 \Phi_4 \\ a_1 = \beta \quad a_4 = 3\alpha + \beta \end{array} \right\} \qquad (35)$$

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) is:

$$\left[ \begin{array}{cccccc} 1 & -1 & 1 & 0 & 0 & 0 \\ -1 & 2 & -2 + 2\varepsilon & -1 & 1 - 2\varepsilon & 0 \\ 1 & -2 + 2\varepsilon & 2 & 1 - 2\varepsilon & -1 & 0 \\ 0 & -1 & 1 - 2\varepsilon & 2 & 2\varepsilon & -1 \\ 0 & 1 - 2\varepsilon & -1 & 2\varepsilon & 2 & -1 \\ 0 & 0 & 0 & -1 & -1 & 1 \end{array} \right] \left\{ \begin{array}{l} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{array} \right\} \qquad (36)$$

With the domain formulation of the second column of Table 1, the discrete matrix is:

$$\left[ \begin{array}{cccccc} 1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 2 - \varepsilon & 0 & -1 + \varepsilon & 0 & 0 \\ 0 & 0 & \varepsilon & 0 & -\varepsilon & 0 \\ 0 & -1 + \varepsilon & 0 & 1 - \varepsilon & 0 & 0 \\ 0 & 0 & -\varepsilon & 0 & 1 + \varepsilon & -1 \\ 0 & 0 & 0 & 0 & -1 & 1 \end{array} \right] \left\{ \begin{array}{l} \alpha_1 \\ \alpha_{2,1} \\ \alpha_{2,2} \\ \alpha_{3,1} \\ \alpha_{3,2} \\ \alpha_4 \end{array} \right\} \qquad (37)$$
With the shifted formulation of the third column of Table 1, the discrete matrix is:

$$ \begin{bmatrix} 1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 2 & 2\epsilon & -1 & 2-2\epsilon & 0 \\ 0 & 2\epsilon & 4\epsilon & -2\epsilon & 0 & 0 \\ 0 & -1 & -2\epsilon & 2 & -2+2\epsilon & -1 \\ 0 & 2-2\epsilon & 0 & -2+2\epsilon & 4-4\epsilon & 0 \\ 0 & 0 & 0 & -1 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ d_2 \\ c_3 \\ d_3 \\ c_4 \end{bmatrix} \quad (38) $$

We enforce boundary conditions directly on the discrete problem, as equality constraints through Lagrange multipliers. The position of the interface depends only on $\epsilon$. Nonetheless, the solution does not depend on the $\epsilon$ parameter, so that the relative error does not depend on $\epsilon$, which makes our analysis relevant.

Let us consider the following arbitrary numerical values: $\alpha = 10/9$ and $\beta = 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¹-norm:

$$ \|u - u_h\|_{H^1} = \sqrt{\int_0^3 \left[ (u - u_h)^2 + (u' - u_h')^2 \right] dx} \quad (39) $$

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¹-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 \cdot 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:

$$ Ku = f \rightarrow K'u' = f' \text{ with } \begin{cases} K' = D_c K D_c \\ u = D_c u' \\ f' = D_c f \end{cases} \quad (40) $$

In linear elasticity, $D_c$ is usually defined as,

$$ [D_c]_{i,i} = \frac{1}{\sqrt{K_{i,i}}} \sqrt{\frac{\max(K_{i,i}) + \min(K_{i,i})}{2}} \quad (41) $$

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\Phi_2 \approx -\Phi_2$). The difference between those equations lies in the terms $1 - 2\epsilon$ and $2 - 2\epsilon$. Those terms imply the sum of heterogeneous quantities that leads to a tremendous loss of accuracy on the difference information $2\epsilon$. For instance, let us consider the following string of calculations in “double precision” arithmetic (8-octet storage):

$$ 1 - (1 - 6.3847497084E-12) \approx 6.348E-12 \quad (42) $$
With the shifted formulation of the third column of Table 1, the discrete matrix is:

$$ \begin{bmatrix} 1 & -1 & 0 & 0 & 0 & 0 \\ -1 & 2 & 2\epsilon & -1 & 2-2\epsilon & 0 \\ 0 & 2\epsilon & 4\epsilon & -2\epsilon & 0 & 0 \\ 0 & -1 & -2\epsilon & 2 & -2+2\epsilon & -1 \\ 0 & 2-2\epsilon & 0 & -2+2\epsilon & 4-4\epsilon & 0 \\ 0 & 0 & 0 & -1 & 0 & 1 \end{bmatrix} \begin{bmatrix} c_1 \\ c_2 \\ d_2 \\ c_3 \\ d_3 \\ c_4 \end{bmatrix} \quad (38) $$

We enforce boundary conditions directly on the discrete problem, as equality constraints through Lagrange multipliers. The position of the interface depends only on $\epsilon$. Nonetheless, the solution does not depend on the $\epsilon$ parameter, so that the relative error does not depend on $\epsilon$, which makes our analysis relevant.

Let us consider the following arbitrary numerical values: $\alpha = 10/9$ and $\beta = 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¹-norm:

$$ \|u - u_h\|_{H^1} = \sqrt{\int_0^3 \left[ (u - u_h)^2 + (u' - u_h')^2 \right] dx} \quad (39) $$

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¹-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 \cdot 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:

$$ Ku = f \rightarrow K'u' = f' \text{ with } \begin{cases} K' = D_c K D_c \\ u = D_c u' \\ f' = D_c f \end{cases} \quad (40) $$

In linear elasticity, $D_c$ is usually defined as,

$$ [D_c]_{i,i} = \frac{1}{\sqrt{K_{i,i}}} \sqrt{\frac{\max(K_{i,i}) + \min(K_{i,i})}{2}} \quad (41) $$

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\Phi_2 \approx -\Phi_2$). The difference between those equations lies in the terms $1 - 2\epsilon$ and $2 - 2\epsilon$. Those terms imply the sum of heterogeneous quantities that leads to a tremendous loss of accuracy on the difference information $2\epsilon$. For instance, let us consider the following string of calculations in “double precision” arithmetic (8-octet storage):

$$ 1 - (1 - 6.3847497084E-12) \approx 6.348E-12 \quad (42) $$
<figure><img src="image_13.png" /><figcaption>**Fig. 9** H1-error of strongly discontinuous approximations with linear elements. Shifted and domain enrichments outperform X-FEM formulation when "eps" (ε) goes to zero</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Conditioning of strongly discontinuous approximations without the use of pre-conditioners. The three approximation spaces yield the same condition number (as they are equivalent)</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. 9** H1-error of strongly discontinuous approximations with linear elements. Shifted and domain enrichments outperform X-FEM formulation when "eps" (ε) goes to zero</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 10** Conditioning of strongly discontinuous approximations without the use of pre-conditioners. The three approximation spaces yield the same condition number (as they are equivalent)</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. 11** 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</figcaption></figure>

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.
<figure><img src="image_15.png" /><figcaption>**Fig. 11** 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</figcaption></figure>

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.
<figure><img src="image_16.png" /><figcaption>**Fig. 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</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**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</figcaption></figure>

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*(Φ<sub>S</sub>)} → 0 for a given vertex node, the middle node satisfies likewise *measure*({Ω⁺ ∩ *Supp*(Φ<sub>M</sub>)} → 0, on the vertex patch. For such a configuration, the enrichment shape functions χ<sub>Ω⁺</sub>Φ<sub>S</sub> and χ<sub>Ω⁺</sub>Φ<sub>M</sub> become almost homothetic (Fig. 16).
<figure><img src="image_16.png" /><figcaption>**Fig. 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</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**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</figcaption></figure>

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(\{\Omega_+ \cap Supp(\Phi_S)\}) \to 0$ for a given vertex node, the middle node satisfies likewise $measure(\{\Omega_+ \cap Supp(\Phi_M)\}) \to 0$, on the vertex patch. For such a configuration, the enrichment shape functions $\chi_{\Omega_+}\Phi_S$ and $\chi_{\Omega_+}\Phi_M$ become almost homothetic (Fig. 16).
<figure><img src="image_18.png" /><figcaption>Fig. 14 H1-error and conditioning analysis with quadratic elements (Lagrange). The conditioning slopes are about three as shown is “Numerical analysis of strongly discontinuous approximations with quadratic elements” section</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15 The use of linear element pre-conditioners barely improves the accuracy and condition number of quadratic elements (Lagrange)</figcaption></figure>

* In case of Lagrange polynomials, we have:

$$ 
\begin{aligned}
\Phi_S(\xi) &= (1 - \xi)(1 - 2\xi), \quad \forall \xi \in [0, 1] \\
\Phi_M(\xi) &= 4\xi(1 - \xi), \quad \forall \xi \in [0, 1]
\end{aligned}
\qquad (43)
$$

Both functions can be related with the following expression,

$$ \Phi_S(\xi) = -4\Phi_M(\xi) + 4(1 - \xi)^2 \qquad (44) $$
<figure><img src="image_18.png" /><figcaption>Fig. 14 H1-error and conditioning analysis with quadratic elements (Lagrange). The conditioning slopes are about three as shown is “Numerical analysis of strongly discontinuous approximations with quadratic elements” section</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 15 The use of linear element pre-conditioners barely improves the accuracy and condition number of quadratic elements (Lagrange)</figcaption></figure>

* In case of Lagrange polynomials, we have:

$$ 
\begin{aligned}
\Phi_S(\xi) &= (1 - \xi)(1 - 2\xi), \quad \forall \xi \in [0, 1] \\
\Phi_M(\xi) &= 4\xi(1 - \xi), \quad \forall \xi \in [0, 1]
\end{aligned}
\qquad (43) $$

Both functions can be related with the following expression,

$$ \Phi_S(\xi) = -4\Phi_M(\xi) + 4(1 - \xi)^2 \qquad (44) $$
<figure><img src="image_20.png" /><figcaption>Fig. 16 Coupling between Lagrange quadratic shape functions. When the crack passes close to the nodes, the discontinuous shape functions of the middle node and one vertex node become almost colinear</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 17 Graphs of polynomial functions $e_1$ (a) and $e_2$ (b)</figcaption></figure>

when $\xi$ goes to 1, the difference information between those functions decreases accordingly to $(1 - \xi)^2$. This generates a vanishing sub-space (consisting of polynomials only defined on the small support $\{\Omega_+ \cap \text{Supp}(\Phi_S)\}$), which implies a conditioning slope of at least 3 with quadratic shape functions (Fig. 15). In other words, the condition number $\kappa(K)$ is bounded by:

$$ \kappa(K) \geq C\epsilon^{-3} \qquad (45) $$

*Proof* Let $\kappa(K)$ be the condition number of the stiffness matrix:

$$ \kappa(K) = \frac{\max(\{e^T K e\}_{\|e\|_2=1})}{\min(\{e^T K e\}_{\|e\|_2=1})} \geq \frac{e_1^T K e_1}{e_2^T K e_2} \quad \forall e_1, e_2/\|e_1\|_2 = 1, \|e_2\|_2 = 1 \qquad (46) $$

where $\|\cdot\|_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:
<figure><img src="image_20.png" /><figcaption>Fig. 16 Coupling between Lagrange quadratic shape functions. When the crack passes close to the nodes, the discontinuous shape functions of the middle node and one vertex node become almost colinear</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 17 Graphs of polynomial functions e₁ (a) and e₂ (b)</figcaption></figure>

when $\xi$ goes to 1, the difference information between those functions decreases accordingly to $(1 - \xi)^2$. This generates a vanishing sub-space (consisting of polynomials only defined on the small support $\{\Omega_+ \cap \text{Supp}(\Phi_S)\}$), which implies a conditioning slope of at least 3 with quadratic shape functions (Fig. 15). In other words, the condition number $\kappa(K)$ is bounded by:

$$ \kappa(K) \geq C\epsilon^{-3} \qquad (45) $$

*Proof* Let $\kappa(K)$ be the condition number of the stiffness matrix:

$$ \kappa(K) = \frac{\max(\{e^T K e\}_{\|e\|_2=1})}{\min(\{e^T K e\}_{\|e\|_2=1})} \geq \frac{e_1^T K e_1}{e_2^T K e_2} \quad \forall e_1, e_2/\|e_1\|_2 = 1, \|e_2\|_2 = 1 \qquad (46) $$

where $\|\cdot\|_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:
$$ egin{cases} e_1 = -\frac{1}{\sqrt{5}}(2\Phi_1 + \Phi_2) \\ e_2 = \frac{1}{\sqrt{17}}(4\chi_{\Omega} - \Phi_3 + \chi_{\Omega} - \Phi_4) \end{cases} \qquad (47) $$

In the discrete space of [33], we have:

$$ \|u\|_2 = \sqrt{\sum_{i \in I} a_i^2 + \sum_{j \in I_H} \alpha_{j,1}^2 + \sum_{j \in I_H} \alpha_{j,2}^2} \qquad (48) $$

so that we have $\|e_1\|_2 = \|e_2\|_2 = 1$.

As $\frac{1}{2}e_1^T Ke_1$ and $\frac{1}{2}e_2^T Ke_2$ represent the elastic energy of these solutions, we have:

$$ \frac{e_1^T Ke_1}{e_2^T Ke_2} = \frac{\int_0^1 (e_1')^2 dx}{\int_{2-\epsilon}^2 (e_2')^2 dx} = \frac{51}{80\epsilon^3} $$

Thus,

$$ \kappa(K) \geq C\epsilon^{-3} $$

The asymptotic behavior of $\Phi_M(\xi)$ [at the first order of $(1 - \xi)$] is:

$$ \Phi_M(\xi) \approx -\Phi_S(\xi)/4 $$

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:

$$ \begin{aligned} \Phi_S(\xi) &= (1 - \xi)^2, \quad \forall \xi \in [0, 1] \\ \Phi_M(\xi) &= 2\xi(1 - \xi), \quad \forall \xi \in [0, 1] \end{aligned} \qquad (49) $$

when $\xi$ goes to 1, there is no coupling between Bernstein polynomials as observed with Lagrange polynomials. However, the polynomial $\Phi_S$ cancels out accordingly to $(1 - \xi)^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 behavior of singular approximations at the crack tip

### 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,
$$ \begin{cases} e_1 = -\frac{1}{\sqrt{5}}(2\Phi_1 + \Phi_2) \\ e_2 = \frac{1}{\sqrt{17}}(4\chi_{\Omega} - \Phi_3 + \chi_{\Omega} - \Phi_4) \end{cases} \qquad (47) $$

In the discrete space of [33], we have:

$$ \|u\|_2 = \sqrt{\sum_{i \in I} a_i^2 + \sum_{j \in I_H} \alpha_{j,1}^2 + \sum_{j \in I_H} \alpha_{j,2}^2} \qquad (48) $$

so that we have $\|e_1\|_2 = \|e_2\|_2 = 1$.

As $\frac{1}{2}e_1^T Ke_1$ and $\frac{1}{2}e_2^T Ke_2$ represent the elastic energy of these solutions, we have:

$$ \frac{e_1^T Ke_1}{e_2^T Ke_2} = \frac{\int_0^1 (e_1')^2 dx}{\int_{2-\epsilon}^2 (e_2')^2 dx} = \frac{51}{80\epsilon^3} $$

Thus,

$$ \kappa(K) \geq C\epsilon^{-3} $$

The asymptotic behavior of $\Phi_M(\xi)$ [at the first order of $(1 - \xi)$] is:

$$ \Phi_M(\xi) \approx -\Phi_S(\xi)/4 $$

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:

$$ \begin{aligned} \Phi_S(\xi) &= (1 - \xi)^2, \quad \forall \xi \in [0, 1] \\ \Phi_M(\xi) &= 2\xi(1 - \xi), \quad \forall \xi \in [0, 1] \end{aligned} \qquad (49) $$

when $\xi$ goes to 1, there is no coupling between Bernstein polynomials as observed with Lagrange polynomials. However, the polynomial $\Phi_S$ cancels out accordingly to $(1 - \xi)^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 behavior of singular approximations at the crack tip

### 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,
<figure><img src="image_22.png" /><figcaption>**Fig. 18** Asymptotic behavior with Bernstein polynomials. Bernstein polynomials allow emphasizing once more that conditioning and accuracy are separated issues: the condition numbers of all three approximation spaces are almost the same, but there is a great difference in accuracy when "eps" (ε) goes to zero</figcaption></figure>

- Cut-off enrichment is not very convenient, as assembling global d.o.f. is out of the scope of the industrial software we used (Code_Aster). Moreover, the method does not extend properly in 3D [29].

*Example 1* (Horizontal crack opening in mode I with linear elements). The domain geometry is a square defined by $\Omega = [-0.5, +0.5] \times [-0.5, +0.5]$. The meshing procedure subdivides the domain into regular sized cells, which edges are parallel to the crack. The size of the enrichment zone is $r = 0.1$.

The analytical solution in displacement corresponds to the exact mode I limited to the $\sqrt{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:

$$ \|u - u^h\|_{\text{energy}} = \sqrt{\int_{\Omega} (\sigma^h - \sigma) : (\varepsilon^h - \varepsilon)d\Omega / \int_{\Omega} \sigma : \varepsilon d\Omega} \quad (50) $$
<figure><img src="image_22.png" /><figcaption>Fig. 18 Asymptotic behavior with Bernstein polynomials. Bernstein polynomials allow emphasizing once more that conditioning and accuracy are separated issues: the condition numbers of all three approximation spaces are almost the same, but there is a great difference in accuracy when "eps" (ε) goes to zero</figcaption></figure>

- Cut-off enrichment is not very convenient, as assembling global d.o.f. is out of the scope of the industrial software we used (Code_Aster). Moreover, the method does not extend properly in 3D [29].

Example 1 (Horizontal crack opening in mode I with linear elements). The domain geometry is a square defined by $\Omega = [-0.5, +0.5] \times [-0.5, +0.5]$. The meshing procedure subdivides the domain into regular sized cells, which edges are parallel to the crack. The size of the enrichment zone is $r = 0.1$.

The analytical solution in displacement corresponds to the exact mode I limited to the $\sqrt{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:

$$ \|u - u^h\|_{\text{energy}} = \sqrt{\int_{\Omega} (\sigma^h - \sigma) : (\varepsilon^h - \varepsilon)d\Omega / \int_{\Omega} \sigma : \varepsilon d\Omega} \quad (50) $$
<figure><img src="image_23.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 20** Boundary conditions around the patch-test. Mixed boundary conditions are applied to avoid enforcing displacement on the nodes of the boundary near the crack interface which would add only a technical difficulty beyond the scope of the paper</figcaption></figure>

On 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 $\sqrt{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 Dirichlet boundary conditions
<figure><img src="image_23.png" /><figcaption>**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</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 20** Boundary conditions around the patch-test. Mixed boundary conditions are applied to avoid enforcing displacement on the nodes of the boundary near the crack interface which would add only a technical difficulty beyond the scope of the paper</figcaption></figure>

On 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 $\sqrt{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 Dirichlet boundary conditions
<figure><img src="image_25.png" /><figcaption>**Fig. 21** Conditioning of crack approximations. The condition number increases swiftly with the geometrical enrichment strategy i.e. the introduction of singular functions on many layers of elements around the crack-tip</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 22** Convergence analysis of crack approximation spaces. The single-layer enrichment (topological) convergence curves are added to recall the benefit of the geometrical enrichment strategy</figcaption></figure>

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).
<figure><img src="image_25.png" /><figcaption>**Fig. 21** Conditioning of crack approximations. The condition number increases swiftly with the geometrical enrichment strategy i.e. the introduction of singular functions on many layers of elements around the crack-tip</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 22** Convergence analysis of crack approximation spaces. The single-layer enrichment (topological) convergence curves are added to recall the benefit of the geometrical enrichment strategy</figcaption></figure>

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).
<figure><img src="image_27.png" /><figcaption>**Fig. 23** Analytical solution corresponding to a mode 1 displacement field for the inclined crack</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 24** Boundary conditions for an inclined crack. The boundary conditions are the same with the horizontal crack test</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 25** Conditioning of vectorial and scalar enrichments with a crack orientation at 44.9°</figcaption></figure>
<figure><img src="image_27.png" /><figcaption>**Fig. 23** Analytical solution corresponding to a mode 1 displacement field for the inclined crack</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Fig. 24** Boundary conditions for an inclined crack. The boundary conditions are the same with the horizontal crack test</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 25** Conditioning of vectorial and scalar enrichments with a crack orientation at 44.9°</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Fig. 26** Convergence analysis for vectorial and scalar enrichments with a crack orientation at 44.9°. Both enrichments have optimal convergence order but error levels are quite high. Results are in accordance with the ones of Fig. 22 and X-FEM is more accurate than GFEM</figcaption></figure>

**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-
<figure><img src="image_30.png" /><figcaption>**Fig. 26** Convergence analysis for vectorial and scalar enrichments with a crack orientation at 44.9°. Both enrichments have optimal convergence order but error levels are quite high. Results are in accordance with the ones of Fig. 22 and X-FEM is more accurate than GFEM</figcaption></figure>

**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-
<figure><img src="image_31.png" /><figcaption>**Fig. 27** Numerical behavior of the X-FEM vectorial enrichment with quadratic elements</figcaption></figure>

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.

## Improvement of strongly discontinuous approximations

### 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. $\chi_{\Omega+}\Phi_S$ and $\chi_{\Omega+}\Phi_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. $\chi_{\Omega+}\Phi_M$ for $\varepsilon$ 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 $\varepsilon = 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. ($\chi_{\Omega+}\Phi_S$ and $\chi_{\Omega+}\Phi_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:
<figure><img src="image_31.png" /><figcaption>**Fig. 27** Numerical behavior of the X-FEM vectorial enrichment with quadratic elements</figcaption></figure>

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.

## Improvement of strongly discontinuous approximations

### 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. $\chi_{\Omega+}\Phi_S$ and $\chi_{\Omega+}\Phi_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. $\chi_{\Omega+}\Phi_M$ for $\varepsilon$ 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 $\varepsilon = 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. ($\chi_{\Omega+}\Phi_S$ and $\chi_{\Omega+}\Phi_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:
<figure><img src="image_32.png" /><figcaption>Fig. 28 H1-error with quadratic Lagrange polynomials. The enrichments are alternated through the removal of the middle node beyond a threshold eps of ε = 10⁻³</figcaption></figure>

* 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.

### Extrapolation in 3D

Let us consider the patch-test of [24]. It is a 3D block of side 4m defined by $\Omega = [-2, +2] \times [-2, +2] \times [-2, +2]$ and split by an inclined interface. The domain is meshed with regular quadratic hexahedral elements (4 × 4 × 4 elements).

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 $\delta$ as:

$$x + y + z + \delta = 0 \qquad (51)$$

The error is analyzed with the energy norm:

$$\|u - u^h\|_{\text{energy}} = \sqrt{\int_{\Omega} (\sigma^h - \sigma) : (\varepsilon^h - \varepsilon) d\Omega / \int_{\Omega} \sigma : \varepsilon d\Omega} \qquad (52)$$

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),
<figure><img src="image_32.png" /><figcaption>Fig. 28 H1-error with quadratic Lagrange polynomials. The enrichments are alternated through the removal of the middle node beyond a threshold eps of ε = 10⁻³</figcaption></figure>

* 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.

### Extrapolation in 3D

Let us consider the patch-test of [24]. It is a 3D block of side 4m defined by $\Omega = [-2, +2] \times [-2, +2] \times [-2, +2]$ and split by an inclined interface. The domain is meshed with regular quadratic hexahedral elements (4 × 4 × 4 elements).

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 $\delta$ as:

$$x + y + z + \delta = 0 \qquad (51)$$

The error is analyzed with the energy norm:

$$\|u - u^h\|_{\text{energy}} = \sqrt{\int_{\Omega} (\sigma^h - \sigma) : (\varepsilon^h - \varepsilon) d\Omega / \int_{\Omega} \sigma : \varepsilon d\Omega} \qquad (52)$$

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),
<figure><img src="image_33.png" /><figcaption>**Fig. 29** 3D test-case with quadratic elements. The strategies to improve the condition number are assessed on a 3D patch-test</figcaption></figure>

* 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]:

$$ \min \left( \begin{array}{l} \text{measure}(\{\Omega_+ \cap \text{Supp}(\Phi_M)\}), \\ \text{measure}(\{\Omega_- \cap \text{Supp}(\Phi_M)\}) \end{array} \right) \le 10^{-9} \text{measure}(\{\text{Supp}(\Phi_M)\}) \quad (53) $$

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 $\delta/\sqrt{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

<figure><img src="image_34.png" /><figcaption>**Fig. 30** Behavior of strongly discontinuous enrichments in 3D with the corrections studied previously</figcaption></figure>
<figure><img src="image_33.png" /><figcaption>Fig. 29 3D test-case with quadratic elements. The strategies to improve the condition number are assessed on a 3D patch-test</figcaption></figure>

* 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]:

$$ \min \left( \begin{array}{l} \text{measure}(\{\Omega_+ \cap \text{Supp}(\Phi_M)\}), \\ \text{measure}(\{\Omega_- \cap \text{Supp}(\Phi_M)\}) \end{array} \right) \le 10^{-9} \text{measure}(\{\text{Supp}(\Phi_M)\}) \quad (53) $$

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 $\delta/\sqrt{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

<figure><img src="image_34.png" /><figcaption>Fig. 30 Behavior of strongly discontinuous enrichments in 3D with the corrections studied previously</figcaption></figure>
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 com-
promise, 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.

Singular approximation at a crack tip

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_\alpha$ 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 $\Omega_+$ and the one on $\Omega_-$ to the whole domain $\Omega$ (an example is given Fig. 31). Hence, two continuous extensions of $K_\alpha(r, \theta)$ are introduced:

* $K_\alpha(r, |\theta|)$ overlaps with $K_\alpha(r, \theta)$ on $\Omega_+$,

* $K_\alpha(r, -|\theta|)$ overlaps with $K_\alpha(r, \theta)$ on $\Omega_-$.

In the elements where $K_\alpha$ has a singular point, the interpolation of $K_\alpha$ is smoothened, so that subsequent interpolation errors of the singular function $K_\alpha$ do not pollute the accuracy of the displacement field as observed in [25].

<figure><img src="image_35.png" /><figcaption>**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)</figcaption></figure>
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 com-
promise, 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.

Singular approximation at a crack tip

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_\alpha$ 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 $\Omega_+$ and the one on $\Omega_-$ to the whole domain $\Omega$ (an example is given Fig. 31). Hence, two continuous extensions of $K_\alpha(r, \theta)$ are introduced:

* $K_\alpha(r, |\theta|)$ overlaps with $K_\alpha(r, \theta)$ on $\Omega_+$,

* $K_\alpha(r, -|\theta|)$ overlaps with $K_\alpha(r, \theta)$ on $\Omega_-$.

In the elements where $K_\alpha$ has a singular point, the interpolation of $K_\alpha$ is smoothened, so that subsequent interpolation errors of the singular function $K_\alpha$ do not pollute the accuracy of the displacement field as observed in [25].

<figure><img src="image_35.png" /><figcaption>**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)</figcaption></figure>
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 $\underline{K}_{\alpha}$, so that the interpolation of $\underline{K}_{\alpha}$ is smooth when the polar coordinate $r$ goes to zero. Let $m$ be the order of the interpolation of $\underline{K}_{\alpha}$, which can be different from $n$. The new interpolation operator reads:

$$ \tilde{\Pi}^m \underline{K}_{\alpha} = \sum_{k \in \{I_{CT,m} \cup \delta I_{CT,m} \} / I_H} \underline{K}_{\alpha}(r_k, \theta_k) \Gamma_{k,m} + \sum_{k \in I_H \cap \{I_{CT,m} \cup \delta I_{CT,m} \} / I_T} \left[ \underline{K}_{\alpha}(r_k, |\theta_k|) \Gamma_{k,m} \chi_{\Omega^+} + \underline{K}_{\alpha}(r_k, -|\theta_k|) \Gamma_{k,m} \chi_{\Omega^-} \right] \quad (54) $$

where $\Gamma_{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} \cup \delta 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 $\Omega_+$,

$$ \tilde{\Pi}^m(\underline{K}_{\alpha}) = \sum_{k \in \{I_{CT,m} \cup \delta I_{CT,m}\} / I_T} \underline{K}_{\alpha}(r_k, |\theta_k|) \Gamma_{k,m} \quad (55) $$

* Assuming the evaluation point (gauss point) is located on $\Omega_-$,

$$ \tilde{\Pi}^m(\underline{K}_{\alpha}) = \sum_{k \in \{I_{CT,m} \cup \delta I_{CT,m}\} / I_T} \underline{K}_{\alpha}(r_k, -|\theta_k|) \Gamma_{k,m} \quad (56) $$

The new “bubble” approximation space, based on the new interpolation operator, is:

$$ V_{bub,m,n}^h = \left\{ \begin{array}{l} v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} (H - H(x_i)) \Phi_i E_j \\ \qquad + \sum_{i \in I_{CT,n}} \sum_{\alpha=1}^{d} c_{i,\alpha} \Gamma_{i,n} (\underline{K}_{\alpha} - \tilde{\Pi}^m \underline{K}_{\alpha}) \end{array} \right\} \quad (57) $$

where $\Phi_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 $\Gamma_{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_{bub,m,n}^h$ 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
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 $\underline{K}_{\alpha}$, so that the interpolation of $\underline{K}_{\alpha}$ is smooth when the polar coordinate $r$ goes to zero. Let $m$ be the order of the interpolation of $\underline{K}_{\alpha}$, which can be different from $n$. The new interpolation operator reads:

$$ \tilde{\Pi}^m \underline{K}_{\alpha} = \sum_{k \in \{I_{CT,m} \cup \delta I_{CT,m} \} / I_H} \underline{K}_{\alpha}(r_k, \theta_k) \Gamma_{k,m} + \sum_{k \in I_H \cap \{I_{CT,m} \cup \delta I_{CT,m} \} / I_T} \left[ \underline{K}_{\alpha}(r_k, |\theta_k|) \Gamma_{k,m} \chi_{\Omega^+} + \underline{K}_{\alpha}(r_k, -|\theta_k|) \Gamma_{k,m} \chi_{\Omega^-} \right] \quad (54) $$

where $\Gamma_{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} \cup \delta 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 $\Omega_+$,

$$ \tilde{\Pi}^m(\underline{K}_{\alpha}) = \sum_{k \in \{I_{CT,m} \cup \delta I_{CT,m}\} / I_T} \underline{K}_{\alpha}(r_k, |\theta_k|) \Gamma_{k,m} \quad (55) $$

* Assuming the evaluation point (gauss point) is located on $\Omega_-$,

$$ \tilde{\Pi}^m(\underline{K}_{\alpha}) = \sum_{k \in \{I_{CT,m} \cup \delta I_{CT,m}\} / I_T} \underline{K}_{\alpha}(r_k, -|\theta_k|) \Gamma_{k,m} \quad (56) $$

The new “bubble” approximation space, based on the new interpolation operator, is:

$$ V_{bub,m,n}^h = \left\{ \begin{array}{l} v^h = \sum_{i \in I} \sum_{j=1}^{d} a_{i,j} \Phi_i E_j + \sum_{i \in I_H} \sum_{j=1}^{d} b_{i,j} (H - H(x_i)) \Phi_i E_j \\ + \sum_{i \in I_{CT,n}} \sum_{\alpha=1}^{d} c_{i,\alpha} \Gamma_{i,n} (\underline{K}_{\alpha} - \tilde{\Pi}^m \underline{K}_{\alpha}) \end{array} \right\} \quad (57) $$

where $\Phi_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 $\Gamma_{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_{bub,m,n}^h$ 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 func-
tion 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_{bub,m=2,n=1}^h$ and $V_{bub,m=2,n=2}^h$ are about five times more accurate than the bubble space $V_{bub,m=1,n=1}^h$. The bubble space $V_{bub,m=1,n=1}^h$ has the same accuracy than X-FEM vectorial enrichment. Clearly, the bubble spaces with

<figure><img src="image_36.png" /><figcaption>**Fig. 32** Conditioning of “bubble” enrichment vs X-FEM vectorial enrichment</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. 33** convergence analysis of “bubble” enrichment strategies with vectorial enrichment</figcaption></figure>
continuous extensions in the bandwidth of elements where the enrichment func-
tion 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_{bub,m=2,n=1}^h$ and $V_{bub,m=2,n=2}^h$ are about five times more accurate than the bubble space $V_{bub,m=1,n=1}^h$. The bubble space $V_{bub,m=1,n=1}^h$ has the same accuracy than X-FEM vectorial enrichment. Clearly, the bubble spaces with

<figure><img src="image_36.png" /><figcaption>**Fig. 32** Conditioning of “bubble” enrichment vs X-FEM vectorial enrichment</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Fig. 33** convergence analysis of “bubble” enrichment strategies with vectorial enrichment</figcaption></figure>
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_{bub, m=2, n=1}^h$ and $V_{bub, m=2, n=2}^h$, 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_{bub, m=1, n=2}^h$ 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].

• The good results detailed in this section with quadratic elements, are confirmed with linear elements. In this case the condition number for the bubble space $V_{bub, m=1, n=1}^h$ evolves with an order with respect to mesh refinement around 2 which corresponds to the optimal expected value (see Fig. 34 related to Figs. 21, 22 and Fig. 35 related to Figs. 25, 26 of "Numerical analysis of crack approximations with linear elements: crack opening in mode 1" section). 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).

<figure><img src="image_38.png" /><figcaption>Fig. 34 Convergence and conditioning analysis with linear elements for a horizontal crack</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. 35 Convergence and conditioning analysis with linear elements for inclined crack at an angle of 44.9°</figcaption></figure>
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_{bub, m=2, n=1}^h$ and $V_{bub, m=2, n=2}^h$, 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_{bub, m=1, n=2}^h$ 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].

• The good results detailed in this section with quadratic elements, are confirmed with linear elements. In this case the condition number for the bubble space $V_{bub, m=1, n=1}^h$ evolves with an order with respect to mesh refinement around 2 which corresponds to the optimal expected value (see Fig. 34 related to Figs. 21, 22 and Fig. 35 related to Figs. 25, 26 of "Numerical analysis of crack approximations with linear elements: crack opening in mode 1" section). 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).

<figure><img src="image_38.png" /><figcaption>Fig. 34 Convergence and conditioning analysis with linear elements for a horizontal crack</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. 35 Convergence and conditioning analysis with linear elements for inclined crack at an angle of 44.9°</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>**Fig. 36** Convergence analysis of stress intensity factor for the mode I problem. Stress Intensity Factors (SIF) 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</figcaption></figure>

## 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:

$$u^h(\mathbf{x}) = \sum_{i \in I} a_i \Phi_i(\mathbf{x}) + \sum_{j \in I_H} b_j F(\mathbf{x}) \Phi_j(\mathbf{x}) \qquad (58)$$

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(\mathbf{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].

<table><caption>**Table 3** Weak discontinuity enrichments comparison</caption><thead><tr><th></th><th>Sukumar [14]</th><th>Moës [27]</th><th>Belytschko [28]</th></tr></thead><tbody><tr><td>F(x)</td><td>Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x)</td><td>Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x) - Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x)</td><td>Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x) - |lsn<sub>i</sub>|</td></tr></tbody></table>
<figure><img src="image_40.png" /><figcaption>**Fig. 36** Convergence analysis of stress intensity factor for the mode I problem. Stress Intensity Factors (SIF) 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</figcaption></figure>

## 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:

$$u^h(\mathbf{x}) = \sum_{i \in I} a_i \Phi_i(\mathbf{x}) + \sum_{j \in I_H} b_j F(\mathbf{x}) \Phi_j(\mathbf{x}) \quad (58)$$

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(\mathbf{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].

<table><caption>**Table 3** Weak discontinuity enrichments comparison</caption><thead><tr><th></th><th>Sukumar [14]</th><th>Moës [27]</th><th>Belytschko [28]</th></tr></thead><tbody><tr><td>F(x)</td><td>Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x)</td><td>Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x) - Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x)</td><td>Σ<sub>i∈I<sub>H</sub></sub> |lsn<sub>i</sub>| Φ<sub>i</sub>(x) - |lsn<sub>i</sub>|</td></tr></tbody></table>
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 Belytschko'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 \in I_H$

$$ \left\{ \begin{array}{l} a'_j = a_j - |\ln s_j| b_j \\ b'_j = b_j \end{array} \right. \qquad (59) $$

Thus, in the next sections both approximations may display similar numerical properties.

Note also that the zero of the level set term $\sum_{i \in I_H} \ln s_i \Phi_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

$$ \frac{d^2 u}{dx^2} = 0, \quad \forall x \in \Omega = [0, x_0[ \cup [x_0, L], \text{ so that,} \\ u(x) = \begin{cases} \alpha^{-}x + \beta, & \text{if } x \le x_0 \\ \alpha^{+}x + \beta', & \text{if } x \ge x_0 \end{cases}, \text{ where } x_0 \text{ 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,

$$ u(0) = 0, \quad u(L) = \bar{u} $$

* Material properties are given by:

$$ E_i(x) = \begin{cases} E_1, & \text{if } x \le x_0 \\ E_2, & \text{if } x \ge x_0 \end{cases}, \text{ where } x_0 \text{ is the location of the interface,} $$
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 Belytschko'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 \in I_H$

$$ \left\{ \begin{array}{l} a'_j = a_j - |\operatorname{lsln}_j| b_j \\ b'_j = b_j \end{array} \right. \qquad (59) $$

Thus, in the next sections both approximations may display similar numerical properties.

Note also that the zero of the level set term $\sum_{i \in I_H} lsln_i \Phi_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

$$ \begin{gathered} \frac{d^2 u}{dx^2} = 0, \quad \forall x \in \Omega = [0, x_0[ \cup [x_0, L], \text{ so that,} \\ u(x) = \begin{cases} \alpha^{-}x + \beta, & \text{if } x \le x_0 \\ \alpha^{+}x + \beta', & \text{if } x \ge x_0 \end{cases}, \text{ where } x_0 \text{ is the location of the interface,} \end{gathered} $$

* With arbitrary Dirichlet boundary conditions we choose to impose on the central element nodes for the sake of simplicity,

$$ u(0) = 0, \quad u(L) = \bar{u} $$

* Material properties are given by:

$$ E_i(x) = \begin{cases} E_1, & \text{if } x \le x_0 \\ E_2, & \text{if } x \ge x_0 \end{cases}, \text{ where } x_0 \text{ is the location of the interface,} $$
* Assuming small strains, the expected solution satisfying the boundary conditions above is:

$$u(x) = \begin{cases} \frac{E_2 \bar{u}}{E_1 L + (E_2 - E_1) x_0} x, & \text{if } x \le x_0 \quad \text{so that } \alpha^- = \frac{E_2 \bar{u}}{E_1 L + (E_2 - E_1) x_0} \\ \frac{E_1 x + (E_2 - E_1) x_0}{E_1 L + (E_2 - E_1) x_0} \bar{u}, & \text{if } x \ge x_0 \quad \text{so that } \alpha^+ = \frac{E_1 \bar{u}}{E_1 L + (E_2 - E_1) x_0} \end{cases} \quad (60)$$

## 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 $\Omega = [-1, 2]$ with the following Dirichlet boundary conditions $u(-1) = 0$, $u(2) = \bar{u}$ so that $L = 3$ m. The interface is $\Gamma = \{x_0\}$.

The weak form of the problem becomes:

Find:

$$\left. \begin{array}{l} u \in \prod : \{ w \in H^1(\Omega) / w(-1) = 0, w(2) = \bar{u}, w(x_0^-) = w(x_0^+) \} \\ \forall v \in \prod_0 : \{ w \in H^1(\Omega) / w(-1) = 0, w(2) = 0, w(x_0^-) = w(x_0^+) \} \end{array} \right\} \\ \mapsto \begin{cases} \int_{-1}^{2} E_i u' v' dx = 0 \\ E_1 u'(x_0^-) = E_2 u'(x_0^+) \end{cases} \quad (61)$$

The X-FEM discrete linear space is:

$$\prod_h : \left\{ \begin{array}{l} w/w = a_1 \Phi_1 + a_2 \Phi_2 + b_2 F \Phi_2 + a_3 \Phi_3 + b_3 F \Phi_3 + a_4 \Phi_4 \\ a_1 = 0 \quad a_4 = \bar{u} \end{array} \right\} \quad (62)$$

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:

$$\left( \begin{array}{cccccc} E_1 & -E_1 & -E_1(1-\varepsilon) & 0 & 0 & 0 \\ -E_1 & E_1(2-\varepsilon) + E_2\varepsilon & 2E_1(1-\varepsilon) & -E_1(1-\varepsilon) - E_2\varepsilon & -E_2\varepsilon & 0 \\ -E_1(1-\varepsilon) & 2E_1(1-\varepsilon) & E_1\left(\frac{2-\varepsilon}{3}\right)^3 + E_2\left(\frac{\varepsilon}{3}\right)^3 & -E_1(1-\varepsilon) & -E_1\left(\frac{1-\varepsilon}{3}\right)^3 - E_2\left(\frac{\varepsilon}{3}\right)^3 & 0 \\ 0 & -E_1(1-\varepsilon) - E_2\varepsilon & -E_1(1-\varepsilon) & E_1(1-\varepsilon) + E_2(\varepsilon+1) & 2E_2\varepsilon & -E_2 \\ 0 & -E_2\varepsilon & -E_1\left(\frac{1-\varepsilon}{3}\right)^3 & -E_2\left(\frac{\varepsilon}{3}\right)^3 & 2E_2\varepsilon & -E_2\varepsilon \\ 0 & 0 & 0 & -E_2 & -E_2\varepsilon & E_2 \end{array} \right) \begin{pmatrix} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{pmatrix}$$

With the formulation of Belytschko et al. [28] it becomes:

$$\left( \begin{array}{cccccc} E_1 & -E_1 & 0 & 0 & 0 & 0 \\ -E_1 & E_1(2-\varepsilon) + E_2\varepsilon & (E_1 - E_2)(1-\varepsilon)\varepsilon & -E_1(1-\varepsilon) - E_2\varepsilon & (E_1 - E_2)(1-\varepsilon)\varepsilon & 0 \\ 0 & (E_1 - E_2)(1-\varepsilon)\varepsilon & \frac{E_1}{3}(u_1 + 2) + \frac{E_2}{3}v_1 & (E_2 - E_1)(1-\varepsilon)\varepsilon & \frac{E_1}{3}u_2 + \frac{E_2}{3}v_2 & 0 \\ 0 & -E_1(1-\varepsilon) - E_2\varepsilon & (E_2 - E_1)(1-\varepsilon)\varepsilon & E_1(1-\varepsilon) + E_2(\varepsilon+1) & (E_2 - E_1)(1-\varepsilon)\varepsilon & -E_2 \\ 0 & (E_1 - E_2)(1-\varepsilon)\varepsilon & \frac{E_1}{3}u_2 + \frac{E_2}{3}v_2 & (E_2 - E_1)(1-\varepsilon)\varepsilon & \frac{E_1}{3}(u_1 + 1) + \frac{E_2}{3}(v_1 + 1) & 0 \\ 0 & 0 & 0 & -E_2 & 0 & E_2 \end{array} \right) \begin{pmatrix} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{pmatrix}$$

$$\text{where: } \begin{cases} u_1(\varepsilon) = -\varepsilon[\varepsilon(4\varepsilon - 6) + 3] & v_1(\varepsilon) = \varepsilon[\varepsilon(4\varepsilon - 6) + 3] \\ u_2(\varepsilon) = [-2\varepsilon^3 + 3\varepsilon - 1] & v_2(\varepsilon) = \varepsilon[\varepsilon(2\varepsilon - 6) + 3] \end{cases}$$

Finally, with Moës et al. formulation [27], the discrete matrix is:

$$\left( \begin{array}{cccccc} E_1 & -E_1 & 0 & 0 & 0 & 0 \\ -E_1 & E_1(2-\varepsilon) + E_2\varepsilon & 2(E_2 - E_1)(1-\varepsilon)\varepsilon^2 & -E_1(1-\varepsilon) - E_2\varepsilon & 2(E_2 - E_1)(1-\varepsilon)^2\varepsilon & 0 \\ 0 & 2(E_2 - E_1)(1-\varepsilon)\varepsilon^2 & E_1f_1 + E_2f_2 & 2(E_1 - E_2)(1-\varepsilon)\varepsilon^2 & E_1h_1 + E_2h_2 & 0 \\ 0 & -E_1(1-\varepsilon) - E_2\varepsilon & 2(E_1 - E_2)(1-\varepsilon)\varepsilon^2 & E_1(1-\varepsilon) + E_2(\varepsilon+1) & 2(E_1 - E_2)(1-\varepsilon)^2\varepsilon & -E_2 \\ 0 & 2(E_2 - E_1)(1-\varepsilon)^2\varepsilon & E_1h_1 + E_2h_2 & 2(E_1 - E_2)(1-\varepsilon)^2\varepsilon & E_1g_1 + E_2g_2 & 0 \\ 0 & 0 & 0 & -E_2 & 0 & E_2 \end{array} \right) \begin{pmatrix} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{pmatrix}$$
* Assuming small strains, the expected solution satisfying the boundary conditions above is:

$$u(x) = \begin{cases} \frac{E_2 \bar{u}}{E_1 L + (E_2 - E_1) x_0} x, & \text{if } x \le x_0 \quad \text{so that } \alpha^- = \frac{E_2 \bar{u}}{E_1 L + (E_2 - E_1) x_0} \\ \frac{E_1 x + (E_2 - E_1) x_0}{E_1 L + (E_2 - E_1) x_0} \bar{u}, & \text{if } x \ge x_0 \quad \text{so that } \alpha^+ = \frac{E_1 \bar{u}}{E_1 L + (E_2 - E_1) x_0} \end{cases} \quad (60)$$

### 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 $\Omega = [-1, 2]$ with the following Dirichlet boundary conditions $u(-1) = 0$, $u(2) = \bar{u}$ so that $L = 3$ m. The interface is $\Gamma = \{x_0\}$.

The weak form of the problem becomes:

Find:

$$\left. \begin{array}{l} u \in \prod : \{ w \in H^1(\Omega) / w(-1) = 0, w(2) = \bar{u}, w(x_0^-) = w(x_0^+) \} \\ \forall v \in \prod_0 : \{ w \in H^1(\Omega) / w(-1) = 0, w(2) = 0, w(x_0^-) = w(x_0^+) \} \end{array} \right\} \\ \mapsto \left\{ \begin{array}{l} \int_{-1}^{2} E_i u' v' dx = 0 \\ E_1 u' (x_0^-) = E_2 u' (x_0^+) \end{array} \right. \quad (61)$$

The X-FEM discrete linear space is:

$$\prod_h : \left\{ \begin{array}{l} w/w = a_1 \Phi_1 + a_2 \Phi_2 + b_2 F \Phi_2 + a_3 \Phi_3 + b_3 F \Phi_3 + a_4 \Phi_4 \\ a_1 = 0 \quad a_4 = \bar{u} \end{array} \right\} \quad (62)$$

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:

$$\left( \begin{array}{cccccc} E_1 & -E_1 & -E_1(1-\varepsilon) & 0 & 0 & 0 \\ -E_1 & E_1(2-\varepsilon) + E_2\varepsilon & 2E_1(1-\varepsilon) & -E_1(1-\varepsilon) - E_2\varepsilon & -E_2\varepsilon & 0 \\ -E_1(1-\varepsilon) & 2E_1(1-\varepsilon) & E_1\left(\frac{2-\varepsilon}{3}\right)^3 + E_2\left(\frac{\varepsilon}{3}\right)^3 & -E_1(1-\varepsilon) & -E_1\left(\frac{1-\varepsilon}{3}\right)^3 - E_2\left(\frac{\varepsilon}{3}\right)^3 & 0 \\ 0 & -E_1(1-\varepsilon) - E_2\varepsilon & -E_1(1-\varepsilon) & E_1(1-\varepsilon) + E_2(\varepsilon+1) & 2E_2\varepsilon & -E_2 \\ 0 & -E_2\varepsilon & -E_1\left(\frac{1-\varepsilon}{3}\right)^3 & -E_2\left(\frac{\varepsilon}{3}\right)^3 & 2E_2\varepsilon & E_1\left(\frac{1-\varepsilon}{3}\right)^3 + E_2\left(\varepsilon+1\right)^3 & -E_2\varepsilon \\ 0 & 0 & 0 & -E_2 & -E_2\varepsilon & -E_2\varepsilon & E_2 \end{array} \right) \begin{pmatrix} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{pmatrix}$$

With the formulation of Belytschko et al. [28] it becomes:

$$\left( \begin{array}{cccccc} E_1 & -E_1 & 0 & 0 & 0 & 0 \\ -E_1 & E_1(2-\varepsilon) + E_2\varepsilon & (E_1 - E_2)(1-\varepsilon)\varepsilon & -E_1(1-\varepsilon) - E_2\varepsilon & (E_1 - E_2)(1-\varepsilon)\varepsilon & 0 \\ 0 & (E_1 - E_2)(1-\varepsilon)\varepsilon & \frac{E_1}{3}(u_1 + 2) + \frac{E_2}{3}v_1 & (E_2 - E_1)(1-\varepsilon)\varepsilon & \frac{E_1}{3}u_2 + \frac{E_2}{3}v_2 & 0 \\ 0 & -E_1(1-\varepsilon) - E_2\varepsilon & (E_2 - E_1)(1-\varepsilon)\varepsilon & E_1(1-\varepsilon) + E_2(\varepsilon+1) & (E_2 - E_1)(1-\varepsilon)\varepsilon & -E_2 \\ 0 & (E_1 - E_2)(1-\varepsilon)\varepsilon & \frac{E_1}{3}u_2 + \frac{E_2}{3}v_2 & (E_2 - E_1)(1-\varepsilon)\varepsilon & \frac{E_1}{3}(u_1 + 1) + \frac{E_2}{3}(v_1 + 1) & 0 \\ 0 & 0 & 0 & -E_2 & 0 & E_2 \end{array} \right) \begin{pmatrix} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{pmatrix}$$

$$\text{where: } \begin{cases} u_1(\varepsilon) = -\varepsilon[\varepsilon(4\varepsilon - 6) + 3] & v_1(\varepsilon) = \varepsilon[\varepsilon(4\varepsilon - 6) + 3] \\ u_2(\varepsilon) = [-2\varepsilon^3 + 3\varepsilon - 1] & v_2(\varepsilon) = \varepsilon[\varepsilon(2\varepsilon - 6) + 3] \end{cases}$$

Finally, with Moës et al. formulation [27], the discrete matrix is:

$$\left( \begin{array}{cccccc} E_1 & -E_1 & 0 & 0 & 0 & 0 \\ -E_1 & E_1(2-\varepsilon) + E_2\varepsilon & 2(E_2 - E_1)(1-\varepsilon)\varepsilon^2 & -E_1(1-\varepsilon) - E_2\varepsilon & 2(E_2 - E_1)(1-\varepsilon)^2\varepsilon & 0 \\ 0 & 2(E_2 - E_1)(1-\varepsilon)\varepsilon^2 & E_1f_1 + E_2f_2 & 2(E_1 - E_2)(1-\varepsilon)\varepsilon^2 & E_1h_1 + E_2h_2 & 0 \\ 0 & -E_1(1-\varepsilon) - E_2\varepsilon & 2(E_1 - E_2)(1-\varepsilon)\varepsilon^2 & E_1(1-\varepsilon) + E_2(\varepsilon+1) & 2(E_1 - E_2)(1-\varepsilon)^2\varepsilon & -E_2 \\ 0 & 2(E_2 - E_1)(1-\varepsilon)^2\varepsilon & E_1h_1 + E_2h_2 & 2(E_1 - E_2)(1-\varepsilon)^2\varepsilon & E_1g_1 + E_2g_2 & 0 \\ 0 & 0 & 0 & -E_2 & 0 & E_2 \end{array} \right) \begin{pmatrix} a_1 \\ a_2 \\ b_2 \\ a_3 \\ b_3 \\ a_4 \end{pmatrix}$$
where:

$$f_1 = \frac{4}{3}(1 - \varepsilon)\varepsilon^2[\varepsilon(4\varepsilon - 2) + 1]; \quad f_2 = \frac{16}{3}(1 - \varepsilon)^2\varepsilon^3;$$

$$h_1 = \frac{4}{3}(1 - \varepsilon)^2\varepsilon^2(-1 + 4\varepsilon)$$

$$h_2 = \frac{4}{3}(1 - \varepsilon)^2\varepsilon^2(3 - 4\varepsilon); \quad g_1 = \frac{16}{3}(1 - \varepsilon)^3\varepsilon^2; \quad g_2 = \frac{4}{3}(1 - \varepsilon)^2\varepsilon(\varepsilon(4\varepsilon - 6) + 3)$$

For all cases, we have taken $\varepsilon = 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 $\varepsilon = 0$ or $\varepsilon = 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 one-dimensional 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 \le x \le L$. Dirichlet boundary conditions are applied on this central element such that $u(0) = 0, u(L) = \bar{u}$. The linear shape functions we use in the central element are of the form $\Phi_1 = 1 - \frac{x}{L}, \Phi_2 = \frac{x}{L}$, while the quadratic ones are taken as:

<figure><img src="image_41.png" /><figcaption>**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</figcaption></figure>
where:

$$f_1 = \frac{4}{3}(1 - \varepsilon)\varepsilon^2[\varepsilon(4\varepsilon - 2) + 1]; \quad f_2 = \frac{16}{3}(1 - \varepsilon)^2\varepsilon^3;$$

$$h_1 = \frac{4}{3}(1 - \varepsilon)^2\varepsilon^2(-1 + 4\varepsilon)$$

$$h_2 = \frac{4}{3}(1 - \varepsilon)^2\varepsilon^2(3 - 4\varepsilon); \quad g_1 = \frac{16}{3}(1 - \varepsilon)^3\varepsilon^2; \quad g_2 = \frac{4}{3}(1 - \varepsilon)^2\varepsilon(\varepsilon(4\varepsilon - 6) + 3)$$

For all cases, we have taken $\varepsilon = 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 $\varepsilon = 0$ or $\varepsilon = 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 one-dimensional 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 \le x \le L$. Dirichlet boundary conditions are applied on this central element such that $u(0) = 0, u(L) = \bar{u}$. The linear shape functions we use in the central element are of the form $\Phi_1 = 1 - \frac{x}{L}, \Phi_2 = \frac{x}{L}$, while the quadratic ones are taken as:

<figure><img src="image_41.png" /><figcaption>**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</figcaption></figure>
$$ \Phi_1 = \left(1 - 2\frac{x}{L}\right) \left(1 - \frac{x}{L}\right), \quad \Phi_2 = \frac{4x}{L} \left(1 - \frac{x}{L}\right), \quad \Phi_3 = -\frac{x}{L} \left(1 - \frac{2x}{L}\right) $$

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 $\Phi(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 $v_1 = v_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 $\vec{x}$ direction, given by $\|\vec{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><caption>**Table 4** Expressions of $F(x)$ for weakly discontinuous enrichments with linear and quadratic approximations resulting from Table 3 for the one dimensional bi-material bar problem</caption><thead><tr><th></th><th>Sukumar [14]</th><th>Moës [27]</th><th>Belytschko [28]</th></tr></thead><tbody><tr><td>Linear F(x)</td><td>|x - x<sub>0</sub>|</td><td>x + x<sub>0</sub> - 2(x<sub>0</sub>/L) - |x - x<sub>0</sub>|</td><td>|x - x<sub>0</sub>| - |x<sub>J</sub> - x<sub>0</sub>|</td></tr><tr><td>Quadratic F(x)</td><td>|x - x<sub>0</sub>|</td><td>• If L/2 ≥ x<sub>0</sub><br>x<sub>0</sub> + x/L(L - 6x<sub>0</sub>)<br>+ 4x<sub>0</sub>(x/L)<sup>2</sup> - |x - x<sub>0</sub>|<br>• If L/2 ≤ x<sub>0</sub><br>x<sub>0</sub> + x/L(2x<sub>0</sub> - 3L)<br>+ 4(L - x<sub>0</sub>)(x/L)<sup>2</sup> - |x - x<sub>0</sub>|</td><td>|x - x<sub>0</sub>| - |x<sub>J</sub> - x<sub>0</sub>|</td></tr></tbody></table>
$$ \Phi_1 = \left(1 - 2\frac{x}{L}\right) \left(1 - \frac{x}{L}\right), \quad \Phi_2 = \frac{4x}{L} \left(1 - \frac{x}{L}\right), \quad \Phi_3 = -\frac{x}{L} \left(1 - \frac{2x}{L}\right) $$

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 $\Phi(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 $v_1 = v_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 $\vec{x}$ direction, given by $\|\vec{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><caption>**Table 4** Expressions of $F(x)$ for weakly discontinuous enrichments with linear and quadratic approximations resulting from Table 3 for the one dimensional bi-material bar problem</caption><thead><tr><th></th><th>Sukumar [14]</th><th>Moës [27]</th><th>Belytschko [28]</th></tr></thead><tbody><tr><td>Linear F(x)</td><td>|x - x<sub>0</sub>|</td><td>x + x<sub>0</sub> - 2(x<sub>0</sub>/L) - |x - x<sub>0</sub>|</td><td>|x - x<sub>0</sub>| - |x<sub>J</sub> - x<sub>0</sub>|</td></tr><tr><td>Quadratic F(x)</td><td>|x - x<sub>0</sub>|</td><td>• If L/2 ≥ x<sub>0</sub><br>x<sub>0</sub> + x/L(L - 6x<sub>0</sub>)<br>+ 4x<sub>0</sub>(x/L)<sup>2</sup> - |x - x<sub>0</sub>|<br>• If L/2 ≤ x<sub>0</sub><br>x<sub>0</sub> + x/L(2x<sub>0</sub> - 3L)<br>+ 4(L - x<sub>0</sub>)(x/L)<sup>2</sup> - |x - x<sub>0</sub>|</td><td>|x - x<sub>0</sub>| - |x<sub>J</sub> - x<sub>0</sub>|</td></tr></tbody></table>
<table><caption>**Table 5** Coefficients used to capture the solution with the finite element approximation for the linear patch test</caption><thead><tr><th></th><th>Sukumar [14]</th><th>Moës [27]</th><th>Belytchko [28]</th></tr></thead><tbody><tr><td>Linear φ(x) and<br>linear F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with<br>blending elements.</td><td><table><tr><td>a<sub>1</sub> = 0</td><td>a<sub>2</sub> = ū</td><td>b<sub>1</sub> = b<sub>2</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr><tr><td>Quadratic φ(x) and<br>quadratic F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with<br>blending elements.</td><td><table><tr><td>a<sub>1</sub> = 0</td><td>a<sub>2</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>ū</td></tr><tr><td>2</td></tr></tbody></table> + <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr><tr><td>Quadratic φ(x) and<br>linear F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with<br>blending elements.</td><td><table><tr><td>a<sub>1</sub> = 0</td><td>a<sub>2</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>ū</td></tr><tr><td>2</td></tr></tbody></table> + <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr><tr><td></td><td></td><td><table><tr><td>a<sub>1</sub> = 0</td><td>a<sub>2</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>ū</td></tr><tr><td>2</td></tr></tbody></table> + <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr><tr><td></td><td></td><td><table><tr><td>a<sub>1</sub> = b<sub>2</sub> = b<sub>3</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr><tr><td></td><td></td><td><table><tr><td>a<sub>1</sub> = b<sub>2</sub> = b<sub>3</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr><tr><td></td><td></td><td><table><tr><td>a<sub>1</sub> = b<sub>2</sub> = b<sub>3</sub> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>(E<sub>2</sub> - E<sub>1</sub>)ū</td></tr><tr><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td></tr></tbody></table></td></tr></tbody></table>