# Partition of unity-based discontinuous finite elements:
GFEM, PUFEM, XFEM

Angelo Simone

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

## To cite this version:

Angelo Simone. Partition of unity-based discontinuous finite elements: GFEM, PUFEM, XFEM. Revue Européenne de Génie Civil, 2007, 10.1080/17747120.2007.9692976. hal-01371130

HAL Id: hal-01371130

https://hal.science/hal-01371130

Submitted on 24 Sep 2016

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_2.png" /></figure>
# Partition of unity-based discontinuous finite elements: GFEM, PUFEM, XFEM

Angelo Simone

Faculty of Civil Engineering and Geosciences
Delft University of Technology
The Netherlands
a.simone@tudelft.nl

**RÉSUMÉ.** Dans cet article, nous présentons les principes élémentaires des éléments finis discontinus basés sur la méthode de la partition de l'unité en montrant leur relation avec la méthode des éléments finis généralisée. Les aspects liés à l'implantation dans un code ainsi que la simplicité de cette méthode sont illustrés sur un exemple unidimensionnel. La capacité de cette approche à représenter des discontinuités en déplacement indépendamment du maillage élément fini est montrée sur la propagation d'une fissure dans un milieu élastique. Nous montrons aussi certaines limitations de cette méthode quand elle est utilisée avec l'approche "dummy stiffness".

**ABSTRACT.** In this paper we review some basic notions of partition of unity-based discontinuous finite elements showing their relation to the Generalized Finite Element Method. A minimal one-dimensional example illustrates some of the issues related to the computer implementation of the method and highlights the relative simplicity of the approach. The ability of the approach in describing displacement discontinuities independently of the finite element mesh is shown in a classical crack propagation problem in an elastic medium. We also illustrate some limitations of this method when used in conjunction with the dummy stiffness approach.

**MOTS-CLÉS :** discontinuités en déplacement, GFEM, PUFEM, XFEM.

**KEYWORDS:** displacement discontinuities, GFEM, PUFEM, XFEM.
# 1. Introduction

Displacement discontinuities are incorporated in standard finite elements for the numerical treatment of geometrical discontinuities, like e.g. rock joints (Goodman et al., 1968), or for the description of problems with evolving boundaries, like e.g. cracks in quasi-brittle materials (Saouma et al., 1990) or delamination in composite materials (Alfano et al., 2001). The incorporation of a discontinuity in the displacement field is conventionally pursued by using interface elements. With interface elements, a discontinuous displacement field is described by the relative displacement of a double set of nodes at inter-element boundaries. This requires special mesh generators and, in the case of evolving boundaries, remeshing procedures (Bittencourt et al., 1996). Alternatively, a discontinuous displacement field can be described by means of embedded discontinuity elements in which the displacement discontinuity, described by using the Heaviside function, is incorporated into the finite-element formulation as an incompatible strain mode and can pass through solids elements (Simo et al., 1993; Jirásek, 2000).

A different approach employs the partition of unity property of finite-element shape functions (the sum of the shape functions must equal unity at each spatial point). Within this approach, the standard approximation basis is enriched locally with special functions. This enrichment results in extra degrees of freedom for the nodes in the domain subjected to the enrichment, without modification of the mesh topology (Babuška et al., 1994; Melenk et al., 1996; Babuška et al., 1997; Duarte et al., 1996b; Duarte et al., 1996a; Oden et al., 1998; Oden et al., 1997). The effect of a displacement discontinuity is described by enriching the standard finite-element polynomial basis with a step function (Moës et al., 1999). The outcome of this approach is a class of elements, known in the literature under the names “generalized finite element method” (GFEM) (Simone et al., 2006), “partition of unity finite element method” (PUFEM) (Gasser et al., 2005) and “extended finite element method” (XFEM) (Moës et al., 1999), which is kinematically equivalent to the class of conventional interface elements, the key difference being the possibility of arbitrarily locating the discontinuity within the domain of an element. When PU-based discontinuous elements are used, the interface response is described by an additional set of global degrees of freedom and by a constitutive law at the discontinuity.

# 2. Beyond the classical finite element method: Exploiting partitions of unity

In this section, we recall some basic notions on the partition of unity paradigm and its relation to finite element methods. A partition of unity is defined as a collection of global functions $N_I(x)$ whose value sums up to unity at each point $x$ in the solution domain $\bar{V}$ of a boundary value problem:

$$ \sum_{I=1}^{n} N_{I}(x) = 1 \quad \forall x \in \bar{V}. \qquad [1] $$
where $n$ is the number of nodes of the discretization. In the finite element method (FEM) [1] holds by construction for any Lagrangian FEM global basis (i.e. shape function). Using linear partitions of unity $N_I$ as basis functions in a FEM framework (Taylor et al., 1998), a scalar field $q(x)$ can be approximated by

$$ q_h(x) = \underbrace{\sum_{I=1}^{n} N_I(x) \hat{q}_I}_{\text{regular interpolation}} + \underbrace{\sum_{I=1}^{n} N_I(x) \sum_{j=1}^{m} E_j(x) \tilde{q}_{Ij}}_{\text{enrichment}}, \quad [2] $$

where the subscript “h” denotes the approximate field, $\hat{q}_I$ (termed “regular” or “standard” degrees of freedom (dofs)) are nodal values related to the basis $N_I$, and $\tilde{q}_{Ij}$ (termed “enrichment” or “extra” degrees of freedom) are the global nodal parameters connected to the basis $E_j$ pertinent to the enrichment ($m$ is the number of terms in $E_j$). The terms “regular” and “enrichment” make reference to the fact that the “regular” interpolation field is considered as the background field upon which the “enrichment” is added. To avoid linear dependency, the order of any polynomial terms in the enrichment basis must be greater than the order of the partitions of unity $N_I$. This is, however, not a sufficient condition; algorithms to deal with the linear dependent case are described in Duarte et al. (2000).

Shape functions associated with [2] are termed GFEM shape functions. These shape functions are built from the product of a partition of unity $N_I$ and enrichment functions $E_j$:

$$ N_I \times \{1, E_j\}. \quad [3] $$

In the GFEM, the partition of unity is in general provided by linear Lagrangian FEM shape functions. The partition of unity property of the functions $N_I$ implies that

$$ \sum_{I=1}^{n} N_I(x) E_j(x) = E_j(x). \quad [4] $$

Consequently, any enrichment function $E_j$ can be represented exactly through linear combinations of GFEM shape functions, and, if the enrichment functions $E_j$ can approximate the solution of a boundary value problem, the corresponding GFEM shape functions also will.

