The Discontinuity-Enriched Finite Element Method

We introduce a new methodology for modeling problems with both weak and strong discontinuities independently of the finite element discretization. At variance with the eXtended/Generalized Finite Element Method (X/GFEM), the new method, named the Discontinuity-Enriched Finite Element Method (DE-FEM), adds enriched degrees of freedom only to nodes created at the intersection between a discontinuity and edges of elements in the mesh. Although general, the method is demonstrated in the context of fracture mechanics, and its versatility is illustrated with a set of traction-free and cohesive crack examples. We show that DE-FEM recovers the same rate of convergence as the standard FEM with matching meshes, and we also compare the new approach to X/GFEM.


INTRODUCTION
This article is concerned with problems whose solution is C 0 -and C −1 -continuous, i.e., problems containing weak and strong discontinuities, respectively. The eXtended/Generalized Finite Element Method (X/GFEM) has revolutionized the modeling of such discontinuities. The most significant asset of X/GFEM is the capability to decouple the representation of a discontinuity from the underlying finite element discretization. Here, we introduce a novel discontinuity-enriched finite element formulation that has several advantages with respect to X/GFEM. This approach promises a new modeling paradigm by placing enriched nodes only along discontinuities. This has important implications in terms of computer implementation and, in particular, in the imposition of Dirichlet boundary conditions. The standard Finite Element Method (FEM) has been traditionally used to model weak and strong discontinuities by designing conforming or matching meshes in which the edges of finite elements align to such discontinuities. For weak discontinuities, the jump in the field gradient is readily accommodated by the C 0 continuity at finite element edges. Strong discontinuities such as cracks require the definition of so-called double nodes along edges of finite elements aligned with the crack geometry. To date, the creation of matching discretizations remains an error-prone process that can even result in an inadequate mesh because finite elements with bad aspect ratios can be a source of large numerical errors. Creating matching meshes also requires the use of robust (and expensive) software. Furthermore, modeling problems with evolving geometries or discontinuities renders the cost of creating such a matching mesh prohibitive (specially in 3-D simulations).
X/GFEM 1,2 tackles such problems in an elegant manner as discontinuities need not conform to the finite element mesh. This flexibility in the analysis, however, comes at a price because the complexity of the conformal mesh generation is transferred

FORMULATION
The method is first demonstrated by means of a one-dimensional problem whose solution is contrasted to that obtained with X/GFEM. A formal formulation for higher dimensions is given in Section 2.2. Figure 1 shows a 1-D bar of length L clamped at the left end and subjected to a load P at the right end. The bar is composed by two parts, Ω I and Ω II , that are connected at x = x Γ by an elastic spring of stiffness k. In addition to the strong discontinuity at x Γ , the mismatch in Young's modulus between the parts also introduces a weak discontinuity. The solution of the corresponding boundary value problem yields the displacement at the free end expressed as the sum of the elastic deformation contributions of the parts and the spring:

