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

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.


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.

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 domainV of a boundary value problem: n I=1 N I (x) = 1 ∀x ∈V . [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 where the subscript "h" denotes the approximate field,q I (termed "regular" or "standard" degrees of freedom (dofs)) are nodal values related to the basis N I , andq 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 : [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 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 where N is a l × (l × n) matrix containing the standard finite element shape functions, N E is a (l×n)×(l×m×n) matrix containing the extra basis terms for the enrichment, q is a (l × n) × 1 vector containing standard degrees of freedom andq is a (l × m × n) × 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 freedomq when needed (a hierarchical finite element formulation based on the partition of unity method).

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.

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 bodyV 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 + ∪ V − ). Displacementsū are prescribed on S u , while tractionsT are prescribed on S t . The displacement field can be decomposed as where H(x) is the Heaviside function centered at the discontinuity surface S d (H = 1 if x ∈ V + , H = 0 if x ∈ V − ) andû andũ are continuous functions onV . The discontinuity is introduced by the Heaviside function H at the discontinuity surface S d and the value of the displacement jump u across the discontinuity surface is given byũ 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 where (·) s refers to the symmetric part of (·). The strain field is defined everywhere in the body except at the discontinuity surface.  Figure 2. BodyV crossed by a discontinuity S d

Describing a displacement discontinuity
In describing a displacement discontinuity, and considering [5] now rewritten for convenience as if the standard displacement field is interpolated by the regular interpolation N a, the displacement jump can be described through the enrichment N N E b = HN Hb (Moës et al., 1999), where 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 × n) × 1 vector of extra degrees of freedom b (same dimensions as the vector of standard degrees of freedom a) and in the (l × n) × (l × n) matrix HH. 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 N H. Hence, the discretized format of [6] reads, for points located in an element with enriched nodes, 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 [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 = N a 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: [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)). crack standard nodes, set I discontinuous enrichment, set J near-tip enrichment, set K r θ

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 where (r, θ) 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) 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;Mittelstedt et al., 2005). Complex crack front geometries, curved crack surfaces or the intersection 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 errors 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 resolving the non-smooth displacement field near a crack tip/front by local mesh refinements combined with high-order approximations (Szabó, 1986) (hp discretization). This approach avoids the difficulties and limitations related to near crack front enrichment 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.

High-order approximations
At variance with the established enrichment of the displacement field with discontinuous 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), , 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 computations. 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 twodimensions. 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û h andũ h , respectively. 1) p = 1, q = 1: the shape functions at a node x I are given by 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 [15] These shape functions, of course, can not represent the discontinuous components of the solutionũ h .
3) p = 2, q = 2: the shape functions at a node x I are given by [16] 4) p = 2, q = 1: the shape functions at a node x I are given by [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).

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 bodyV without body forces depicted in Figure 2 can be summarized by where σ 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 whereū 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ŵ u andw 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 with D e the tangent matrix for the bulk material, and the traction rate at a discontinuitẏ where T relates traction rateṫ and displacement jump rate 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)) where In the above relations, B is a matrix containing derivatives of shape functions, σ is an array containing the components of the stress tensor in engineering notation and 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).

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-discontinuity

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)
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 . 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 interaction 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 discontinuity extension direction is aligned with the previous discontinuity segment to avoid an incorrect direction determination due to the bias produced by an unsymmetric domain 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.

A minimal one-dimensional example
Some of the issues related to the implementation of PU-based discontinuous elements are illustrated by means of a simple one-dimensional example of a bar in tendiscontinuity tip 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.

Stiffness matrix computation.
The stiffness matrix of the one-dimensional truss element of length L depicted in Figure 9a is 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 truss . The remaining sub-matrices are expanded as In K bb , the volume integral equals K ab while the surface integral is evaluated on x = L/2 and yields Assembly of the sub-matrices into the element stiffness matrix yields

Figure 9. One-dimensional (a) truss element and (b) PU-based discontinuous truss element
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 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).

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 which yields and its solution yields 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 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.

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

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-  right support. The beam is in concrete with Young's modulus E = 35 GPa, Poisson's ratio ν = 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 platens 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.  (2006)) In this application the bulk material is considered as linear elastic whereas the normal traction-separation law at the discontinuity is defined by where f t is the tensile strength of the material and G f is the fracture energy. κ is a history parameter equal to the largest value of the normal separation at the interface. with h s = ln d κ=1.0 d int , d int the initial crack shear stiffness, u s the crack sliding displacement, and d κ=1.0 the crack shear stiffness when κ = 1.0. The material tangent T (ṫ = T u ) can be expressed by differentiating the expressions of the normal and tangential traction: 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 10. Applied load plotted against crack mouth sliding displacement (cmsd) for experimental results Figure 15. Experimental (Schlangen, 1993) versus numerical (Wells et al., 2001) load-crack mouth sliding displacement (cmsd) for variable shear stiffness (d int = 1.0 × 10 2 N/mm, d κ=1.0 = 1.0 × 10 −4 N/mm) and constant shear stiffness (d int = 1.0 N/mm)

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 adjustment 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 reference 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 Babuš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).