In finite element notation, the approximation for a vector field $q(x)$ of order $l$ over an $n$-node enriched element, with all the nodes enriched, can be written as

$$ q_h(x) = N(x)\hat{q} + N(x)N_E(x)\tilde{q}, \quad [5] $$

where $N$ is a $l \times (l \times n)$ matrix containing the standard finite element shape functions, $N_E$ is a $(l \times n) \times (l \times m \times n)$ matrix containing the extra basis terms for the enrichment, $\hat{q}$ is a $(l \times n) \times 1$ vector containing standard degrees of freedom and $\tilde{q}$ is a $(l \times m \times n) \times 1$ vector containing the extra degrees of freedom. The number of extra degrees of freedom per node is equal to the number of terms in the basis for the enrichment
multiplied by the number of nodal unknowns. In standard finite elements, the matrix $N_E$ is empty.

As noted in Oden et al. (1998) and Taylor et al. (1998), the approach of Equation [5] allows the enrichment to be performed from node to node in a mesh by activating the extra degrees of freedom $\tilde{q}$ when needed (a hierarchical finite element formulation based on the partition of unity method).

## 2.1. Enrichment functions

The essence of the generalized FEM is condensed in [2]. This special approximation allows to paste a-priori knowledge about the behavior of the solution of a problem directly into the shape functions. Using this feature, polynomial functions, discontinuous functions and even singular functions have been exploited as enrichments. These customized enrichment functions have been used to model local features like cracks, voids, inclusions, microstructures etc.

Figure 1 shows the construction of a generalized FEM shape function as the product of a partition of unity and an enrichment function for two different enrichment functions. Next, we concentrate on the discontinuous enrichment.

## 2.2. Kinematics for a body crossed by a discontinuity

To incorporate a discontinuity in the displacement field it is necessary to characterize the kinematics of a body crossed by a discontinuity. A body $\bar{V}$ bounded by $S$ and crossed by a single discontinuity surface $S_d$ is considered (see Figure 2). The internal discontinuity surface $S_d$ divides the body into two sub-domains, $V^+$ and $V^-$ ($V = V^+ \cup V^-$). Displacements $\bar{u}$ are prescribed on $S_u$, while tractions $\bar{T}$ are prescribed on $S_t$. The displacement field can be decomposed as

$$u(x, t) = \hat{u}(x, t) + \mathcal{H}(x)\tilde{u}(x, t) \quad \text{if } x \notin S_d, \qquad [6]$$

where $\mathcal{H}(x)$ is the Heaviside function centered at the discontinuity surface $S_d$ ($\mathcal{H} = 1$ if $x \in V^+$, $\mathcal{H} = 0$ if $x \in V^-$) and $\hat{u}$ and $\tilde{u}$ are continuous functions on $\bar{V}$. The discontinuity is introduced by the Heaviside function $\mathcal{H}$ at the discontinuity surface $S_d$ and the value of the displacement jump $[u]$ across the discontinuity surface is given by $\tilde{u}$ at $S_d$. In the geometrically linear case, the strain field in $V$ is computed as the symmetric part of the gradient of the displacement field and reads

$$\varepsilon(x, t) = \nabla^s \hat{u}(x, t) + \mathcal{H}(x) \nabla^s \tilde{u}(x, t) \quad \text{if } x \notin S_d, \qquad [7]$$

where $(\cdot)^s$ refers to the symmetric part of $(\cdot)$. The strain field is defined everywhere in the body except at the discontinuity surface.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** A GFEM shape function $\phi_{Ii}$ as the product of a partition of unity $N_I$ and an enrichment function $E_i$ (Duarte et al., 2006) (the enrichment takes the form of a polynomial function in (a) and of a discontinuous function in (b))</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Body $\bar{V}$ crossed by a discontinuity $S_d$</figcaption></figure>
### 3. Describing a displacement discontinuity

In describing a displacement discontinuity, and considering [5] now rewritten for convenience as

$$u_h = Na + NN_Eb, \quad [8]$$

if the standard displacement field is interpolated by the regular interpolation $Na$, the displacement jump can be described through the enrichment $NN_Eb = \mathcal{H}NHb$ (Moës et al., 1999), where $\mathcal{H}$ is the scalar valued Heaviside function ($m = 1$) and $H$ is a diagonal matrix identifying through a 1/0 switch which degrees of freedom to enrich ($H$ is the identity matrix if all the degrees of freedom of the element are enriched). This results in the $(l \times n) \times 1$ vector of extra degrees of freedom $b$ (same dimensions as the vector of standard degrees of freedom $a$) and in the $(l \times n) \times (l \times n)$ matrix $\mathcal{H}H$. The enrichment concerns only nodes whose support is crossed by a discontinuity. For nodes whose support is not crossed by a discontinuity, the basis responsible for the enrichment is empty since the Heaviside function is a constant function over their supports and can be neglected. Since element stiffness matrices for elements with extra degrees of freedom are assembled only for active extra degrees of freedom, the standard matrix $N$ can be used in place of $NH$. Hence, the discretized format of [6] reads, for points located in an element with enriched nodes,

$$u_h = Na + \mathcal{H}Nb, \quad [9]$$

where the global nodal degrees of freedom $a$ and $b$ represent, in the arrangement of [9], the total displacement field. The displacement jump across the discontinuity $S_d$ is given by

$$[u_h] = Nb \Big|_{S_d}. \quad [10]$$

For points located in an element without enriched nodes, the Heaviside function is a constant function over their supports and therefore it is not considered. Consequently, since there is no enhancement, the standard finite-element interpolation $u_h = Na$ is retrieved.

It is worth noting that in this case the GFEM shape functions [3] are built from the product of a partition of unity $N_I$ and the step function $H$ as shown in Figure 1b:

$$N_I \times \{1, H\}. \quad [11]$$

The case of multiple discontinuities is easily dealt with by adding the corresponding enrichment functions (a detailed description of the procedure in case of intersecting discontinuities is reported in Simone et al. (2006) and Simone et al. (2006)).
<figure><img src="image_5.png" /><figcaption>Figure 3. Nodes affected by enrichment in a crack problem</figcaption></figure>

## 4. Near-tip enrichment

So far, we have considered the enrichment across a crack surface. However, with reference to the crack geometry depicted in Figure 3 care should be taken to properly describe the singularities at the crack tip in case of a linear elastic sample. As proposed by Belytschko et al. (1999), a set of near-tip functions that spans the exact asymptotic displacement field near a crack-tip in a two dimensional linear elastic body subjected to general mixed-mode loading can be used as enrichment around the crack tip. One such a set is given by the functions

