<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><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>• if L/2 ≥ x<sub>0</sub></td><td><table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><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><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><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>• if L/2 ≤ x<sub>0</sub></td><td><table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><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><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>)(L - x<sub>0</sub>)ū</td><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 style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><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><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table></td><td>The analytical solution<br>cannot be represented.<br>Incompatible with<br>blending elements.</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 style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><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>(E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>ū + E<sub>1</sub>Lū</td><td>2(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>ū</td><td>2</td></tr></tbody></table></td></tr><tr><td></td><td></td><td><table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>a<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>)x<sub>0</sub>ū + E<sub>1</sub>Lū</td><td>(E<sub>1</sub>L + (E<sub>2</sub> - E<sub>1</sub>)x<sub>0</sub>)</td></tr></tbody></table> = <table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>ū</td><td>(E<sub>2</sub> - E<sub>1</sub>)ū</td></tr></tbody></table></td></tr><tr><td></td><td></td><td><table style="display: inline-block; vertical-align: middle; text-align: center;"><tbody><tr><td>b<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><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><td>The analytical solution<br>cannot be represented.<br>Incompatible with<br>blending elements.</td></tr></tbody></table>
<table><caption>**Table 6** Comparison of relative error in energy norm for different interface locations between the different enrichment functions for linear elements used in the bi-material bar problem (results obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and $\vec{u} = 3 \cdot 10^{-6}$m)</caption><thead><tr><th>Interface location</th><th>Sukumar</th><th>Moës</th><th>Belytschko</th></tr></thead><tbody><tr><td>12.3</td><td>2.38688E-01</td><td>3.26389E-15</td><td>2.38688E-01</td></tr><tr><td>12.4</td><td>2.42403E-01</td><td>4.30236E-15</td><td>2.42403E-01</td></tr><tr><td>12.5</td><td>2.45626E-01</td><td>2.82096E-15</td><td>2.45626E-01</td></tr><tr><td>12.6</td><td>2.48370E-01</td><td>5.19794E-15</td><td>2.48370E-01</td></tr><tr><td>12.7</td><td>2.50650E-01</td><td>6.25056E-15</td><td>2.50650E-01</td></tr></tbody></table>

<table><caption>**Table 7** Comparison of relative error in energy norm for different interface locations between the different enrichment functions for quadratic elements used in the bi-material bar problem (results obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and $\vec{u} = 3 \cdot 10^{-6}$m)</caption><thead><tr><th>Interface location</th><th>Sukumar</th><th>Moës</th><th>Belytschko</th></tr></thead><tbody><tr><td>12.3</td><td>8.37168E-02</td><td>2.19329E-14</td><td>8.37168E-02</td></tr><tr><td>12.4</td><td>8.02193E-02</td><td>2.72702E-14</td><td>8.02193E-02</td></tr><tr><td>12.5</td><td>7.62151E-02</td><td>3.39596E-12</td><td>7.62151E-02</td></tr><tr><td>12.6</td><td>7.17862E-02</td><td>1.85735E-14</td><td>7.17862E-02</td></tr><tr><td>12.7</td><td>6.70550E-02</td><td>1.43340E-14</td><td>6.70550E-02</td></tr></tbody></table>

The exact displacements in the two materials are given by:

$$ \vec{u}_1(x) = -\frac{f}{E_1} \left( \frac{x^2}{2} - Lx \right) \vec{x} \quad \text{if } x \le x_0 $$

$$ \vec{u}_2(x) = - \left[ \frac{f}{E_2} \left( \frac{x^2}{2} - Lx \right) + f \left( \frac{x_0^2}{2} - Lx_0 \right) \left( \frac{1}{E_1} - \frac{1}{E_2} \right) \right] \vec{x} \quad \text{if } x \ge x_0 $$

No linear element with linear $F(x)$ is able to capture this quadratic solution directly. For Sukumar and Belytschko formulations with quadratic elements, the xfem blending elements which are adjacent to those completely cut by the interface give a cubic interpolation (linear $F(x)$) or order four interpolation (quadratic $F(x)$) with just one $b_i$ degree of freedom (see Fig. 8, for the left and right elements). Since the patch test is quadratic, these degrees of freedom must be equal to zero which cancels out all discontinuity. Hence

<figure><img src="image_42.png" /><figcaption>**Fig. 38** Convergence in energy norm at a rate of 0.5 [44] for Sukumar and Belytschko enrichments. These results are obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and $\vec{u} = 3 \cdot 10^{-6}$m</figcaption></figure>
<table><caption>**Table 6** Comparison of relative error in energy norm for different interface locations between the different enrichment functions for linear elements used in the bi-material bar problem (results obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and $\vec{u} = 3 \cdot 10^{-6}$m)</caption><thead><tr><th>Interface location</th><th>Sukumar</th><th>Moës</th><th>Belytschko</th></tr></thead><tbody><tr><td>12.3</td><td>2.38688E-01</td><td>3.26389E-15</td><td>2.38688E-01</td></tr><tr><td>12.4</td><td>2.42403E-01</td><td>4.30236E-15</td><td>2.42403E-01</td></tr><tr><td>12.5</td><td>2.45626E-01</td><td>2.82096E-15</td><td>2.45626E-01</td></tr><tr><td>12.6</td><td>2.48370E-01</td><td>5.19794E-15</td><td>2.48370E-01</td></tr><tr><td>12.7</td><td>2.50650E-01</td><td>6.25056E-15</td><td>2.50650E-01</td></tr></tbody></table>

<table><caption>**Table 7** Comparison of relative error in energy norm for different interface locations between the different enrichment functions for quadratic elements used in the bi-material bar problem (results obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and $\vec{u} = 3 \cdot 10^{-6}$m)</caption><thead><tr><th>Interface location</th><th>Sukumar</th><th>Moës</th><th>Belytschko</th></tr></thead><tbody><tr><td>12.3</td><td>8.37168E-02</td><td>2.19329E-14</td><td>8.37168E-02</td></tr><tr><td>12.4</td><td>8.02193E-02</td><td>2.72702E-14</td><td>8.02193E-02</td></tr><tr><td>12.5</td><td>7.62151E-02</td><td>3.39596E-12</td><td>7.62151E-02</td></tr><tr><td>12.6</td><td>7.17862E-02</td><td>1.85735E-14</td><td>7.17862E-02</td></tr><tr><td>12.7</td><td>6.70550E-02</td><td>1.43340E-14</td><td>6.70550E-02</td></tr></tbody></table>

The exact displacements in the two materials are given by:

$$ \vec{u}_1(x) = -\frac{f}{E_1} \left( \frac{x^2}{2} - Lx \right) \vec{x} \quad \text{if } x \le x_0 $$

$$ \vec{u}_2(x) = - \left[ \frac{f}{E_2} \left( \frac{x^2}{2} - Lx \right) + f \left( \frac{x_0^2}{2} - Lx_0 \right) \left( \frac{1}{E_1} - \frac{1}{E_2} \right) \right] \vec{x} \quad \text{if } x \ge x_0 $$

No linear element with linear $F(x)$ is able to capture this quadratic solution directly. For Sukumar and Belytschko formulations with quadratic elements, the xfem blending elements which are adjacent to those completely cut by the interface give a cubic interpolation (linear $F(x)$) or order four interpolation (quadratic $F(x)$) with just one $b_i$ degree of freedom (see Fig. 8, for the left and right elements). Since the patch test is quadratic, these degrees of freedom must be equal to zero which cancels out all discontinuity. Hence

<figure><img src="image_42.png" /><figcaption>**Fig. 38** Convergence in energy norm at a rate of 0.5 [44] for Sukumar and Belytschko enrichments. These results are obtained with E1 = 2.05 MPa, E2 = 20.5 MPa, L = 25 m and $\vec{u} = 3 \cdot 10^{-6}$m</figcaption></figure>
<figure><img src="image_43.png" /><figcaption>**Fig. 39** Plate with a material interface and a variable volumic load</figcaption></figure>

quadratic patch tests cannot be satisfied and are just approximated by these formulations, which will affect the energy norm as noticed by [14]. Note that Moës formulation with quadratic shape functions $\Phi(x)$ and linear $F(x)$ is able to catch the quadratic patch test solution, and that Moës formulation with quadratic shape functions $\Phi(x)$ and quadratic $F(x)$ is not able to catch it. Table 7 gives a summary of the principal results obtained for the quadratic patch test.

Figure 40 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the quadratic patch test, using linear shape functions $\Phi(x)$ and linear $F(x)$ and quadratic shape functions $\Phi(x)$ and quadratic $F(x)$. Using linear elements we obtain the optimal convergence rates of 1. Using quadratic elements we obtain 1.5. In this later case, a convergence rate of 2 would be expected since the analytical solution does not lie in the finite element space (cf. Table 8). We will not prove formally that the obtained convergence rates are in fact the expected ones, nevertheless we can give some arguments that those rates are indeed the expected ones.

In the quadratic patch test case, the first and second derivatives of the analytical solution are discontinuous. The finite element space with enrichment we use with quadratic $F(x)$ and quadratic shape functions $\Phi(x)$ does not allow to capture this solution. However it allows to capture the linear patch test solution with a discontinuous first derivative. We can deduce from this analysis that its rate of convergence on energy will be greater than 1 but lower than 2. Then the loss of 0.5 in convergence rate observed numerically for quadratic elements, with a value of 1.5 observed numerically instead of an expected 2, is coherent with the analysis for linear elements (0.5 instead of 1 when elements do not catch the discontinuity on the first derivative) performed by Moës et al. [44].

From this analysis, we recommend to use Moës formulation with quadratic elements and linear $F(x)$.

### Convergence analysis for weak discontinuities

In this paragraph we will investigate the convergence of our X-FEM formulations dedicated to weak discontinuities when the solution is not in the approximation space.
<figure><img src="image_43.png" /><figcaption>**Fig. 39** Plate with a material interface and a variable volumic load</figcaption></figure>

quadratic patch tests cannot be satisfied and are just approximated by these formulations, which will affect the energy norm as noticed by [14]. Note that Moës formulation with quadratic shape functions $\Phi(x)$ and linear $F(x)$ is able to catch the quadratic patch test solution, and that Moës formulation with quadratic shape functions $\Phi(x)$ and quadratic $F(x)$ is not able to catch it. Table 7 gives a summary of the principal results obtained for the quadratic patch test.

Figure 40 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the quadratic patch test, using linear shape functions $\Phi(x)$ and linear $F(x)$ and quadratic shape functions $\Phi(x)$ and quadratic $F(x)$. Using linear elements we obtain the optimal convergence rates of 1. Using quadratic elements we obtain 1.5. In this later case, a convergence rate of 2 would be expected since the analytical solution does not lie in the finite element space (cf. Table 8). We will not prove formally that the obtained convergence rates are in fact the expected ones, nevertheless we can give some arguments that those rates are indeed the expected ones.

In the quadratic patch test case, the first and second derivatives of the analytical solution are discontinuous. The finite element space with enrichment we use with quadratic $F(x)$ and quadratic shape functions $\Phi(x)$ does not allow to capture this solution. However it allows to capture the linear patch test solution with a discontinuous first derivative. We can deduce from this analysis that its rate of convergence on energy will be greater than 1 but lower than 2. Then the loss of 0.5 in convergence rate observed numerically for quadratic elements, with a value of 1.5 observed numerically instead of an expected 2, is coherent with the analysis for linear elements (0.5 instead of 1 when elements do not catch the discontinuity on the first derivative) performed by Moës et al. [44].

From this analysis, we recommend to use Moës formulation with quadratic elements and linear $F(x)$.

### Convergence analysis for weak discontinuities

In this paragraph we will investigate the convergence of our X-FEM formulations dedicated to weak discontinuities when the solution is not in the approximation space.
<table><caption>**Table 8 coefficients used to capture the solution with the finite element approximation for the quadratic patch test**</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>Linear φ(x) and<br>linear F(x)</td><td>Sukumar [14]</td><td>Moës [27]</td><td>Belytchko [28]</td></tr><tr><td>Quadratic φ(x) and<br>quadratic F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td><td>No solution.<br>Incompatibility for the central element.</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td></tr><tr><td>Quadratic φ(x) and<br>linear F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td><td>a<sub>1</sub> = 0<br>a<sub>2</sub> = <sup>f</sup>&frasl;<sub>4E<sub>1</sub></sub>(<sup>5</sup>&frasl;<sub>2</sub>x<sub>0</sub>L - x<sub>0</sub><sup>2</sup>) + <sup>f</sup>&frasl;<sub>4E<sub>2</sub></sub>(<sup>3</sup>&frasl;<sub>2</sub>L<sup>2</sup> - <sup>5</sup>&frasl;<sub>2</sub>x<sub>0</sub>L + x<sub>0</sub><sup>2</sup>)<br>a<sub>3</sub> = <sup>f</sup>&frasl;<sub>2E<sub>1</sub></sub>x<sub>0</sub>(2L - x<sub>0</sub>) + <sup>f</sup>&frasl;<sub>2E<sub>2</sub></sub>(x<sub>0</sub> - L)<sup>2</sup></td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td></tr></tbody></table>
<table><caption>**Table 8 coefficients used to capture the solution with the finite element approximation for the quadratic patch test**</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>Linear φ(x) and<br>linear F(x)</td><td>Sukumar [14]</td><td>Moës [27]</td><td>Belytchko [28]</td></tr><tr><td>Quadratic φ(x) and<br>quadratic F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td><td>No solution.<br>Incompatibility for the central element.</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td></tr><tr><td>Quadratic φ(x) and<br>linear F(x)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td><td>a<sub>1</sub> = 0<br>a<sub>2</sub> = <sup>f</sup>&frasl;<sub>4E<sub>1</sub></sub>(<sup>5</sup>&frasl;<sub>2</sub>x<sub>0</sub>L - x<sub>0</sub><sup>2</sup>) + <sup>f</sup>&frasl;<sub>4E<sub>2</sub></sub>(<sup>3</sup>&frasl;<sub>2</sub>L<sup>2</sup> - <sup>5</sup>&frasl;<sub>2</sub>x<sub>0</sub>L + x<sub>0</sub><sup>2</sup>)<br>a<sub>3</sub> = <sup>f</sup>&frasl;<sub>2E<sub>1</sub></sub>x<sub>0</sub>(2L - x<sub>0</sub>) + <sup>f</sup>&frasl;<sub>2E<sub>2</sub></sub>(x<sub>0</sub> - L)<sup>2</sup><br>b<sub>1</sub> = <sup>f</sup>&frasl;<sub>4</sub>(<sup>1</sup>&frasl;<sub>E<sub>1</sub></sub> - <sup>1</sup>&frasl;<sub>E<sub>2</sub></sub>)(2L - x<sub>0</sub>)<br>b<sub>2</sub> = <sup>f</sup>&frasl;<sub>8</sub>(<sup>1</sup>&frasl;<sub>E<sub>1</sub></sub> - <sup>1</sup>&frasl;<sub>E<sub>2</sub></sub>)(3L - 2x<sub>0</sub>)<br>b<sub>3</sub> = <sup>f</sup>&frasl;<sub>4</sub>(<sup>1</sup>&frasl;<sub>E<sub>1</sub></sub> - <sup>1</sup>&frasl;<sub>E<sub>2</sub></sub>)(L - x<sub>0</sub>)</td><td>The analytical solution<br>cannot be represented.<br>Incompatible with blending elements.</td></tr></tbody></table>
<figure><img src="image_44.png" /><figcaption>**Fig. 40** Convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the quadratic patch test, using linear shape functions and a linear enrichment function (linear, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic)</figcaption></figure>

### Volumic load problem with straight interface

The test we propose in the following was already investigated by Moës et al. [18]. Let us consider again the squared plate of Fig. 39 with side $L = 2$ m separated in its middle by two materials $E_1 = 1$ MPa and $E_2 = 10$ MPa and with $v_1 = v_2 = 0$. This time the variable volumic force is $||\vec{f}|| = f(\frac{x}{L})^2$ imposed on the plate along the $\vec{x}$ direction.

The exact displacements in the two materials are:

$$ \vec{u}_1(x) = \frac{f}{3E_1L^2} \left(L^3x - \frac{x^4}{4}\right) \vec{x} $$

$$ \vec{u}_2(x) = \left[ \frac{f}{3E_2L^2} \left(L^3x - \frac{x^4}{4}\right) + \frac{31fL^2}{192} \left(\frac{1}{E_1} - \frac{1}{E_2}\right) \right] \vec{x} $$

We considered an unstructured mesh with 5, 10, 20 and 40 triangular elements on each side. Figure 41 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the volumic load problem, using linear shape functions $\Phi(x)$ and linear $F(x)$, quadratic shape functions $\Phi(x)$ and linear $F(x)$ and quadratic shape functions $\Phi(x)$ and quadratic $F(x)$. Using linear elements or quadratic elements and linear enrichment function we obtain the optimal convergence rates. Using quadratic elements and quadratic enrichment function we obtain 1.5. Note that the second derivative of the analytical displacement is discontinuous in this case too.

### Circular inclusion with imposed displacement

To investigate the behavior of weak discontinuity enrichments in presence of curved interfaces, Moës et al. [18] suggested the following test case. A circular inclusion $\Omega_1$ of radius $a = 0.4$ m is placed in the center of a material $\Omega_2$. A linear displacement is imposed on the boundary $\Gamma_2$ ($r = b = 2$ m). The same material parameters are chosen: $E_1 = 1$ MPa, $v_1 = 0.25$ and $E_2 = 10$ MPa, $v_2 = 0.3$. This problem has an analytical solution provided in [18]. The analytical displacement reads:

$$ u_r(r) = \begin{cases} \left[ \left(1 - \frac{b^2}{a^2}\right) \right] \alpha + \frac{b^2}{a^2}, & \text{if } 0 \le r \le a \\ \left( r - \frac{b^2}{a^2} \right) \alpha + \frac{b^2}{a^2}, & \text{if } a \le r \le b \end{cases} $$
<figure><img src="image_44.png" /><figcaption>**Fig. 40** Convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the quadratic patch test, using linear shape functions and a linear enrichment function (linear, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic)</figcaption></figure>

### Volumic load problem with straight interface

The test we propose in the following was already investigated by Moës et al. [18]. Let us consider again the squared plate of Fig. 39 with side $L = 2$ m separated in its middle by two materials $E_1 = 1$ MPa and $E_2 = 10$ MPa and with $v_1 = v_2 = 0$. This time the variable volumic force is $\|\vec{f}\| = f(\frac{x}{L})^2$ imposed on the plate along the $\vec{x}$ direction.

The exact displacements in the two materials are:

$$ \vec{u}_1(x) = \frac{f}{3E_1L^2} \left(L^3x - \frac{x^4}{4}\right) \vec{x} $$

$$ \vec{u}_2(x) = \left[ \frac{f}{3E_2L^2} \left(L^3x - \frac{x^4}{4}\right) + \frac{31fL^2}{192} \left(\frac{1}{E_1} - \frac{1}{E_2}\right) \right] \vec{x} $$

We considered an unstructured mesh with 5, 10, 20 and 40 triangular elements on each side. Figure 41 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for the volumic load problem, using linear shape functions $\Phi(x)$ and linear $F(x)$, quadratic shape functions $\Phi(x)$ and linear $F(x)$ and quadratic shape functions $\Phi(x)$ and quadratic $F(x)$. Using linear elements or quadratic elements and linear enrichment function we obtain the optimal convergence rates. Using quadratic elements and quadratic enrichment function we obtain 1.5. Note that the second derivative of the analytical displacement is discontinuous in this case too.

### Circular inclusion with imposed displacement

To investigate the behavior of weak discontinuity enrichments in presence of curved interfaces, Moës et al. [18] suggested the following test case. A circular inclusion $\Omega_1$ of radius $a = 0.4$ m is placed in the center of a material $\Omega_2$. A linear displacement is imposed on the boundary $\Gamma_2$ ($r = b = 2$ m). The same material parameters are chosen: $E_1 = 1$ MPa, $v_1 = 0.25$ and $E_2 = 10$ MPa, $v_2 = 0.3$. This problem has an analytical solution provided in [18]. The analytical displacement reads:

$$ u_r(r) = \begin{cases} \left[ \left(1 - \frac{b^2}{a^2}\right) \alpha + \frac{b^2}{a^2}, & \text{if } 0 \le r \le a \\ \left(r - \frac{b^2}{a^2}\right) \alpha + \frac{b^2}{a^2}, & \text{if } a \le r \le b \end{cases}, $$
<figure><img src="image_45.png" /><figcaption>**Fig. 41** Convergence rates on the relative error in terms of the energy for the volumic load problem. The enrichment strategies are: linear shape functions and linear enrichment function (linear, linear), quadratic shape functions and a linear enrichment function (quadratic, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic)</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Fig. 42** Convergence rates on the relative errors with respect to mesh refinement for the circular inclusion problem. The enrichment strategies are: linear shape functions and linear enrichment function (weak, linear, linear), quadratic shape functions and linear enrichment (weak, quadratic, linear) and quadratic shape functions and quadratic enrichment function (weak, quadratic, quadratic). Two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface have also been investigated: one with quadratic shape functions and a linear approximation of the level set (strong, quadratic, linear) and one with quadratic shape functions and a quadratic approximation of the level set (strong, quadratic, quadratic)</figcaption></figure>

with:

$$ \alpha = \frac{(\lambda_1 + \mu_1 + \mu_2) b^2}{(\lambda_2 + \mu_2) a^2 + (\lambda_1 + \mu_1)(b^2 - a^2) + \mu_2 b^2} $$

The domain of computation is a square $\Omega_2$ of side $L = 2$ m. We considered an unstructured mesh with 5, 10, 20, 40, 80 and 160 triangular elements by side. Figure 42 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for circular inclusion, using linear shape functions $\Phi(x)$ and linear $F(x)$, quadratic shape functions $\Phi(x)$ and linear $F(x)$ and quadratic shape functions $\Phi(x)$ and quadratic $F(x)$. Using linear elements we obtain the optimal convergence rate. Using quadratic elements we obtain 1.5. In the case of quadratic $F(x)$, this behavior arises from the discontinuity of the second derivative of the analytical displacement. The case of linear
<figure><img src="image_45.png" /><figcaption>**Fig. 41** Convergence rates on the relative error in terms of the energy for the volumic load problem. The enrichment strategies are: linear shape functions and linear enrichment function (linear, linear), quadratic shape functions and a linear enrichment function (quadratic, linear) and quadratic shape functions and a quadratic enrichment function (quadratic, quadratic)</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Fig. 42** Convergence rates on the relative errors with respect to mesh refinement for the circular inclusion problem. The enrichment strategies are: linear shape functions and linear enrichment function (weak, linear, linear), quadratic shape functions and linear enrichment (weak, quadratic, linear) and quadratic shape functions and quadratic enrichment function (weak, quadratic, quadratic). Two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface have also been investigated: one with quadratic shape functions and a linear approximation of the level set (strong, quadratic, linear) and one with quadratic shape functions and a quadratic approximation of the level set (strong, quadratic, quadratic)</figcaption></figure>

with:

$$ \alpha = \frac{(\lambda_1 + \mu_1 + \mu_2) b^2}{(\lambda_2 + \mu_2) a^2 + (\lambda_1 + \mu_1)(b^2 - a^2) + \mu_2 b^2} $$

The domain of computation is a square $\Omega_2$ of side $L = 2$ m. We considered an unstructured mesh with 5, 10, 20, 40, 80 and 160 triangular elements by side. Figure 42 shows the convergence rates on the relative error in terms of the energy norm with respect to mesh refinement for circular inclusion, using linear shape functions $\Phi(x)$ and linear $F(x)$, quadratic shape functions $\Phi(x)$ and linear $F(x)$ and quadratic shape functions $\Phi(x)$ and quadratic $F(x)$. Using linear elements we obtain the optimal convergence rate. Using quadratic elements we obtain 1.5. In the case of quadratic $F(x)$, this behavior arises from the discontinuity of the second derivative of the analytical displacement. The case of linear
$F(x)$ is slightly different, due to its link with the geometrical approximation of the level sets and more particularly of the interface established at the end of "Enrichment functions for weak discontinuities" section. Ferté et al. [45, 46] showed for strong discontinuities that the approximation of the level set impacts the expected theoretical convergence rates, in the case of quadratic elements: 2 is expected if a quadratic approximation of the level set is used (equivalently quadratic $F(x)$), while 1.5 is expected if a linear representation of the level set is used (equivalently linear $F(x)$). To illustrate this behavior, we also considered two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface: one with quadratic shape functions and a linear approximation of the level set and one with quadratic shape functions and a quadratic approximation of the level set. The convergence rates obtained for these two models are also shown on Fig. 42.

Finally for weak discontinuities, optimal rates of convergence in energy norm are difficult to obtain for quadratic elements. If linear approximations of the level sets (equivalently linear $F(x)$) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the geometry. If quadratic approximations of the level sets (equivalently quadratic $F(x)$) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the approximation spaces. It appears that both convergence rates end up being the same, with a loss of 0.5 in energy norm with respect to optimal orders of convergence of 2. Optimal orders of convergence are recovered when using a strong discontinuity approximation space with Lagrangian multipliers to represent weak discontinuities.

## Conclusion

This paper outlines numerical issues around popular approximation spaces to model weak discontinuities or strong discontinuities and cracks, particularly with quadratic elements. We analyzed those issues on a few benchmarks. Those benchmarks revealed that popular enrichment strategies do not behave well with asymptotic configurations (when the discontinuity gets close to a node of the mesh, typically for a distance below 10% of the length of the cut edge) for strongly discontinuous enrichments and that the singular enrichment could be improved with respect to conditioning and accuracy.

Instead of working on complex and expensive strategies to solve those issues for strongly discontinuous enrichments, we preferred to apply a slight reshape of the approximation spaces, as it is in our sense, the most effective approach to deal with those issues for industrial applications. In fact, the new approximation spaces get through the proposed benchmarks, with significant improvement of the numerical behavior and accuracy of the formulations studied.

Nonetheless, more theoretical work is needed to understand why very unsmooth "bubble" functions tend to give better results than "straightforward" singular functions. At least, our results confirmed the numerical properties of "bubble" spaces denoted in the literature [25, 31, 32]. Moreover, the integration procedure needs to be improved, particularly in 3D with quadratic elements, to further the analysis on complex crack surface and crack front geometries.

At last, in case of weakly discontinuous enrichments, it was shown that most approximation spaces of the literature did not exhibit conditioning issues. However, for quadratic elements, orders of convergence in energy norm were 0.5 lower than those obtained with
$F(x)$ is slightly different, due to its link with the geometrical approximation of the level sets and more particularly of the interface established at the end of "Enrichment functions for weak discontinuities" section. Ferté et al. [45, 46] showed for strong discontinuities that the approximation of the level set impacts the expected theoretical convergence rates, in the case of quadratic elements: 2 is expected if a quadratic approximation of the level set is used (equivalently quadratic $F(x)$), while 1.5 is expected if a linear representation of the level set is used (equivalently linear $F(x)$). To illustrate this behavior, we also considered two models based on the enrichment designed for strong discontinuities, using Lagrange multipliers to enforce the continuity of the displacement through the interface: one with quadratic shape functions and a linear approximation of the level set and one with quadratic shape functions and a quadratic approximation of the level set. The convergence rates obtained for these two models are also shown on Fig. 42.

Finally for weak discontinuities, optimal rates of convergence in energy norm are difficult to obtain for quadratic elements. If linear approximations of the level sets (equivalently linear $F(x)$) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the geometry. If quadratic approximations of the level sets (equivalently quadratic $F(x)$) are used associated to a quadratic approximation of the enrichment, convergence rates are limited due to errors on the approximation spaces. It appears that both convergence rates end up being the same, with a loss of 0.5 in energy norm with respect to optimal orders of convergence of 2. Optimal orders of convergence are recovered when using a strong discontinuity approximation space with Lagrangian multipliers to represent weak discontinuities.

## Conclusion

This paper outlines numerical issues around popular approximation spaces to model weak discontinuities or strong discontinuities and cracks, particularly with quadratic elements. We analyzed those issues on a few benchmarks. Those benchmarks revealed that popular enrichment strategies do not behave well with asymptotic configurations (when the discontinuity gets close to a node of the mesh, typically for a distance below 10% of the length of the cut edge) for strongly discontinuous enrichments and that the singular enrichment could be improved with respect to conditioning and accuracy.

Instead of working on complex and expensive strategies to solve those issues for strongly discontinuous enrichments, we preferred to apply a slight reshape of the approximation spaces, as it is in our sense, the most effective approach to deal with those issues for industrial applications. In fact, the new approximation spaces get through the proposed benchmarks, with significant improvement of the numerical behavior and accuracy of the formulations studied.

Nonetheless, more theoretical work is needed to understand why very unsmooth "bubble" functions tend to give better results than "straightforward" singular functions. At least, our results confirmed the numerical properties of "bubble" spaces denoted in the literature [25, 31, 32]. Moreover, the integration procedure needs to be improved, particularly in 3D with quadratic elements, to further the analysis on complex crack surface and crack front geometries.

At last, in case of weakly discontinuous enrichments, it was shown that most approximation spaces of the literature did not exhibit conditioning issues. However, for quadratic elements, orders of convergence in energy norm were 0.5 lower than those obtained with
an equivalent strongly discontinuous enrichment, so that we recommend using the latest approach even in the case of weak discontinuities, associated to the strategy exposed in the present work to avoid conditioning issues.

**Authors' contributions**

MN conducted the numerical implementation of strong discontinuities approximations exposed in the paper. PM contributed to draft the manuscript and helped assessing the numerical properties of approximations spaces. PM helped strengthening the theoretical background and he provided numerous ideas to shape new approximations spaces. AM and SG participated in the study of weak discontinuities implementation and benchmarks. All authors read and approved the final manuscript.

**Author details**

¹IMSIA, UMR EDF/CNRS/CEA/ENSTA ParisTech 9219, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France, ²Ecole Centrale de Nantes, Laboratoire de Mécanique et Matériaux, 1 Rue de La Noe, 44321 Nantes, France, ³Rheinisch Westfalische Technische Hochschule Aachen, German Research School for Simulation Sciences GmbH, Schinkelstrasse 2, 52062 Aachen, Nordrhein-Westfalen, Germany.

**Acknowledgements**

Not applicable.

**Competing interests**

The authors declare that they have no competing interests.

**Availability of data and materials**

Not applicable.

**Consent for publication**

Not applicable.

**Ethics approval and consent to participate**

Not applicable.

**Funding**

Not applicable.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 26 October 2016 Accepted: 18 September 2017

Published online: 11 December 2017

**References**

1. Rots JG, Nauta P, Kusters GMA, Blaauwendraad J. Smeared crack approach and fracture localization in concrete. Heron. 1985;30(1):1–48.

2. Oliver J. A consistent characteristic length for smeared cracking models. Int J Num Methods Eng. 1989;28(2):461–74.

3. Ortiz M, Leroy Y, Needleman A. A finite element method for localized failure analysis. Comput Methods Appl Mech Eng. 1987;61(2):189–214.

4. Dvorkin EN, Cuitiño AM, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. Int J Num Methods Eng. 1990;30:541–64.

5. Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids. 1994;42:1397–434.

6. Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Num Methods Eng. 1999;44:1267–82.

7. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Num Methods Eng. 1999;46:135–50.

8. Duarte CA, Babuška I, Oden JT. Generalized finite element methods for three dimensional structural mechanics problems. Comput Struct. 2000;77:215–32.

9. Melenk JM, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng. 1996;139:289–314.

10. Westergaard HMW. Bearing pressures and cracks. J Appl Mech. 1939;6:A49–53.

11. Williams ML. On the stress distribution at the base of a stationary crack. ASME J Appl Mech. 1956;24:109–14.

12. Legrain G, Moës N, Verron E. Stress analysis around crack tips in finite strain problems using the extended finite element method. Int J Num Methods Eng. 2005;63(2):290–314.

13. Elguedj T, Gravouil A, Comberscure A. Appropriate extended functions for X-FEM simulation of plastic fracture mechanics. Comput Methods Appl Mech Eng. 2006;195:501–15.

14. Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. Comput Methods Appl Mech Eng. 2001;190:6183–200.

15. Gravouil A, Moës N, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets—part II: level set update. Int J Num Methods Eng. 2002;53:2569–86.
an equivalent strongly discontinuous enrichment, so that we recommend using the latest approach even in the case of weak discontinuities, associated to the strategy exposed in the present work to avoid conditioning issues.

**Authors' contributions**

MN conducted the numerical implementation of strong discontinuities approximations exposed in the paper. PM contributed to draft the manuscript and helped assessing the numerical properties of approximations spaces. PM helped strengthening the theoretical background and he provided numerous ideas to shape new approximations spaces. AM and SG participated in the study of weak discontinuities implementation and benchmarks. All authors read and approved the final manuscript.

**Author details**

¹IMSIA, UMR EDF/CNRS/CEA/ENSTA ParisTech 9219, Université Paris Saclay, 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France, ²Ecole Centrale de Nantes, Laboratoire de Mécanique et Matériaux, 1 Rue de La Noe, 44321 Nantes, France, ³Rheinisch Westfalische Technische Hochschule Aachen, German Research School for Simulation Sciences GmbH, Schinkelstrasse 2, 52062 Aachen, Nordrhein-Westfalen, Germany.

**Acknowledgements**

Not applicable.

**Competing interests**

The authors declare that they have no competing interests.

**Availability of data and materials**

Not applicable.

**Consent for publication**

Not applicable.

**Ethics approval and consent to participate**

Not applicable.

**Funding**

Not applicable.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 26 October 2016 Accepted: 18 September 2017

Published online: 11 December 2017

**References**

1. Rots JG, Nauta P, Kusters GMA, Blaauwendraad J. Smeared crack approach and fracture localization in concrete. Heron. 1985;30(1):1–48.

2. Oliver J. A consistent characteristic length for smeared cracking models. Int J Num Methods Eng. 1989;28(2):461–74.

3. Ortiz M, Leroy Y, Needleman A. A finite element method for localized failure analysis. Comput Methods Appl Mech Eng. 1987;61(2):189–214.

4. Dvorkin EN, Cuitiño AM, Gioia G. Finite elements with displacement interpolated embedded localization lines insensitive to mesh size and distortions. Int J Num Methods Eng. 1990;30:541–64.

5. Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids. 1994;42:1397–434.

6. Ortiz M, Pandolfi A. Finite-deformation irreversible cohesive elements for three-dimensional crack-propagation analysis. Int J Num Methods Eng. 1999;44:1267–82.

7. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Num Methods Eng. 1999;46:135–50.

8. Duarte CA, Babuška I, Oden JT. Generalized finite element methods for three dimensional structural mechanics problems. Comput Struct. 2000;77:215–32.

9. Melenk JM, Babuška I. The partition of unity finite element method: basic theory and applications. Comput Methods Appl Mech Eng. 1996;139:289–314.

10. Westergaard HMW. Bearing pressures and cracks. J Appl Mech. 1939;6:A49–53.

11. Williams ML. On the stress distribution at the base of a stationary crack. ASME J Appl Mech. 1956;24:109–14.

12. Legrain G, Moës N, Verron E. Stress analysis around crack tips in finite strain problems using the extended finite element method. Int J Num Methods Eng. 2005;63(2):290–314.

13. Elguedj T, Gravouil A, Comberscure A. Appropriate extended functions for X-FEM simulation of plastic fracture mechanics. Comput Methods Appl Mech Eng. 2006;195:501–15.

14. Sukumar N, Chopp DL, Moës N, Belytschko T. Modeling holes and inclusions by level sets in the extended finite element method. Comput Methods Appl Mech Eng. 2001;190:6183–200.

15. Gravouil A, Moës N, Belytschko T. Non-planar 3D crack growth by the extended finite element and level sets—part II: level set update. Int J Num Methods Eng. 2002;53:2569–86.
16. Duarte CA, Hamzeh ON, Liszka TJ, Tworzydlo WW. Generalized finite element method for the simulation of three dimensional crack propagation. Comput Methods Appl Mech Eng. 2001;190:2227-62.

17. Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for cracked domains. Int J Num Methods Eng. 2005;64:354-81.

18. Dréau K, Chevageon N, Moës N. Studied X-FEM enrichment to handle material interfaces with higher order finite element. Comput Methods Appl Mech Eng. 2010;199:1922-36.

19. Cheng KW, Fries T-P. Higher-order XFEM for curved strong and weak discontinuities. Int J Num Methods Eng. 2010;82:564-90.

20. Duarte CA, Reno LG, Simone A. A high-order generalized FEM for through-the thickness branched cracks. Int J Num Methods Eng. 2007;72(3):325-51.

21. Bechet E, Minnebo H, Moes N, Burgardt B. Improved implementation and robustness study of the X-FEM method for stress analysis around cracks. Int J Num Methods Eng. 2005;64:1033-56.

22. Strouboulis T, Babuška I, Copps K. The design and analysis of Generalized Finite Element Method. Comput Methods Appl Mech Eng. 2000;181:43-69.

23. Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Num Methods Eng. 2000;48:1741-60.

24. Siavelis M, Guiton MLE, Massin P, Moës N. Large sliding contact along branched discontinuities with X-FEM. Comput Mech. 2013;52:201-19.

25. Guptaa V, Duarte CA, Babuška I, Banerjee U. A stable and optimally convergent generalized FEM (SGFEM) for linear elastic fracture mechanics. Comput Methods Appl Mech Eng. 2013;266:23-39.

26. Fries T-P. A corrected XFEM approximation without problems in blending elements. Int J Num Methods Eng. 2008;75:503-32.

27. Moes N, Cloirec M, Cartraud P, Remacle J-F. A computational approach to handle complex microstructure geometries. Comput Methods Appl Mech Eng. 2003;192:3163-77.

28. Belytschko T, Black T, Moës N, Sukumar N, Usui S. Structured extended finite element methods of solids defined by implicit surfaces. Int J Num Methods Eng. 2003;56:609-35.

29. Nicaise S, Renard Y, Chahine E. Optimal convergence analysis for the extended finite element method. Int J Num Methods Eng. 2011;86:528-48.

30. Chevageon N, Moës N, Minnebo H. Improved crack tip enrichment functions and integration for crack modeling using the extended finite element method. J Multiscale Comput Eng. 2013;11:597-631.

31. Guptaa V, Duarte CA, Babuška I, Banerjee U. Stable GFEM (SGFEM), improved conditioning and accuracy of GFEM/XFEM for three dimensional fracture mechanics. Comput Methods Appl Mech Eng. 2015;289:355-86.

32. Babuška I, Banerjee U. Stable generalized finite element method (SGFEM). Comput Methods Appl Mech Eng. 2012;201-204:91-111.

33. Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput Methods Appl Mech Eng. 2004;193:3523-40.

34. Ling D, Yang Q, Cox B. An augmented finite element method for modeling arbitrary discontinuities in composite materials. Int J Fract. 2009;156:53-73.

35. Simone A, Duarte CA, Van der Giessen E. A generalized finite element method for polycrystals with discontinuous grain boundaries. Int J Num Methods Eng. 2006;67:1122-45.

36. Belytschko T, Black T. Elastic crack growth with minimal remeshing. Int J Num Methods Eng. 1999;45:601-20.

37. Areias P, Belytschko T. Letter to editor. Comput Methods Appl Mech Eng. 2006;195:1275-6.

38. Du Q, Wang D, Zhu L. On mesh geometry and stiffness matrix conditioning for general finite element spaces. SIAM J Num Anal. 2009;47(2):1421-44.

39. Rathod HT, Hiremath SV. Boundary integration of polynomials over an arbitrary linear tetrahedron in Euclidean three-dimensional space. Comput Methods Appl Mech Eng. 1998;153:81-106.

40. Mousavi SE, Sukumar N. Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons. Comput Mech. 2011;47:535-54.

41. Nistor I, Guiton MLE, Massin P, Moes N, Geniaut S. An X-FEM approach for large sliding contact along discontinuities. Int J Num Methods Eng. 2009;78:1407-35.

42. Destuynder P, Djaoua M, Lescure S. Quelques remarques sur la mécanique de la rupture élastique. J de Mécanique Théorique Appliquée. 1983;2(1):113-35.

43. Destuynder P, Djaoua M. Sur une interprétation mathématique de l'intégrale de Rice en théorie de la rupture fragile. Math Methods Appl Sci. 1981;3:70-87.

44. Moës N, Oden JT, Zohdi TI. Investigation of the interactions between the numerical and the modeling errors in the homogenized Dirichlet projection method. Comput Methods Appl Mech Eng. 1998;159:79-101.

45. Ferté G, Massin P, Moës N. Convergence analysis of linear or quadratic X-FEM for curved free boundaries. Comput Methods Appl Mech Eng. 2014;278:794-827.

46. Ferté G, Massin P, Moës N. Interface problems with linear or quadratic X-FEM: design of a stable multiplier space and error analysis. Int J Num Methods Eng. 2014;100(11):834-70.
16. Duarte CA, Hamzeh ON, Liszka TJ, Tworzydło WW. Generalized finite element method for the simulation of three dimensional crack propagation. Comput Methods Appl Mech Eng. 2001;190:2227-62.

17. Laborde P, Pommier J, Renard Y, Salaün M. High-order extended finite element method for cracked domains. Int J Num Methods Eng. 2005;64:354-81.

18. Dréau K, Chevageon N, Moës N. Studied X-FEM enrichment to handle material interfaces with higher order finite element. Comput Methods Appl Mech Eng. 2010;199:1922-36.

19. Cheng KW, Fries T-P. Higher-order XFEM for curved strong and weak discontinuities. Int J Num Methods Eng. 2010;82:564-90.

20. Duarte CA, Reno LG, Simone A. A high-order generalized FEM for through-the thickness branched cracks. Int J Num Methods Eng. 2007;72(3):325-51.

21. Bechet E, Minnebo H, Moes N, Burgardt B. Improved implementation and robustness study of the X-FEM method for stress analysis around cracks. Int J Num Methods Eng. 2005;64:1033-56.

22. Strouboulis T, Babuška I, Copps K. The design and analysis of Generalized Finite Element Method. Comput Methods Appl Mech Eng. 2000;181:43-69.

23. Daux C, Moës N, Dolbow J, Sukumar N, Belytschko T. Arbitrary branched and intersecting cracks with the extended finite element method. Int J Num Methods Eng. 2000;48:1741-60.

24. Siavelis M, Guiton MLE, Massin P, Moës N. Large sliding contact along branched discontinuities with X-FEM. Comput Mech. 2013;52:201-19.

25. Guptaa V, Duarte CA, Babuška I, Banerjee U. A stable and optimally convergent generalized FEM (SGFEM) for linear elastic fracture mechanics. Comput Methods Appl Mech Eng. 2013;266:23-39.

26. Fries T-P. A corrected XFEM approximation without problems in blending elements. Int J Num Methods Eng. 2008;75:503-32.

27. Moes N, Cloirec M, Cartraud P, Remacle J-F. A computational approach to handle complex microstructure geometries. Comput Methods Appl Mech Eng. 2003;192:3163-77.

28. Belytschko T, Black T, Moës N, Sukumar N, Usui S. Structured extended finite element methods of solids defined by implicit surfaces. Int J Num Methods Eng. 2003;56:609-35.

29. Nicaise S, Renard Y, Chahine E. Optimal convergence analysis for the extended finite element method. Int J Num Methods Eng. 2011;86:528-48.

30. Chevageon N, Moës N, Minnebo H. Improved crack tip enrichment functions and integration for crack modeling using the extended finite element method. J Multiscale Comput Eng. 2013;11:597-631.

31. Guptaa V, Duarte CA, Babuška I, Banerjee U. Stable GFEM (SGFEM), improved conditioning and accuracy of GFEM/XFEM for three dimensional fracture mechanics. Comput Methods Appl Mech Eng. 2015;289:355-86.

32. Babuška I, Banerjee U. Stable generalized finite element method (SGFEM). Comput Methods Appl Mech Eng. 2012;201-204:91-111.

33. Hansbo A, Hansbo P. A finite element method for the simulation of strong and weak discontinuities in solid mechanics. Comput Methods Appl Mech Eng. 2004;193:3523-40.

34. Ling D, Yang Q, Cox B. An augmented finite element method for modeling arbitrary discontinuities in composite materials. Int J Fract. 2009;156:53-73.

35. Simone A, Duarte CA, Van der Giessen E. A generalized finite element method for polycrystals with discontinuous grain boundaries. Int J Num Methods Eng. 2006;67:1122-45.

36. Belytschko T, Black T. Elastic crack growth with minimal remeshing. Int J Num Methods Eng. 1999;45:601-20.

37. Areias P, Belytschko T. Letter to editor. Comput Methods Appl Mech Eng. 2006;195:1275-6.

38. Du Q, Wang D, Zhu L. On mesh geometry and stiffness matrix conditioning for general finite element spaces. SIAM J Num Anal. 2009;47(2):1421-44.

39. Rathod HT, Hiremath SV. Boundary integration of polynomials over an arbitrary linear tetrahedron in Euclidean three-dimensional space. Comput Methods Appl Mech Eng. 1998;153:81-106.

40. Mousavi SE, Sukumar N. Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons. Comput Mech. 2011;47:535-54.

41. Nistor I, Guiton MLE, Massin P, Moes N, Geniaut S. An X-FEM approach for large sliding contact along discontinuities. Int J Num Methods Eng. 2009;78:1407-35.

42. Destuynder P, Djaoua M, Lescure S. Quelques remarques sur la mécanique de la rupture élastique. J de Mécanique Théorique Appliquée. 1983;2(1):113-35.

43. Destuynder P, Djaoua M. Sur une interprétation mathématique de l'intégrale de Rice en théorie de la rupture fragile. Math Methods Appl Sci. 1981;3:70-87.

44. Moës N, Oden JT, Zohdi TI. Investigation of the interactions between the numerical and the modeling errors in the homogenized Dirichlet projection method. Comput Methods Appl Mech Eng. 1998;159:79-101.

45. Ferté G, Massin P, Moës N. Convergence analysis of linear or quadratic X-FEM for curved free boundaries. Comput Methods Appl Mech Eng. 2014;278:794-827.

46. Ferté G, Massin P, Moës N. Interface problems with linear or quadratic X-FEM: design of a stable multiplier space and error analysis. Int J Num Methods Eng. 2014;100(11):834-70.