A 1-D example with both weak and strong discontinuities
where w = x Γ ∕L is the location of the discontinuity relative to the length of the bar. We now solve the problem with a single finite element, denoted Ω, employing a solution of the form One may immediately recognize that Equation (2) shows a structure similar to that of X/GFEM. The first term corresponds to the standard FEM component of the approximation, where the Lagrange shape functions N j and the standard degrees of freedom (DOFs) U j are associated with the end nodes (those identified with the symbol in Figure 2A). The other terms enrich the approximation. The functions and are enrichment functions for a weak and a strong discontinuity, respectively, and and are their corresponding enriched DOFs. These enrichment functions are created from Lagrange shape functions in sub-domains Ω I , Ω II (see Figure 2B). In Equation (2), is the IGFEM enrichment function that is used to introduce a jump in the gradient of the field. 8,10 For this simple example, we will just use Therefore, the weak enrichment attains a maximum value of one at x Γ , and it ramps linearly to zero at the ends of the bar. Note that the IGFEM weak enrichment uses a scale factor to improve the conditioning as the interface is moved closer to a node; 10 to ease the presentation, such scaling is omitted here. The information about the discontinuous field, here understood as a kinematic enrichment, is incorporated in . To construct this enrichment, we seek to satisfy the following two conditions: 1. In the spirit of IGFEM, 8 the enrichment function is created through a linear combination of Lagrange shape functions in sub-domains or integration elements, i.e., elements Ω I , Ω II created to the left and right side of the discontinuity; and 2. we require a jump of unit magnitude, i.e., [[ ] By ensuring this condition, the enriched DOF in Equation (2) represents the jump in the displacement field at x = x Γ , thus attributing an immediate physical interpretation to the enriched DOF. The first condition can be expressed as where w 1 and w 2 are still unknown weights, and, as shown in Figure 2B, To satisfy the second condition, we remark that the Lagrange shape functions N j form a partition of unity, i.e., ∑ j N j = 1. By noting that |∓N 1 − (±N 2 )| = 1, ∀x ∈ Ω, the weights are obtained from Equation (5) yields w 1 = −x Γ ∕L = −w while Equation (6) yields w 2 = 1 + w 1 = 1 − w. Thus, the enrichment function shown in Figure 2C takes the final form Equation (7) can be interpreted as having the Lagrange shape functions (associated with the discontinuity location) in integration elements being weighted by the location of the discontinuity relative to the length. Observe that, even though Equation (7) seems to degenerate as x Γ → 0 or x Γ → L, the enrichment function is bounded. In fact, noting that w = x Γ ∕L, we find that Having determined the discontinuous enrichment, we can now proceed with the solution of this problem following standard procedures.
, the stiffness matrix and force vector are given by respectively, where ⟦E⟧ ≡ E II − E I . It is worth noting that both weak and strong enrichment functions are equal to zero at the left end; therefore, imposing the essential boundary condition, i.e., u (0) = 0, is as straightforward as in the standard FEM. The system KU = F has solution and thus, the approximation recovers the exact displacement field We can infer the following by inspecting Equation (9): with an inactive weak enrichment; • the perfectly bonded case is recovered when k → ∞, the solution becomes with an inactive strong enrichment; and • when both conditions E II = E I = E and k → ∞ are met (perfectly bonded parts of the same material), the solution vector is and both weak and strong enrichments are inactive.
We would like to stress that an independent linear field at each side of x Γ can only be obtained by using both weak and strong enrichments. In other words, using the strong discontinuity enrichment alone would fail to produce two independent linear fields because d ∕dx is the same at either side of the discontinuity.
We now turn our attention to the solution obtained with X/GFEM, where we use the Heaviside step function in Figure 2D, as enrichment. This enrichment, equal to zero at the left end of the bar, has been chosen to simplify the enforcement of the essential boundary condition; therefore, as in standard FEM, this relates to the standard DOF only. A step function with a value different from zero at the clamped boundary would have required the enforcement of a multi freedom constraint through master-slave elimination, Lagrange multiplier adjunction, or penalty augmentation, 17,18 as the displacement field is now a function of both standard and enriched DOFs. This is at variance with DE-FEM, where the enrichment function is always zero at standard mesh nodes; this implies that standard FEM procedures can always be used in DE-FEM to enforce Dirichlet boundary conditions. The bar is again discretized with a single finite element, and the vector of nodal unknowns contains the standard DOFs followed by those corresponding to the enrichment function. The stiffness matrix and the force vector are expressed as and The system KU = F has solution yielding the same displacement field as in Equation (10).

Formulation in higher dimensions
Consider the Euclidean vector space R d (d = 2, 3) spanned by an orthonormal basis {e i }. A solid body (cf. Figure 3   Γ is, in turn, subdivided into two disjoint regions, Γ u ≠ ∅ and Γ t , where essential (Dirichlet) and natural (Neumann) boundary conditions are prescribed, respectively. Consider a crack Γ c ⊆ Γ t parameterized by a coordinate s along which a displacement jump ≡ ⟦u⟧ ∀x ∈ Γ c is defined. The two regions adjacent to the crack are conveniently labeled as positive and negative as shown in Figure 3. To uniquely label these domains, the crack is endowed with a direction given by the unit vector n c normal to the crack. Given a point p ∈ Ω, let y = arg min x∈Γ c ‖x − p‖ be the closest point on the crack. Point p lies on the positive side if n c · (p − y) > 0, where n c is evaluated at point y.
Consider now the vector-valued function spaces where H 1 (Ω) denotes the first-order Sobolev function space on Ω. Takingũ ∈  ∶ũ| Γ u =ū, we define the linear variety  ⋆ ≡ũ +  0 so that every element in  ⋆ satisfies the Dirichlet boundary condition. The weak form of the static boundary value problem for a deformable cracked body is expressed as: Given the prescribed displacementū ∶ Γ u → R d , the prescribed tractiont ∶ Γ t ∖Γ c → R d , the cohesive traction t c ∶ Γ c → R d , the Cauchy stress ∶ Ω × Ω → R d , and the body force b ∶ Ω → R d , find the displacement field u ∈  ⋆ such that where ≡ 1 2 ( v + v ⊺ ). To solve this equation, the domain is discretized into a series of finite elements Ω i that are not conforming to the crack, such that Ω h = int The functions u h that are used in the finite-dimensional form of Equation (19) have to be chosen from a set that allows for a displacement jump at Γ c . In DE-FEM, the displacement field takes the form In the first term, h is the index set of standard FEM nodes, N i denotes the ith Lagrange shape function, and U i is the vector of corresponding DOFs. In the second and third terms, which represent the enriched portion of the approximation, w and s denote index sets of nodes associated with weak ( i ) and strong ( i ) enriched DOFs. The corresponding enrichment functions for weak and strong discontinuities are i and i , respectively. These enrichment functions are constructed by the combination of Lagrange shape functions in elements (the so-called integration elements) that are created by partitioning the mesh element intersected by a discontinuity. These functions are local to their parent elements by construction, and thus the multiplication by the partition of unity of the underlying mesh is not required as it is done with X/GFEM. Figure    nodes created along the crack instead to those of the original mesh. In what follows, we will focus on the construction of the strong discontinuity enrichment, although a similar reasoning can be applied to the construction of weak enrichments; the latter follows previous studies 8, 10 and will not be reported here.
For the construction of the strong discontinuity enrichment for 2-D elements, we consider the three-node triangular element shown in Figure 5. The strong discontinuity Γ c intersects the element edges at locations where we create the new nodes i and j, to which we associate enriched DOFs i and j , respectively. We identify the node with a plus (minus) sign superscript when it is considered as part of the positive (negative) side of the discontinuity. These newly created nodes are used to identify the integration elements Ω I , Ω II , and Ω III shown in Figure 5. The objective is now to determine the strong enrichment function by obtaining the weights that are associated with nodes i and j. Following the same procedure described in Section 2.1, the equations are used to determine w i1 , w i2 . Next, we define w i = −w i1 so that w i2 = 1 − w i . It is worth noting that the weight is bounded, i.e., 0 ⩽ w i ⩽ 1. It can be shown that obtaining the weight w i in this way is equivalent to Figure 5) are the lengths of the sides of the integration elements Ω II and Ω I , respectively. Once the weight is determined, the discontinuous enrichment function associated with i is Similarly, for the jth node, the equations and (24) determine w j1 , w j2 so that w j = −w j1 and 0 ⩽ w j ⩽ 1. As before, the weight can be simply determined as w j = a j ∕(a j + b j ), where this time a j and b j correspond to the sides of integration elements Ω III and Ω I , respectively. The discontinuous enrichment function associated with j is Appendix D presents a Mathematica ® script that computes the weights and plots the discontinuous enrichment functions for the simple case shown in Figure 5. Figure 6 illustrates the creation of the discontinuous enrichment functions associated with enriched DOFs i and j for a bilinear quadrilateral element. The discontinuity splits the element into integration elements Ω I , Ω II that are connected to newly created nodes i, j. Following the procedure described above, the strong discontinuity enrichment functions i and j associated with these nodes are given by where the weights w i and w j are simply determined by finding the distance to the discontinuity relative to the corresponding edge length, i.e., w i = a i ∕(a i + b i ) and w j = a j ∕(a j + b j ) (cf. Figure 6). Note that obtaining the weights requires consistency insofar as to which nodes are used to compute the distance to the discontinuity. In other words, only nodes lying on one side of the crack (whichever that side is) are used for the computation of the weights.
In the foregoing discussion, integration elements are also quadrilateral elements. There are cases, however, when a discontinuity traverses a quadrilateral element in such a way that at least one triangular integration element is needed. It is then important to note that quadrilateral elements can also be split in triangular integration elements as long as the resulting enrichment functions are able to recover the bilinear behavior of the parent element, i.e., the xy term. This can be easily accomplished by using six-node quadratic triangles as integration elements. In either case, the resulting enrichment functions in DE-FEM ramp to zero at nodes of the original mesh, and thus, they retain local support by construction. Figure 7 shows the creation of integration elements when dealing with a crack tip j. In this situation, the crack tip is endowed with just a weak discontinuity enrichment j because, by definition, there is no jump in the displacement at a crack tip. Retaining the full approximation, with strong and weak enrichment, would require to prescribe j = 0. This and other computer implementation details are discussed in Appendix A.  Following standard procedures, the local system k i u i = f i for the ith element is expressed as where the subscript u refers to the standard part of the approximation, subscripts and indicate the weak and strong enriched contributions, respectively, and In addition, D represents the material constitutive matrix, here assumed linear elastic. Because of the form of the enrichment functions in DE-FEM, these integrals are evaluated numerically through Gauss quadrature on the interior of integration elements where the functions involved are well behaved. The global stiffness matrix and force vector are obtained by considering the contribution of all elements in the discretization. In other words, denoting with A the finite element assembly operator, K = A i k i and F = A i f i .

Cohesive formulation
When an element is crossed by a strong discontinuity, a corresponding cohesive element is defined with the purpose of simplifying the generation of its stiffness matrix contribution and assembly procedures. In the rest of the paper, we therefore use the term elements to identify those used to discretize the continuum and cohesive elements for those that identify the portion of the discontinuity defined by the intersection of an element and the discontinuity itself. Clearly, a cohesive discontinuity is defined as the union of all its cohesive elements.
In the discretized form of Equation (19), the cohesive traction contribution is obtained by the individual contribution of the cohesive elements, i.e., Π c = ∑ e Π e c . In the following, all quantities are assumed at the element level, and thus, we omit the superscript e. Figure 8 shows element Ω i crossed by a cohesive crack. The position along the discontinuity segment, defined by points i and j, is described by means of a master coordinate following the traditional isoparametric formulation, i.e., x = x ( ) , −1 ⩽ ⩽ 1. The displacement jump at the discontinuity is expressed as The jump along the discontinuity is interpolated by using standard hat shape functions and the values of the strong enriched DOFs as where N = , I is the d × d identity matrix, and ⊗ denotes the Kronecker product. Noteworthy, the displacement jump at the discontinuity is determined by considering quantities that are local to the discontinuity: the hat shape functions and the enriched DOFs. This fact has major consequences in the implementation of the cohesive formulation. Indeed, the evaluation of the cohesive element contributions can be done by only integrating along the line elements that define the discontinuity, as shown by the cohesive contributions k c to the tangent stiffness matrix k and f c to the force vector f : where t c is the vector of cohesive tractions that depends on the cohesive law. As before, the contribution of all cohesive elements is assembled as A i k c,i and A i f c,i for the global stiffness matrix and force vector, respectively. Noteworthy, the computer implementation is the same as that of traction-free cracks, the difference only being in the creation of cohesive elements when cohesive cracks are present.

Comparison between DE-FEM and X/GFEM
Important differences between X/GFEM and DE-FEM are listed below.

Stress intensity factors
In some of the examples found later in this manuscript, DE-FEM is used to extract stress intensity factors (SIFs). Mixed-mode SIFs are obtained following the formulation by Yau et al. 21 : where K I , K II denote SIFs for modes I and II, respectively, the superscripts refer to two independent equilibrium states, and M (1,2) denotes the conservation integral considering their superposition. In this work, we use the domain representation of such integral 22 : where the terms are taken with respect to a coordinate system that is local to the crack tip (see Figure 9A). The weight function q 1 has unit magnitude in finite elements completely enclosed by the circle of radius , ramping down to zero at nodes of intersected elements (see Figure 9B). Note that elements fully enclosed by the circle do not contribute to the interaction integral since q 1 ∕ x j = 0. By choosing an appropriate equilibrium state, it is possible to obtain the mixed-mode SIFs from Equations (33) and (34). If, for example, we assume an auxiliary equilibrium state of pure mode I, then K

NUMERICAL EXAMPLES
Unless otherwise specified, in the following examples we use plane strain conditions and consistent units when these are not specified. Both weak and strong enrichment functions as specified by Equation (20) are used throughout to ensure independent linear fields at each side of the discontinuities. In addition, a quadrature rule that exactly integrates integration elements is also used. In the X/GFEM examples the two sub-elements generated when an element is cut by a discontinuity are triangulated, and in each triangular domain we used a three-point Gauss quadrature scheme.

A discontinuous patch test
This example shows the capability of the new formulation to generate two kinematically independent configurations when a domain is arbitrarily cut by a discontinuity. This will be demonstrated by means of a patch test. It bears emphasis that due to the simplicity in the application of essential boundary conditions in DE-FEM, the discontinuous patch test given here is simpler than that proposed for the X/GFEM by Dolbow and Devan. 23 Consider the square plate shown in Figure 10A with side L = 1. The plate is cut into two parts by a crack at x 2 = 0 (solid red line), and tractionst = t e 1 and 2t are applied to the right edge above (Γ + t ) and below (Γ − t ) the crack, respectively. A possible discretization of the problem is given in Figure 10B, where new nodes ( ) have been placed along the crack at the intersection with the element edges. It is worth noting that since both weak and strong discontinuity enrichments are considered, each of these nodes is associated with four DOFs (two for i and two for i ).
Given that the force is acting in the x 1 direction, the force vector acting on nodes 2, 3, and 15 can be readily obtained by performing the following 1-D integration where N = [ N 2 N 3 15 15 ] , and 15 has a unit value at node 15 and ramps linearly to zero at nodes 2 and 3 (no scaling is used here for the weak enrichment). Finally, it is worth noting that a node is also created along the constrained left edge of the plate. Because there is no displacement jump at x 9 , 9 = 0 and 9 can be readily obtained from Equation (20) as Prescribing 9 and 9 is as straightforward as in the standard FEM.
The resulting deformed configuration, together with the stress field, is reported in Figure 10C. It can be seen that the formulation passes this discontinuous patch test, i.e., both parts at either side of the crack have their independent constant stress field. It bears emphasis that, as expected, a zero value was computed for all enriched DOFs concerning the behavior in the x 2 direction. Unneeded enrichments are therefore not active.
Once again, we would like to stress the importance of using both weak and strong enrichment functions. In view of reducing computational requirements, one may be tempted to use only the strong discontinuity enrichment which, after all, is responsible for representing the jump in the displacement field. However, this approach would fail the discontinuous patch test because there are insufficient DOFs to represent independent linear fields at each side of the strong discontinuity.

DE-FEM convergence properties
In this example, we demonstrate the convergence properties of DE-FEM for classical 2-D plane strain mode I and mode II problems. Consider a square plate with length L = 2 discretized using a structured mesh of three-node triangles. A horizontal crack is defined from the left edge to the center of the plate. The displacement at the plate boundary is prescribed by using the crack-tip displacement fieldsū for modes I and II, as given in Appendix C in Equation (C1a) with K I = 1 and Equation (C1b) with K II = 1, respectively. The boundary condition u =ū is applied not only at all mesh nodes along the plate boundary but also to the discontinuity node created at the intersection between the plate boundary and the crack (at x = −e 1 ). Figure 11 shows the deformed configurations for both mode I and mode II obtained by the standard FEM and by DE-FEM. For the latter, the displacement field is accurately represented only after post-processing the displacement field in elements traversed by the crack.
The convergence behavior is now characterized by measuring the error in the energy norm, defined as

FIGURE 12
Convergence results for the error in the energy norm for (A) mode I and (B) mode II problems. Figure 12 shows this error as a function of the number of DOFs for modes I (A) and II (B). These figures show that DE-FEM converges at the same rate as that of the standard FEM with the use of matching meshes and X/GFEM with step enrichment alone. In the case of mode II, DE-FEM is even more accurate than standard FEM, and it is more accurate than X/GFEM with step enrichment in both cases. As mentioned in the previous example, the formulation does not pass the patch test if only the strong discontinuity enrichment is used. In terms of convergence, using the strong enrichment alone results in convergence rates of 0.21 and 0.22 for modes I and II, respectively.

Square plate with double-edge cracks
This example aims at recovering the mode I SIF for a plate containing two edge cracks for different values of the ratio a∕W (cf. Figure 13A). Irwin 24 reported that an approximation of the SIF for this problem can be expressed as K I = K 0 a , where .
(38)  Figure 13B and were obtained by using an unstructured mesh of constant-strain triangles discretizing half of the domain because of symmetry. A non-uniform mesh size was used in the discretization, varying from h∕W = 0.2 at x 2 = ±H to h∕W = 0.002 at x 2 = 0. The use of this fine mesh is needed because results are compared with the accurate technique of Yamamoto and Tokuda mentioned above, and the proposed enrichment functions do not capture the stress singularity at the crack tip. For the computation of the mode I SIF, a radius ∕W = 0.1 was used. The figure shows that DE-FEM is able to achieve a more accurate result than that of Bowie, but not as accurate as that of Yamamoto and Tokuda. Nevertheless, considering the fact that the latter mixes analytic and numerical results, and that DE-FEM does not make use of any singularity enrichment, the accuracy of the prediction is remarkable.

Square plate with centered slanted crack
In this example, we study the extraction of SIFs in mixed-mode loading. Consider a square plate of side 2W = 10 subjected to a traction ±te 2 at x 2 = ±W as shown in Figure 14. A crack of length 2a = 1 is situated at the center of the plate and inclined at an angle . The SIFs for this problem are given by 27

FIGURE 14
Square plate with centered crack slanted at an angle from the horizontal line. Results for this problem are given in Figure 15, where the legends show the type of mesh used and the deformed configuration for = ∕4. It can be seen that even coarse meshes give accurate values of the mixed-mode SIFs. The results for a fine unstructured finite element mesh lie on top of the theoretical values.

Beam with perfectly bonded cohesive interface
This example illustrates the use of DE-FEM with a cohesive crack. Consider a notched beam of length L = 450 mm, height h = 100 mm, and thickness b = 100 mm, with a notch a = 20 mm as illustrated in Figure 16A. The beam is composed by an elastic material with Young's modulus E = 2 × 10 4 N/mm 2 and Poisson's ratio = 0.2. Boundary conditions include a roller  Normal stiffness k n in N/mm 3 . support at its lower-left corner, a fixed support at its lower-right corner, and an applied tractiont = −1 kN/mm 2 e 2 . The beam is crossed by a discontinuity passing through the center line. A traction-free discontinuity represents the notch, while the rest of the discontinuity line is assigned increasing values of the normal stiffness, i.e., k n = { 2 × 10 3 , 2 × 10 4 , 2 × 10 5 } N/mm 3 , to simulate perfect material bonding. Rigid-body rotation is prevented by taking a tangential stiffness k t = 1 N/mm 3 . This problem was first studied by Rots 28 using interface elements and later by Simone 16 in the context of X/GFEM. Simone showed that, unless a structured finite element mesh is used, the cohesive tractions at the interface suffer from high oscillations close to the notch as the normal stiffness coefficient k n increases.
A structured mesh of bilinear quadrangles and an unstructured mesh of linear triangles, as shown in Figure 16B,C, are employed. Cohesive elements are integrated using a two-point Gauss-Lobatto rule. The results for this study are given in Figure 17, where the cohesive traction profile (in abscissas) is plotted against its position along the beam centerline (ordinates). As it was reported in the study of Simone, 16 the traction profile resulting from the X/GFEM simulation is smooth when dealing with a structured mesh of quadrangles ( Figure 17A), but the use of an unstructured triangular mesh produces oscillations in the cohesive traction profile as k n is increased ( Figure 17B). The same discretizations are used with DE-FEM, and while the structured quadrangular mesh gives roughly the same result as X/GFEM ( Figure 17C), it can be seen that the traction profile for the unstructured triangular mesh shows no oscillations ( Figure 17D).

CONCLUDING REMARKS
In this article, we have introduced DE-FEM, a novel enriched finite element formulation for the solution of problems with both weak and strong discontinuities. While emphasis has been placed on the solution of problems in fracture mechanics, the formulation is general and can therefore be applied to solve other partial differential equations where the field is C −1 -and/or C 0 -continuous. When modeling strong discontinuities, both weak and strong enrichment functions are used. We showed that the new formulation converges at the same rate as the standard FEM with meshes that conform to the discontinuities. Furthermore, we showed that DE-FEM can be used to extract SIFs with an accuracy commensurate with the fact that the proposed enrichment functions do not capture the stress singularity at the crack tip.
By design, DE-FEM inherits the approximation properties of IGFEM because of the construction of the finite-dimensional approximation space (see Equation (20)). DE-FEM can therefore seamlessly model weak and strong discontinuities, a desirable feature when a crack develops at the interface between two different materials. Apart from the improved kinematics, a notable feature of DE-FEM is the placement of nodes with enriched DOFs along a discontinuity. This clearly gives DE-FEM an edge over X/GFEM, as the complexity of the latter in terms of computer implementation has hindered its wide adoption in commercial software packages. 29 In fact, implementing DE-FEM in a standard finite element package requires mainly a means to partition elements into integration elements-with their corresponding routine for computing local arrays-and a post-processing engine for visualizing results. Everything else, from the mesh data structures to the enforcement of essential boundary conditions, remains roughly the same.
The most striking advantage of DE-FEM is the enforcement of essential boundary conditions. This was demonstrated in an earlier publication, 14 where the IGFEM formulation (weak enrichments only) coupled to Lagrange multipliers solved the fundamental issue of artificial oscillations in the traction field while imposing Dirichlet boundary conditions on an external boundary that does not match the mesh. Here, we have demonstrated that the same feature holds when constraining enriched DOFs along a line internal to the computational domain using penalty augmentation. It bears emphasis that constraining extra DOFs within the domain in X/GFEM using the simple penalty augmentation approach yields unreliable results. 16 This is clearly a shortcoming of the classical X/GFEM formulation. Even though headway has been made for imposing Dirichlet boundary conditions on non-matching meshes, 15 Mesh-discontinuity interaction The mesh independence of DE-FEM does not come for free. The complexity of the creation of a matching mesh in standard FEM is transferred to the mesh-discontinuity interaction and the corresponding enriched formulation. A robust geometric engine is needed for the modifications to the mesh data structure when discontinuities are present in the computational domain. The pseudo-code for the mesh-discontinuity interaction procedure is listed in Algorithm 1.
The geometric engine first detects the finite elements that are crossed by the discontinuities. Once an element has been tagged, new nodes are created at the element edge-discontinuity intersections. An approach appealing for computer implementation when dealing with strong discontinuities is to create double nodes, one for two weak DOFs and one for two strong DOFs; this simplifies the post-processing step for visualization as each node is correctly connected to integration elements at either side of a crack. If an element is only partially crossed by a crack, only one node is created within the element domain to hold the weak DOFs; since there is no displacement jump at the tip, there is no need to add a node for strong DOFs.
Noteworthy, the same data structures used in standard FEM can be used with DE-FEM with no modification for storing information related to enriched nodes; new nodes are added to the original node array, and enriched DOFs are stored in the same solution vector as standard DOFs. Assuming zero-indexed arrays, any DOF can be accessed in constant time from the solution vector by multiplying the node id (either standard or enriched) by the number of DOFs per node and adding the desired DOF.
Nodes are connected by cohesive elements if the problem deals with cohesive cracks. In addition, the element sub-domains at both sides of the discontinuity are divided into integration elements that will later be used in the integration and assembly procedures. Elements that are partially crossed are also subdivided into integration elements that are created with the inner tip node as one of the vertices. All of this implies that (i) the parent element is no longer used in the assembly process and (ii) the mesh data structure has to be modified accordingly to include the new nodes and integration elements.
The geometric engine is also responsible for storing the weights for weak and strong discontinuities at the newly created nodes (this information will later be used to compute the enrichment functions). As explained earlier in Section 2.2, determining the weights associated with strongly enriched DOFs relies on the discontinuities having a given orientation. As a result, integration elements would lay on either the positive or negative side of the crack, although this distinction is arbitrary and it is only important to keep consistency. In our implementation, given a point p, we find its closest point along the crack y, and we compute the sign (n c · (p − y)) to determine the region to which point p belongs.
Given the mesh element set  and assuming || ≪ || with |·| denoting set cardinality, the time complexity for the pseudo-code given in Algorithm 1 is  (||). An efficient implementation would make use of a space decomposition technique to speedup discontinuity-element intersections. Choices of space decomposition techniques include, among others, tree data structures (binary space partition (BSP) trees, k-d trees, octrees) and grids. 31 Integration All elements that do not interact with discontinuities are treated as in standard FEM. For integration elements, the computation of the local stiffness matrix and force vector is outlined in Algorithm 2.
As shown in the pseudo-code, integration is performed on the integration element following an isoparametric approach where the master element coordinate is denoted . At each Gauss quadrature point, the algorithm computes the shapes functions for the parent element and the enrichment functions, together with their derivatives. Therefore, the code given is split in a part that deals with the integration element and a part that deals with the parent mesh element. Once the functions and their derivatives are stacked into arrays, the contribution to the stiffness matrix and the force vector follows standard procedures.
It was mentioned earlier that one of the advantages of the formulation is to have a completely decoupled assembly of cohesive elements. For completeness, Algorithm 3 provides the assembly routine for cohesive elements.
It can be seen that the algorithm needs only information about the strong DOFs i and j associated with the end points of a cohesive element. It is shown that the orientation of the cohesive matrix is taken into account by using a rotation matrix and that a residual vector is also computed in case the cohesive law used is nonlinear.
Boundary conditions Imposing boundary conditions is as simple as in standard FEM. If a discontinuity intersects a Dirichlet boundary, as demonstrated in Examples 3.1 and 3.2, the corresponding DOFs can be readily prescribed after evaluating Equation (20).
In Section 3.1, we demonstrated a case where a natural boundary condition is prescribed at one of the edges that intersects the crack. We showed how a jump in the BC values can be readily taken into consideration by splitting the boundary integral into two parts at either side of the crack. The pseudo-code for this case is omitted because it would be similar to that given in Algorithm 2, with the main differences being (i) integration is performed over 1-D elements; (ii) there is no need to compute derivatives; and (iii) only a contribution to the force vector is needed.
Post-processing As mentioned earlier, nodes of the original mesh conserve the Kronecker property, and thus post-processing is only needed if a better representation of the displacement field is required, e.g., for visualizing the crack opening. Note that even the displacement jump at nodes along the discontinuity can be obtained without requiring a post-processing stage.
The pseudo-code for the post-processing procedure is listed in Algorithm 4. The main loop over integration elements is performed to change the connectivity table of those elements that reside on one side of the crack. The displacement field is computed at both sides of the crack through Equation (20), and the global displacement vector U is modified accordingly. The data structures in the algorithm with a ± subscript imply two sets (both sides of the crack). Implicit in the algorithm is the fact that enriched nodes are visited once. Since nodes of the original mesh keep their values, no further modification is required there.

APPENDIX B: ELEMENT ARRAYS FOR A ONE-DIMENSIONAL X/GFEM ELEMENT WITH DISCONTINUOUS ENRICHMENT
We report the expression of the relation KU = F for the 1-D bar of Figure 1 using a two-node X/GFEM element enriched with the Heaviside step function Equation (14). This enriched element has four DOFs, two per node, labeled withû (standard) and u (enriched). With reference to the displacement contribution to the discretized weak formulation of the governing equations described by Simone et al., 32  ] . (B1) The entries in the stiffness matrix are computed according to where N is a matrix containing one-dimensional linear shape functions, B is a matrix containing derivatives of these shape functions, and the surface integral has been computed by evaluating the integrand at x = x Γ , the discontinuity location. The final expression of the stiffness is given in Equation (15). The boundary condition at the clamped end (u(0) = 0) poses a constrain only on the standard DOF û 1 (û 1 = 0) as it can be easily determined by expanding the X/GFEM approximation in Simone et al. 32 , Equation (4).
As for the natural boundary conditions, they are directly enforced in F through fû = ∫ Γ t N ⊺t dΓ and (B3a) where Γ + t is the part of the boundary where prescribed tractionst are applied with the enrichment function different from zero. With reference to the 1-D problem in Section 2.1, there is an unknown force R at the clamped end that in principle should be added to both DOFs active at x = 0. However, since the enrichment function is zero at that location, only the entry corresponding to û 1 is filled with R, see Equation (B3a). Analogous considerations hold for the concentrated force of intensity P at x = L: since the enrichment function is equal to one at x = L, both DOF entries are filled with P. The right-hand side of the full system is given in Equation (16).

APPENDIX C: ASYMPTOTIC FIELDS
For completeness, we provide the expressions of the displacement and stress fields in the vicinity of a crack due to a remote tensile stress: In Equation (C1a) and (C1b), r and are the polar coordinates from the crack tip, G is the shear modulus of the material, and is Kosolov's constant given by The stress (in Voigt's notation) is given by

APPENDIX D: COMPUTATION OF WEIGHTS FOR A TRIANGULAR ELEMENT
In this section, we list the Mathematica ® code that computes the weights for the triangle in Figure 5. The code creates an interactive manipulation environment where the location of the triangle vertices and the points that define the discontinuity can be moved, and results are computed in real time. The code is not optimized for speed and does not constrain the line-to-line intersections. Yet it can be used to understand through dynamic manipulation how the weights of the discontinuous enrichment functions are computed. In