$$E(r, \theta) = \left[ \sqrt{r} \sin \frac{\theta}{2}, \quad \sqrt{r} \cos \frac{\theta}{2}, \quad \sqrt{r} \sin \frac{\theta}{2} \sin \theta, \quad \sqrt{r} \cos \frac{\theta}{2} \sin \theta \right], \quad [12]$$

where $(r, \theta)$ are the local polar coordinates at the crack as shown in Figure 3.

With these functions as enrichment at the crack tip, the complete displacement field is expressed as (Moës et al., 1999)

$$\mathbf{u}^h = \underbrace{\sum_{i \in I} \mathbf{a}_i N_i}_{\text{continuous}} + \underbrace{\sum_{j \in J} \mathbf{b}_j N_j H}_{\text{discontinuous}} + \underbrace{\sum_{k \in K} N_k \left( \sum_{l=1}^{4} c_{kl} E_l \right)}_{\text{tip}}, \quad [13]$$

where $I$ is the set of standard nodes, $J$ is the set of circled nodes and $K$ is the set of squared nodes reported in Figure 3. Note that the use of the functions reported in [12] is limited to two-dimensional linear elastic fracture mechanics and should not be used in other situations. Also, the inclusion of the near-tip field in [12] is an important advantage over other approaches in which only the discontinuity in the displacement field is considered (cf. embedded discontinuities models (Jirásek, 2000)).

## Hp discretization

Near-tip enrichment functions are not always available. Indeed, the stress state in the neighborhood of a crack front is not well known in three-dimensions, and analyti-
cal expansions are available only for particular crack geometries (Broberg, 1999; Mit-
telstedt et al., 2005). Complex crack front geometries, curved crack surfaces or the in-
tersection of the crack surface with the boundary, create complex stress distributions
that are, in general, not amenable to closed form expansions. Similarly, non-linear
fracture models generate stress states that are not amenable to analytical expansions
in the neighborhood of a process zone. Two additional difficulties with near crack tip
enrichment functions are their numerical integration and the need to enrich several
layers of elements around the crack front in order to achieve optimal convergence
rates (Laborde et al., 2005; Béchet et al., 2005; Xiao et al., 2006). Integration er-
rors may lead to wrong conclusions in a convergence analysis since they may cancel
discretization errors when computing, for example, the strain energy.

A strategy to cope with the lack of general near-tip enrichments, consists in re-
solving the non-smooth displacement field near a crack tip/front by local mesh refi-
nements combined with high-order approximations (Szabó, 1986) (*hp* discretization).
This approach avoids the difficulties and limitations related to near crack front enrich-
ment functions while being able to deliver high convergence rates as demonstrated
in Duarte *et al.* (2007). The construction of *hp* discretizations in the neighborhood of
a crack front is not difficult in the generalized FEM since crack surfaces are allowed
to cut elements. The mesh refinement can, therefore, be performed *as if there were no*
*cracks in the domain*.

**5. High-order approximations**

At variance with the established enrichment of the displacement field with dis-
continuous functions and the use of linear approximations proposed by Moës et al.
(1999) in the so-called XFEM (see [13]), approximation orders higher than one have
been exploited in various forms by Wells et al. (2001), Stazi et al. (2003), Chessa et
al. (2003), Mariani et al. (2003), and Laborde et al. (2005).

The performance of the method can be improved considerably by consistently
using high-order hierarchical GFEM approximations based on a linear finite element
partition of unity. In this situation, all nodal degrees of freedom are defined at element
vertex nodes which facilitates mesh generation, especially in three-dimensional com-
putations. This is in contrast with the use of high-order non-hierarchical Lagrangian
partition of unity functions with piece-wise constant enrichment functions (Wells et
al., 2001; Stazi et al., 2003; Chessa et al., 2003; Laborde et al., 2005). This approach
leads to nodes on element edges, faces and interior. Hierarchical approximations are
also ideal for the construction of non-uniform polynomial approximations which can
be very effective for some classes of problems (Duarte et al., 2002). Other advantages
of hierarchical GFEM, such as the non-zero structure of the resulting stiffness and
mass matrices, are discussed in Taylor et al. (1998), Oden et al. (1998), and Duarte et
al. (2000). The issue of linear dependence of hierarchical GFEM enrichments is also
discussed in detail in Duarte et al. (2000).
Below, we list a few examples of GFEM approximations for single cracks in two-dimensions. Hereafter, $N_I, I = 1, ..., n$ are standard Lagrangian shape functions, $x_I = (x_I, y_I)$ are the coordinates of the node $I$, $h_I$ is the diameter of the largest finite element sharing the node $I$, $n$ is the number of nodes in the mesh, $p$ and $q$ are used to denote the polynomial degree of $\hat{u}_h$ and $\tilde{u}_h$, respectively.

1) $p = 1, q = 1$: the shape functions at a node $x_I$ are given by

$$N_I \times \{1, \mathcal{H}\} \quad [14]$$

which are the same generalized FE shape functions in [11] i.e. we recover the simplest possible case, that of a piecewise linear approximation of $u$. These shape function will be used in the remainder of the paper.

2) $p = 2, q = 0$: this case corresponds to the generalized FE shape functions introduced in Oden et al. (1998) and Duarte et al. (2000). The shape functions at a node $x_I$ are given by

$$N_I \times \left\{ 1, \frac{(x - x_I)}{h_I}, \frac{(y - y_I)}{h_I} \right\}. \quad [15]$$

These shape functions, of course, can not represent the discontinuous components of the solution $\tilde{u}_h$.

3) $p = 2, q = 2$: the shape functions at a node $x_I$ are given by

$$N_I \times \left\{ 1, \frac{(x - x_I)}{h_I}, \frac{(y - y_I)}{h_I}, \mathcal{H}, \mathcal{H} \frac{(x - x_I)}{h_I}, \mathcal{H} \frac{(y - y_I)}{h_I} \right\}. \quad [16]$$

4) $p = 2, q = 1$: the shape functions at a node $x_I$ are given by

$$N_I \times \left\{ 1, \frac{(x - x_I)}{h_I}, \frac{(y - y_I)}{h_I}, \mathcal{H} \right\}. \quad [17]$$

These shape functions are equivalent to those proposed by Stazi et al. (2003).

A study of the performance of the above shape functions can be found in Duarte et al. (2007).

## 6. Governing equations

Having characterized the displacement field for a body crossed by a discontinuity, we now briefly recall the governing equations for such a body. The equilibrium equations and boundary conditions for the body $\bar{V}$ without body forces depicted in Figure 2 can be summarized by

$$\nabla \cdot \boldsymbol{\sigma} = \mathbf{0} \quad \text{in } V, \quad [18a]$$

$$\boldsymbol{\sigma} \mathbf{n} = \bar{\mathbf{t}} \quad \text{on } S_t, \quad [18b]$$

$$\boldsymbol{\sigma} \mathbf{n}_d = \mathbf{t} \quad \text{on } S_d, \quad [18c]$$
where $\sigma$ is the Cauchy stress tensor, $n$ is the outward unit normal to the body and $n_d$ is the inward unit normal to $V^+$ on $S_d$. Equation [18c] represents traction continuity at the discontinuity surface $S_d$. The strong form is completed by the boundary conditions

$$ u = \bar{u} \quad \text{on } S_u, \qquad [19a] $$

$$ \tilde{u} = 0 \quad \text{on } S_u, \qquad [19b] $$

where $\bar{u}$ is a prescribed displacement. Equation [19a] is a standard essential boundary condition while [19b] has been imposed to simplify the finite-element implementation (Wells et al., 2001). The constitutive relationships for bulk and discontinuity will be specified later.

Following standard procedures, the governing equations [18a] to [18c] can be cast in a weak form. To this end, the equilibrium equations [18a] are multiplied by the weight function $w_u$, which is decomposed into $\hat{w}_u$ and $\tilde{w}_u$ consistent with the displacement decomposition in [6], and integrated over the domain $V$ to obtain a weak equilibrium statement which can be expanded using integration by parts, Gauss' theorem and [18b]. More details on the derivation of the discrete weak governing equations can be found in Wells et al. (2001) and Simone et al. (2006).

## Discretized governing equations

Following a Bubnov-Galerkin approach, the linearized form of the discretized weak governing equation is obtained by substituting in the discrete weak governing equations the stress rate in the bulk

$$ \dot{\sigma} = D_e \dot{\varepsilon} = D_e (B\dot{a} + \mathcal{H}B\dot{b}), \qquad [20] $$

with $D_e$ the tangent matrix for the bulk material, and the traction rate at a discontinuity

$$ \dot{t} = T[\dot{u}] = T(N\dot{b}) \Big|_{S_d}, \qquad [21] $$

where $T$ relates traction rate $\dot{t}$ and displacement jump rate $[\dot{u}]$.

After standard manipulations, the linearized weak form of the governing equations reads (cf e.g. Wells et al. (2001) and Simone et al. (2006))

$$ \begin{bmatrix} K_{aa} & K_{ab} \\ K_{ba} & K_{bb} \end{bmatrix} \begin{bmatrix} \delta a \\ \delta b \end{bmatrix} = \begin{bmatrix} f_{\text{ext},a} \\ f_{\text{ext},b} \end{bmatrix} - \begin{bmatrix} f_{\text{int},a} \\ f_{\text{int},b} \end{bmatrix}, \qquad [22] $$
where

$$K_{aa} = \int_V B^T D_e B \, dV, \qquad [23a]$$

$$K_{ab} = \int_{V^+} B^T D_e B \, dV, \qquad [23b]$$

$$K_{ba} = K_{ab}^{T} = \int_{V^+} B^T D_e B \, dV, \qquad [23c]$$

$$K_{bb} = \int_{V^+} B^T D_e B \, dV + \int_{S_d} N^T T N \, dS, \qquad [23d]$$

and

$$f_{\text{int},a} = \int_V B^T \boldsymbol{\sigma} \, dV, \qquad f_{\text{ext},a} = \int_{S_t} N^T \bar{\boldsymbol{t}} \, dS, \qquad [24a]$$

$$f_{\text{int},b} = \int_{V^+} B^T \boldsymbol{\sigma} \, dV + \int_{S_d} N^T \boldsymbol{t} \, dS, \qquad f_{\text{ext},b} = \int_{S_t^+} N^T \bar{\boldsymbol{t}} \, dS. \qquad [24b]$$

In the above relations, $B$ is a matrix containing derivatives of shape functions, $\boldsymbol{\sigma}$ is an array containing the components of the stress tensor in engineering notation and $\boldsymbol{t}$ is the array of the traction forces across the discontinuity.

Note that symmetry of the stiffness matrix is assured if the material tangent matrices $D_e$ and $T$ are symmetric. As a final remark, the arrays $N$ multiplying $a$ and $b$ in [9] are normally not the same since only part of the extra degrees of freedom in the array $b$ might be activated. Consequently, the matrices $B$ multiplying $a$ and $b$ in [20] are also normally not the same. However, since the system of equations in [22] is assembled only for the active degrees of freedom, it is still possible to use the standard $N$ and $B$ matrices. Care should of course be taken in considering only active degree of freedom in the strain computation. More details on the finite-element implementation can be found in Moës et al. (1999), Wells et al. (2001), Simone et al. (2006), and Sukumar et al. (2003).

## 7. Finite element implementation

PU-based discontinuous elements require some important modifications to existing finite element codes. Some of the issues are illustrated in the example reported in Section 8. Here we make reference to some issues relevant to static and propagating discontinuities.

*Numerical integration:* The integration of element matrices on part of the element domain requires the definition of special integration schemes. A common strategy consists in considering quadrature subcells aligned with the discontinuity line (Moës et al., 1999). When an element is intersected by a discontinuity, the two resulting do-
<figure><img src="image_6.png" /><figcaption>**Figure 4.** *Composite integration scheme for a quadrilateral element crossed by a discontinuity (the crossed circles are integration points on the discontinuity and the black dots are integration points for the continuum)*</figcaption></figure>

mains are triangulated and each triangular sub-domain is mapped to a parent unit triangle over which a three-point symmetric quadrature rule, with interior points within the triangular sub-domain, is considered (see Figure 4). The discontinuity contribution to the stiffness matrix and the internal force vector are integrated on the discontinuity through a two-point Gauss integration scheme. For elements not crossed by a discontinuity, a standard Gauss integration scheme is used.

It is worth noting that the resulting composite integration scheme for the continuum reported in Figure 4, which consists of 24 integration points in eight triangular sub-domains, integrates correctly a quadratic function over a parent quadrilateral element. In the numerical integration of the element matrices, this composite integration scheme results, for the terms containing the strain energy density, in a full integration rule for bilinear quadrilateral elements and in a reduced integration rule for a quadratic quadrilateral elements. The actual construction of the composite integration rule is based on an exact algebraic inverse iso-parametric mapping for bilinear quadrilateral elements. This technique can be used for quadratic quadrilateral elements provided that the element has straight sides and evenly spaced side nodes.

An alternative to the above procedure is to use a high order integration scheme in the elements so that there is at least one integration scheme in each domain. More recently, Ventura (2006) has proposed an approach in which discontinuous shape functions are replaced by equivalent polynomials thus eliminating the need for integration subcells.

*Propagating discontinuities:* Most commonly used criteria for discontinuity extension are based on stress intensity factors (Moës *et al.*, 2002), principal stresses (Jirásek *et al.*, 2001; Wells *et al.*, 2001) and loss of hyperbolicity (Belytschko *et al.*, 2003). Criteria based on damage or plastic strain accumulation can be used (Simone *et al.*, 2003; Simone *et al.*, 2004; Comi *et al.*, 2007) when the bulk material is inelastic.

As for the crack tip location, the two possible alternatives are to locate a crack tip within an element (Moës *et al.*, 1999) or at its boundary (Wells *et al.*, 2001). The latter
approach has the advantage of greatly simplifying the computational implementation
by setting the displacement jump at the discontinuity to zero. This is achieved by
considering only standard **a** dofs for the nodes on an element boundary touched by
a discontinuity (see Figure 5). When a discontinuity is extended in a neighboring
element, the nodes behind the discontinuity tip are enhanced.

The direction of discontinuity extension is computed using the principal directions
of a nonlocal stress tensor which is calculated as a weighted average of stresses using
a Gaussian weight function (Jirásek *et al.*, 2001; Wells *et al.*, 2001) with an the inter-
action radius equal to three times the average element size ahead the discontinuity tip.
The discontinuity propagates in the direction normal to the direction of the maximum
nonlocal principal stress. When the discontinuity is close to a boundary, the disconti-
nuity extension direction is aligned with the previous discontinuity segment to avoid
an incorrect direction determination due to the bias produced by an unsymmetric do-
main in the weighting procedure.

Finally, to preserve the robustness of the Newton-Raphson solution procedure, a discontinuity is introduced as a straight segment at the end of a load increment in the element ahead of a discontinuity tip if the initiation criterion is fulfilled. This procedure is repeated in the elements ahead of this element with the extended discontinuity until the initiation criterion is no longer satisfied.

Many important details have been glossed over for the sake of brevity in this short
summary; the interested reader is referred to Moës *et al.* (1999), Dolbow *et al.* (2000),
Sukumar *et al.* (2003), and Simone *et al.* (2006) for further details.

**8. A minimal one-dimensional example**

Some of the issues related to the implementation of PU-based discontinuous ele-
ments are illustrated by means of a simple one-dimensional example of a bar in ten-

<figure><img src="image_7.png" /><figcaption>**Figure 5.** *Extension of a discontinuity within the body V̄ (the black circles indicate nodes with standard degrees of freedom α while the white circles indicate nodes with extra degrees of freedom β for the nodes whose support is crossed by a discontinuity)*</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6. Geometry and boundary conditions for the tension test</figcaption></figure>

sion as depicted in Figure 6. The solution will be given for standard and for PU-based discontinuous elements. The bar is modelled by means of the two discretizations depicted in Figure 7. The spring is first represented by a conventional interface element (Figure 7a) while the discontinuous interpolation of [9] is exploited for the discretization depicted in Figure 7b. The domain $V^+$ is defined by $0 < x < L$ in Figure 6. In the following, $E$ is the Young's modulus, $A$ is the cross-sectional area and $d$ is the spring stiffness. The element types needed for the finite-element solution are depicted in Figure 8.

## 8.1. Stiffness matrix computation.

The stiffness matrix of the one-dimensional truss element of length $L$ depicted in Figure 9a is

$$ K_{\text{truss}} = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. \qquad [25] $$

The same stiffness matrix describes a one-dimensional conventional interface element which can be conceived as a translational spring element with stiffness $d = EA/L$ in [25].

The sub-matrices in [23] are expanded for the truss element with a discontinuity in the middle section depicted in Figure 9b. Note that the positive part $V^+$ of the domain goes from $x = 0$ to $x = L/2$ ($H = 1$ for $x < L/2$; see Figure 9b) and the presence of extra degrees of freedom $b$ for both nodes. Sub-matrix $K_{aa}$ is the same as $K_{\text{truss}}$. The remaining sub-matrices are expanded as

$$ \begin{align*} K_{ab} = K_{ba} &= \int_{V^+} B^T D_e B \ dV = \int_{0}^{L/2} B^T D_e B \ dx \\ &= \frac{EA}{2L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}. \quad [26] \end{align*} $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Discretization for the tension test with (a) conventional interface elements and (b) PU-based discontinuous elements (the black circles indicates the standard degrees of freedom **a** while the white circles indicates the extra degrees of freedom **b**; the discontinuity is indicated by the vertical dotted line).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Element assembly for the discretization with PU-based discontinuous elements (note the effect of the choice of the domain $V^+$ on the activation of the extra degrees of freedom **b**)</figcaption></figure>

In $K_{bb}$, the volume integral equals $K_{ab}$ while the surface integral is evaluated on $x = L/2$ and yields

$$ \int_{S_d} N^T T N \, dS = (N^T T N) \Big|_{x=L/2} = \frac{d}{4} \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}. \qquad [27] $$

Assembly of the sub-matrices into the element stiffness matrix yields

$$ K = \begin{bmatrix} K_{aa} & K_{ab} \\ K_{ba} & K_{bb} \end{bmatrix} = \begin{bmatrix} \frac{EA}{L} & -\frac{EA}{L} & \frac{EA}{2L} & -\frac{EA}{2L} \\ -\frac{EA}{L} & \frac{EA}{L} & -\frac{EA}{2L} & \frac{EA}{2L} \\ \frac{EA}{2L} & -\frac{EA}{2L} & \frac{EA}{2L} + \frac{d}{4} & -\frac{EA}{2L} + \frac{d}{4} \\ -\frac{EA}{2L} & \frac{EA}{2L} & -\frac{EA}{2L} + \frac{d}{4} & \frac{EA}{2L} + \frac{d}{4} \end{bmatrix}. \quad [28] $$
<figure><img src="image_11.png" /><figcaption>**Figure 9.** One-dimensional (a) truss element and (b) PU-based discontinuous truss element</figcaption></figure>

If the discontinuity does not cross the element but extra degrees of freedom b are active at one of the element nodes (see e.g. Figure 8a), the element stiffness matrix can be derived following the procedure described above and reads as

$$K = \frac{2EA}{L} \begin{bmatrix} 1 & -1 & -1 \\ -1 & 1 & 1 \\ -1 & 1 & 1 \end{bmatrix}, \qquad [29]$$

where the degrees of freedom have been ordered in the sequence $[ a_1 \ a_2 \ b_2 ]$ (note that the truss length is $L/2$).

## 8.2. Assembly and solution.

For the bar with the spring modelled by a conventional interface element (see Figure 7a), assembly of local stiffness matrices into the global stiffness matrix for the active degrees of freedom results in the system of equations

$$ \begin{bmatrix} \frac{EA}{L} + d & -d & 0 \\ -d & \frac{EA}{L} + d & -\frac{EA}{L} \\ 0 & -\frac{EA}{L} & \frac{EA}{L} \end{bmatrix} \begin{bmatrix} a_2 \\ a_3 \\ a_4 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ P \end{bmatrix}, \qquad [30] $$

which yields

$$ \begin{bmatrix} a_2 & a_3 & a_4 \end{bmatrix} = \begin{bmatrix} \frac{PL}{EA} & \frac{PL}{EA} + \frac{P}{d} & \frac{2PL}{EA} + \frac{P}{d} \end{bmatrix}. \qquad [31] $$
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Displacement field for the one-dimensional tension test: (a) total displacement field given by the sum of (b) the $\hat{u}$ field and (c) the $\tilde{u}$ field (see [32] for the numerical values)</figcaption></figure>

For the discretization with the PU-based discontinuous element (see Figure 7b), the global system of equations reads

$$ \begin{bmatrix} \frac{3EA}{L} & -\frac{EA}{L} & \frac{5}{2}\frac{EA}{L} & -\frac{1}{2}\frac{EA}{L} & 0 \\ -\frac{EA}{L} & \frac{3EA}{L} & -\frac{1}{2}\frac{EA}{L} & \frac{1}{2}\frac{EA}{L} & -\frac{2EA}{L} \\ \frac{5}{2}\frac{EA}{L} & -\frac{1}{2}\frac{EA}{L} & \frac{5}{2}\frac{EA}{L} + \frac{d}{4} & -\frac{1}{2}\frac{EA}{L} + \frac{d}{4} & 0 \\ -\frac{1}{2}\frac{EA}{L} & \frac{1}{2}\frac{EA}{L} & -\frac{1}{2}\frac{EA}{L} + \frac{d}{4} & \frac{1}{2}\frac{EA}{L} + \frac{d}{4} & 0 \\ 0 & -\frac{2EA}{L} & 0 & 0 & \frac{2EA}{L} \end{bmatrix} \begin{bmatrix} a_2 \\ a_3 \\ b_2 \\ b_3 \\ a_4 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ P \end{bmatrix} $$

and its solution yields

$$ \begin{bmatrix} a_2 & a_3 & b_2 & b_3 & a_4 \end{bmatrix} = \begin{bmatrix} \frac{1}{2} PL & + \frac{P}{d} & \frac{3}{2} PL & + \frac{P}{d} & -\frac{P}{d} & -\frac{P}{d} & \frac{2PL}{EA} + \frac{P}{d} \end{bmatrix}. \quad [32] $$

For interface elements, the displacement jump across the discontinuity $S_d$ is expressed as the difference $[u_h] = a_3 - a_2 = P/d$ of the displacement of the doubled nodes. In PU-based discontinuous elements, the displacement jump is expressed through [10] as $[u_h] = Nb|_{S_d} = 0.5 \times (b_2 + b_3) = -P/d$ (see Figure 10). Note that the minus sign in the displacement jump value indicates that the local frame of the discontinuity is in the opposite direction with respect to the global frame $x$—in other words, the negative sign is due to the choice of the positive part of the domain.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Occurrence of oscillations in the traction profile for some values of the interface stiffness and with Gauss integration scheme for the interface contribution (from Simone (2004))</figcaption></figure>

**9. PU-based discontinuous elements and the “dummy stiffness” approach**

The “dummy stiffness” approach is often used in combination with conventional interface elements to simulate perfect contact prior to the loss of material coherence along potential crack lines (like *e.g.* the dotted line in Figure 11a). This approach consists in considering a high value of the interface stiffness and is known to produce unrealistic oscillations in the traction profile like those shown in Figure 11b (Schellekens *et al.*, 1993). These oscillations can however be eliminated by an appropriate choice of the interface integration scheme.

PU-based discontinuous elements, when used with a “dummy stiffness” approach, suffer from the same type of problem (Simone, 2004). Similar to classical interface elements, oscillations in the traction profile are linked to the appearance of pathological coupling between degrees of freedom in the part of the stiffness matrix responsible for the “interface” behavior (matrix $K_{bb}$ in [23]b). Unlike interface elements, these oscillations are not removed by using a particular integration scheme; they can be however circumvented by activating the degrees of freedom responsible for the displacement jump when they are required. This problem may still be present if favorable conditions exist, *e.g.* high stress gradient and high discontinuity stiffness in elastic regime or cyclic loading conditions which promote crack closure. In these cases, it is worth considering the solution proposed by Mourad *et al.* (2006).

**10. Crack propagation in a single-edge notched beam**

A single-edge notched beam, depicted in Figure 12, is subjected to an antisymmetric four-point-shear loading (Schlangen, 1993) which results in a curved crack path running from the lower right part of the notch towards a point to the right of the lower-
<figure><img src="image_14.png" /><figcaption>**Figure 12.** *Single-edge notched beam (depth = 100 mm; all dimensions in mm).*</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** *Curved crack path in the SEN beam (from Wells et al. (2001))*</figcaption></figure>

right support. The beam is in concrete with Young's modulus $E = 35$ GPa, Poisson's ratio $\nu = 0.2$, tensile strength $f_t = 3.0$ MPa and fracture energy $G_f = 0.1$ N/mm². The material parameters have been taken from Wells *et al.* (2001). The loading plats have a Young's modulus one order of magnitude larger than the concrete. The load is applied by means of an indirect displacement control procedure with the crack mouth sliding displacement (cmsd), defined as the relative vertical displacement of the opposite faces of the notch, taken as control parameter. The boundary conditions are specified by constraining the displacement at the upper-right support in both directions and by constraining the vertical displacement at the upper-left support. The beam is analyzed under plane strain conditions. The simulation is performed using 1184 six node triangular elements. Crack propagation has been obtained using a principal stress initiation criterion.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Deformation of the SEN beam and $\sigma_{xx}$ stresses. The displacements are amplified by a factor 100 (from Remmers (2006))</figcaption></figure>

In this application the bulk material is considered as linear elastic whereas the normal traction-separation law at the discontinuity is defined by

$$t_n = f_t \exp \left( -\frac{f_t}{G_f} \kappa \right), \qquad [33]$$

where $f_t$ is the tensile strength of the material and $G_f$ is the fracture energy. $\kappa$ is a history parameter equal to the largest value of the normal separation at the interface. Shear traction contributions can be incorporated in the formulation if necessary. In this example Wells et al. (2001) used

$$t_s = d_{int} \exp (h_s \kappa) [[u]]_s \qquad [34]$$

with $h_s = \ln (\frac{d_{\kappa=1.0}}{d_{int}})$, $d_{int}$ the initial crack shear stiffness, $[[u]]_s$ the crack sliding displacement, and $d_{\kappa=1.0}$ the crack shear stiffness when $\kappa = 1.0$. The material tangent $T (\dot{t} = T [[\dot{u}]])$ can be expressed by differentiating the expressions of the normal and tangential traction:

$$\begin{bmatrix} \dot{t}_n \\ \dot{t}_s \end{bmatrix} = \begin{bmatrix} -\frac{f_t^2}{G_f} \exp \left( -\frac{f_t}{G_f} \kappa \right) & 0 \\ h_s d_{int} \exp (h_s \kappa) [[\dot{u}]]_s & d_{int} \exp (h_s \kappa) \end{bmatrix} \begin{bmatrix} [[\dot{u}]]_n \\ [[\dot{u}]]_s \end{bmatrix}. \qquad [35]$$

The predicted crack path shown in Figures 13 and 14 is curved and is in agreement with the experimental results (Schlangen, 1993). Worth noting is that the crack is independent of the mesh (see close-up of crack path in Figure 13).

In terms of global quantities, the numerical load-cmsd curves are compared to the experimental response in Figure 15. The agreement is quite good for both constant and variable shear stiffness up to the first part of the post-peak branch; the numerical curves are slightly too brittle for cmsd > 0.3 mm indicating that a more sophisticated expression of the traction-separation law is needed.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** *Experimental (Schlangen, 1993) versus numerical (Wells et al., 2001) load-crack mouth sliding displacement (cmsd) for variable shear stiffness (d<sub>int</sub> = 1.0 × 10<sup>2</sup> N/mm, d<sub>κ=1.0</sub> = 1.0 × 10<sup>-4</sup> N/mm) and constant shear stiffness (d<sub>int</sub> = 1.0 N/mm)*</figcaption></figure>

**11. Concluding remarks**

Partition of unity-based discontinuous elements are a powerful alternative to
conventional interface elements in problems with evolving boundaries. One of the
obvious advantages, compared to conventional interface elements, is that mesh ad-
justment along the propagating boundary is avoided since a discontinuity can pass
arbitrarily through solid elements. This extremely flexible and powerful approach has
been applied to a wide variety of problems which could not be described here. For
more information on the partition of unity method and its extension to discontinuous
element representation, the interested reader is referred to the publications that refe-
rence Duarte *et al.* (Duarte *et al.*, 1996b; Duarte *et al.*, 2000; Duarte *et al.*, 1996a; Oden
*et al.*, 1998), Melenk *et al.* (1996), and Moës *et al.* (1999).

As a final remark, it bears emphasis that XFEM, PUFEM and GFEM are examples
of the so-called partition of unity method which originated in the works of Ba-
buška et al. (Babuška et al., 1994; Melenk et al., 1996) and Duarte and Oden (Duarte
et al., 1995; Duarte et al., 1996b).

NOTE. — This paper is mainly based on Duarte *et al.* (2007) and Simone (2004).

**12. References**

Alfano G., Crisfield M. A., « Finite element interface models for the delamination analysis of laminated composites: Mechanical and computational issues », *International Journal for Numerical Methods in Engineering*, vol. 50, n° 7, 2001, p. 1701-1736.

Babuška I., Caloz G., Osborn J. E., « Special finite element methods for a class of second order elliptic problems with rough coefficients », *SIAM Journal on Numerical Analysis*, vol. 31, n° 4, 1994, p. 945-981.
Babuška I., Melenk J. M., « The partition of unity finite element method », *International Journal for Numerical Methods in Engineering*, vol. 40, 1997, p. 727-758.

Béchet E., Minnebo H., Moës N., Burgardt B., « Improved implementation and robustness study of the X-FEM for stress analysis around cracks », *International Journal for Numerical Methods in Engineering*, vol. 64, 2005, p. 1033-1056.

Belytschko T., Black T., « Elastic crack growth in finite elements with minimal remeshing », *International Journal for Numerical Methods in Engineering*, vol. 45, n° 5, 1999, p. 601-620.

Belytschko T., Chen H., Xu J., Zi G., « Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment », *International Journal for Numerical Methods in Engineering*, vol. 58, n° 12, 2003, p. 1873-1905.

Bittencourt T. N., Wawrzynek P. A., Ingraffea A. R., Souza J. L., « Quasi-automatic simulation of crack propagation for 2D LEFM problems », *Engineering Fracture Mechanics*, vol. 52, n° 2, 1996, p. 321-334.

Broberg K. B., *Cracks and Fracture*, Academic Press, San Diego, 1999.

Chessa J., Wang H., Belytschko T., « On the construction of blending elements for local partition of unity enriched finite elements », *International Journal for Numerical Methods in Engineering*, vol. 57, 2003, p. 1015-1038.

Comi C., Mariani S., Perego U., « An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation », *International Journal for Numerical and Analytical Methods in Geomechanics*, vol. 31, 2007, p. 213-238.

Dolbow J., Moës N., Belytschko T., « Discontinuous enrichment in finite elements with a partition of unity method », *Finite Elements in Analysis and Design*, vol. 36, 2000, p. 235-260.

Duarte C. A., Babuška I., « Mesh-independent directional p-enrichment using the generalized finite element method », *International Journal for Numerical Methods in Engineering*, vol. 55, n° 12, 2002, p. 1477-1492.

Duarte C. A., Babuška I., Oden J. T., « Generalized Finite Element Methods for Three Dimensional Structural Mechanics Problems », *Computers and Structures*, vol. 77, 2000, p. 215-232.

Duarte C. A. M., Oden J. T., « *Hp Clouds*—A meshless method to solve boundary-value problems », Technical Report 95-06, May, 1995. TICAM, University of Texas at Austin, TX, USA.

Duarte C. A. M., Oden J. T., « An *hp* adaptive method using clouds », *Computer Methods in Applied Mechanics and Engineering*, vol. 139, 1996a, p. 237-262.

Duarte C. A. M., Oden J. T., « *Hp Clouds* – An *hp* Meshless Method », *Numerical Methods for Partial Differential Equations*, vol. 12, n° 6, 1996b, p. 673-705.

Duarte C. A., Reno L. G., Simone A., « A high-order generalized FEM for through-the-thickness branched cracks », *International Journal for Numerical Methods in Engineering*. In the press, http://dx.doi.org/10.1002/nme.2012.

Duarte C., Kim D.-J., « Analysis and Applications of a Generalized Finite Element Method with Global-Local Enrichment Functions », *Computer Methods in Applied Mechanics and Engineering*. Submitted.
Gasser T. C., Holzapfel G. A., « Modeling 3D crack propagation in unreinforced concrete using PUFEM », *Computer Methods in Applied Mechanics and Engineering*, vol. 194, n° 25-26, 2005, p. 2859-2896.

Goodman R. E., Taylor R. L., Brekke T. L., « A model for the mechanics of jointed rocks », *Journal of the Soil Mechanics and Foundations Division*, vol. 94, 1968, p. 637-659.

Jirásek M., « Comparative study on finite elements with embedded cracks », *Computer Methods in Applied Mechanics and Engineering*, vol. 188, 2000, p. 307-330.

Jirásek M., Zimmermann T., « Embedded crack model. Part II: Combination with smeared cracks », *International Journal for Numerical Methods in Engineering*, vol. 50, 2001, p. 1291-1305.

Laborde P., Pommier J., Renard Y., Salaün M., « High-order extended finite element method for cracked domains », *International Journal for Numerical Methods in Engineering*, vol. 64, 2005, p. 354-381.

Mariani S., Perego U., « Extended finite element method for quasi-brittle fracture », *International Journal for Numerical Methods in Engineering*, vol. 58, n° 1, 2003, p. 103-126.

Melenk J. M., Babuška I., « The partition of unity finite element method: Basic theory and applications », *Computer Methods in Applied Mechanics and Engineering*, vol. 139, n° 1-4, 1996, p. 289-314.

Mittelstedt C., Becker W., « Semi-analytical computation of 3D stress singularities in linear elasticity », *Communications in Numerical Methods in Engineering*, vol. 21, 2005, p. 247-257.

Moës N., Belytschko T., « Extended finite element method for cohesive crack growth », *Engineering Fracture Mechanics*, vol. 69, 2002, p. 813-833.

Moës N., Dolbow J., Belytschko T., « A finite element method for crack growth without remeshing », *International Journal for Numerical Methods in Engineering*, vol. 46, n° 1, 1999, p. 131-150.

Mourad H. M., Dolbow J., Harari I., « A bubble-stabilized finite element method for Dirichlet constraints on embedded interfaces », *International Journal for Numerical Methods in Engineering*, vol. 69, n° 4, 2006, p. 772-793.

Oden J. T., Duarte C. A., « Clouds, cracks and FEM's », in B. D. Reddy (ed.), *Recent Developments in Computational and Applied Mechanics*, International Center for Numerical Methods in Engineering, CIMNE, Barcelona, Spain, 1997, p. 302-321.

Oden J. T., Duarte C. A., Zienkiewicz O. C., « A new cloud-based hp finite element method », *Computer Methods in Applied Mechanics and Engineering*, vol. 153, n° 1-2, 1998, p. 117-126.

Remmers J. J. C., Discontinuities in Materials and Structures, PhD thesis, Delft University of Technology, 2006.

Saouma V. E., Brühwiler E., Boggs H. L., « A review of fracture mechanics applied to concrete dams », *International Journal of Dam Engineering*, vol. 1, n° 1, 1990, p. 41-57.

Schellekens J. C. J., de Borst R., « On the numerical integration of interface elements », *International Journal for Numerical Methods in Engineering*, vol. 36, 1993, p. 43-66.

Schlangen E., Experimental and Numerical Analysis of Fracture Processes in Concrete, PhD thesis, Delft University of Technology, 1993. Also published in HERON, 38(2):3-117, 1993.
Simo J. C., Oliver J., Armero F., « An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids », *Computational Mechanics*, vol. 12, 1993, p. 277-296.

Simone A., « Partition of unity-based discontinuous elements for interface phenomena: Computational issues », *Communications in Numerical Methods in Engineering*, vol. 20, n° 6, 2004, p. 465-478.

Simone A., Duarte C. A., der Giessen E. V., « A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries », *International Journal for Numerical Methods in Engineering*, vol. 67, n° 8, 2006, p. 1122-1145.

Simone A., Sluys L. J., « The use of displacement discontinuities in a rate-dependent medium », *Computer Methods in Applied Mechanics and Engineering*, vol. 193, n° 27-29, 2004, p. 3015-3033.

Simone A., Wells G. N., Sluys L. J., « From continuous to discontinuous failure in a gradient-enhanced continuum damage model », *Computer Methods in Applied Mechanics and Engineering*, vol. 192, n° 41-42, 2003, p. 4581-4607.

Stazi F. L., Budyn E., Chessa J., Belytschko T., « An extended finite element method with higher-order elements for curved cracks », *Computational Mechanics*, vol. 31, 2003, p. 38-48.

Sukumar N., Prévost J.-H., « Modeling quasi-static crack growth with the extended finite element method - Part I: Computer implementation », *International Journal of Solids and Structures*, vol. 40, 2003, p. 7513-7537.

Szabó B. A., « Mesh design for the *p*-version of the finite element method », *Computer Methods in Applied Mechanics and Engineering*, vol. 55, 1986, p. 181-197.

Taylor R. L., Zienkiewicz O. C., Oñate E., « A hierarchical finite element method based on the partition of unity », *Computer Methods in Applied Mechanics and Engineering*, vol. 152, n° 1-2, 1998, p. 73-84.

Ventura G., « On the elimination of quadrature subcells for discontinuous functions in the eXtended Finite-Element Method », *International Journal for Numerical Methods in Engineering*, vol. 66, n° 5, 2006, p. 761-795.

Wells G. N., Sluys L. J., « A new method for modelling cohesive cracks using finite elements », *International Journal for Numerical Methods in Engineering*, vol. 50, n° 12, 2001, p. 2667-2682.

Xiao Q. Z., Karihaloo B. L., « Improving the accuracy of XFEM crack tip fields using higher order quadrature and statically admissible stress recovery », *International Journal for Numerical Methods in Engineering*, vol. 66, 2006, p. 1378-1410.