A stable discontinuity-enriched finite element method for 3-D problems containing weak and strong discontinuities

A new enriched finite element technique, named the Discontinuity-Enriched Finite Element Method (DE-FEM), was introduced recently for solving problems with both weak and strong discontinuities in 2-D. In this mesh-independent procedure, enriched degrees of freedom are added to new nodes collocated at the intersections between discontinuities and the sides of finite elements of the background mesh. In this work we extend DE-FEM to 3-D and describe in detail the implementation of a geometric engine capable of handling interactions between discontinuities and the background mesh. Several numerical examples in linear elastic fracture mechanics demonstrate the capability and performance of DE-FEM in handling discontinuities in a fully mesh-independent manner. We compare convergence properties and the ability to extract stress intensity factors with standard FEM. Most importantly, we show DE-FEM provides a stable formulation with regard to the condition number of the resulting system stiffness matrix.


Introduction
Enriched finite element methods have fundamentally changed the modeling of problems containing discontinuities. By means of enrichment functions that incorporate the required field jumps, these methods are able to completely decouple the geometrical description of discontinuities from the underlying finite element (FE) discretization, eliminating the need for creating matching or discontinuity-conforming meshes. The Discontinuity-Enriched Finite Element Method (DE-FEM), recently proposed by Aragón and Simone for 2-D problems [1], uses a single formulation to model problems containing both weak and strong discontinuities (C −1 -continuous field) by placing enriched degrees of freedom (DOFs) only to nodes created along discontinuities. This new versatile procedure, which has several advantages over other commonly used enriched formulations, is demonstrated in this manuscript for 3-D problems in elastostatics.
Standard FEM can be used to model problems with discontinuities, but it requires a matching discretization where the sides of elements align to material interfaces or cracks. When dealing with fracture problems, mesh adaptivity in the neighborhood of crack tips is usually needed to improve the FE approximation at the expense of a more complex mesh creation. Alternatively, the stress singularity can be captured by using special elements arranged concentrically at the crack tip, e.g., quadrilateral finite elements where two nodes in parametric space map to the crack tip location in physical space [2], or quarter-point finite elements, where the mid-size nodes of quadratic Lagrange finite elements are moved closer to the crack tip [3]. Although these special meshes can capture the stress singularity, generating them is a time-consuming process. Creating matching meshes in general can be quite challenging depending on the morphology of the problem, particularly because of the strict conditions on mesh quality that are required-elements with bad aspect ratios could reduce the approximation accuracy [4]; moreover, robustness is still questionable, especially for 3-D FE meshes, as issues remain in generating meshes that correctly match boundaries [5]. Recently, numerical techniques such as Universal Meshes [6] and the Conforming to Interface Structured Adaptive Mesh Refinement (CISAMR) [7] have been proposed to modify meshes locally to discontinuities. These methods can therefore guarantee adequate discretizations for the analysis by ensuring elements with proper aspect ratios-which are crucial for recovering accurate gradient fields-while avoiding the need for global remeshing.
The eXtended/Generalized Finite Element Method (X/GFEM) [8][9][10] solves the aforementioned shortcomings in an elegant manner by enriching the approximation space with functions that can reproduce the discontinuities [11,12]. In X/GFEM, the discretization is then decoupled from discontinuities by detecting cut elements, and by adding enriched DOFs to the standard mesh nodes of cut elements so that discontinuities can be resolved by the formulation. Although X/GFEM provides great flexibility in the choice of underlying FE discretizations, many properties that are usually taken for granted in standard FEM are lost, including the physical meaning of DOFs (unless shifting is used [13]), the straightforward imposition of essential (Dirichlet) boundary conditions [14,15], and most importantly, the stability-here understood as the condition number of the resulting system matrix, which in standard FEM scales with mesh size h as O ( h −2 ) . Indeed, several enrichment functions could result in an unstable formulation where the condition number of the system matrix grows much faster than that of standard FEM. Pursuing Stable GFEM (SGFEM) that is devoid of this issue has been the focus of recent research in the field [16], with SGFEM proposed for weak discontinuities in [17,18] and for strong discontinuities in [19][20][21]. Moreover, depending on the type of enrichment functions used, special techniques are needed in X/GFEM to prevent the loss of accuracy in blending elements, i.e., non-cut elements that share enriched DOFs [22][23][24]. Finally, the computer implementation of X/GFEM is also more involved than that of standard FEM [25] since a variable number of DOFs per mesh node is required, and quadrature rules used depend on the type of enrichment functions-for fracture problems, for instance, elements with singular enrichments require special integration techniques [26].
Within the realm of enriched formulations, the Discontinuity-Enriched Finite Element Method presents a new paradigm for the mesh-independent analysis of weak and strong discontinuities [1]. DE-FEM can be seen as a technique that combines X/GFEM's most salient mesh-independent property-by decoupling the FE discretization from the problem's complex geometric features-while keeping the attractive properties of standard FEM. Indeed, in DE-FEM enrichment functions vanish at mesh nodes by construction, and as a consequence DOFs associated to the latter keep their physical interpretation, e.g., displacement at the node location. In addition, essential boundary conditions can be prescribed directly as in standard FEM-both standard and enriched DOFs [27]. The computer implementation is also straightforward since the new enriched nodes are added to the same data structure used to store the original mesh nodes, and their corresponding enriched DOFs are retrieved from the solution vector in the same way as standard DOFs. Because DE-FEM was designed for simplicity in the formulation (and thus in the implementation), enrichment functions that capture the stress singularities along crack fronts are not used (although they may also be added). As a result, convergence rates for singular problems are not optimal.
In the absence of strong discontinuities, DE-FEM simplifies to the Interface-enriched Generalized Finite Element Method (IGFEM) [28] or to the Hierarchical Interface-enriched Finite Element Method (HIFEM) [29], which were devised for weak discontinuities alone [30][31][32][33][34][35][36]. HIFEM builds on IGFEM to resolve multiple weak discontinuities within a single element-and even intersecting discontinuities-by means of a hierarchical construction of enrichment functions. Given the comprehensive literature on the study of weak discontinuities with IGFEM/HIFEM, here we devote ourselves to problems containing strong discontinuities and to their interactions with weak discontinuities. DE-FEM was first proposed for analyzing 2-D problems in elastostatics. In this paper we extend DE-FEM to 3-D, and we discuss in detail its computer implementation in a displacement-based finite element code. In particular, we discuss algorithmic considerations for a 3-D geometric engine that handles the interactions between the background (original) mesh and material interfaces and/or cracks. By means of a discontinuous patch test, we demonstrate DE-FEM can generate two independent kinematic fields as long as proper care is taken in the construction of a conforming integration mesh. A convergence study then demonstrates the accuracy and convergence rates of the method, which is then used to extract stress intensity factors (SIFs). Most importantly, we show that DE-FEM is a stable formulation, i.e., the condition number of the stiffness matrix grows at roughly the same rate as that of standard FEM. Finally, we showcase DE-FEM in the challenging problem of intersecting weak and strong discontinuities, which is handled effortlessly by means of the hierarchical implementationà la HIFEM.

Formulation
Consider the elastostatics boundary value problem for a cracked body ⊂ R 3 with closure and boundary ∂ ≡ = \ . The body is subjected to a non-homogeneous essential boundary condition (BC)ū : u → R 3 and to a traction BCt : t → R 3 , the latter having a zero value along a crack c ⊂ t , as shown in Fig. 1. The abstract form of the weak (variational) formulation is: where L (·) and B (·) are the linear and bilinear forms, respectively,ũ is a vector-valued function that satisfies u| u =ū, and V 0 is the vector-valued function space with components in H 1 0 ( ) (Sobolev space that satisfies homogeneous essential boundary conditions on u ). For elastostatics, the linear and bilinear forms are given, respectively, by and In (2), b : → R 3 is the body force, and in (3) the stress tensor σ : → R 3 × R 3 obeys Hooke's law on the linearized strain, i.e., σ = Cϵ, with ϵ(v) = 1 2 ( ∇v + ∇v ⊤ ) .
In order to solve the finite-dimensional form of (1), the domain is discretized into finite elements e i (simplexes in particular) such that ∪ i e i = h ≈ . Following a Bubnov-Galerkin approach, both the trial solution and the weight function are chosen from the DE-FEM space: where the standard FE space (first term) is augmented by an enriched space with special enrichment functions that can capture the jump in the gradient (weak term) and the primal (strong term) fields. In the standard FEM term, ι h is the index set containing all original mesh nodes; N i denotes the ith Lagrange shape function and U i are the corresponding standard DOFs. In the enriched term, ι w and ι s are index sets of enriched nodes created along discontinuities. For weak discontinuities, these nodes are associated with enrichment functions ψ i and corresponding enriched DOFs α i . In addition to these, strong discontinuities add a new set of nodes associated with enrichment functions χ i and enriched DOFs β i . In order to have an arbitrary number of discontinuities interact with a single element, we follow the work of Soghrati [29]. The procedure is schematically shown in Fig. 2, where an element e i is split by two discontinuities. The first discontinuity creates a first level of integration subdomains, and these in turn become the parent elements for a second discontinuity. These relationships are stored in an ordered tree data structure, which is used during numerical quadrature to build the enrichment functions in a hierarchical manner. Therefore, for the eth cut element, the approximation u h ∈ S h e can be written as where h = {1, . . . , D} is the index set used to represent D discontinuities (both weak and strong) that interact with the element, and h s ⊆ h only includes the indices for strong discontinuities. The formulation in uncut elements remains the same as in standard FEM, i.e., only the first term in (5) is used. The hierarchical information is stored in an ordered tree that is also used for constructing enrichment functions during the numerical quadrature of local element arrays.
For the sake of simplicity, we denote by ϕ ≡ [ a vector that stacks shape and enrichment functions that are non-zero on the eth element and by B ≡ ∂ϕ = [ ∂ N ∂ψ ∂χ ] the strain-displacement matrix, which is obtained by applying to ϕ the differential operator ∂: The element local stiffness matrix and the nodal load vector are given, respectively, by where C is the constitutive matrix given by λ and µ are the Lamé constants, and 1 and I are the unity and identity matrices, respectively. Details about the construction of k e and f e , together with pseudo-code for constructing the enrichment functions in a hierarchical manner, are given later in Section 4. The finite-dimensional form obtained by taking into account all finite elements leads to an augmented system of linear equations In (10) we made explicit the partition of global arrays into standard and enriched parts. Aragón and Simone [1] discuss the formulation in more detail for 2-D elastostatics, considering also the case of cohesive tractions acting on cracks.

Geometric engine
The complexity of mesh generation in standard FEM is transferred in DE-FEM to a geometric engine that takes care of the interaction between a background mesh (usually structured) and the discontinuities. The engine needs to be designed for handling complex scenarios with multiple discontinuities, as shown for instance in Fig. 3, and to handle not only cracks (strong discontinuities) but also material interfaces (weak discontinuities).
The operations of the geometric engine are summarized in the flow chart of Fig. 4. The engine takes as input the background mesh and the set of discontinuities. In this work, FE models are discretized using linear tetrahedra, cracks are described explicitly through planar regions (analogous to line segments in 2-D), and material interfaces are described implicitly through level set functions.
The engine first loops over discontinuities and then loops over mesh elements for each discontinuity. Tetrahedral elements cut by the first discontinuity are detected. The spatial locations of new enriched nodes is then determined by performing intersection tests between the discontinuity and element sides. At these locations we create the corresponding number of enriched nodes depending on the type of discontinuity; one enriched node per intersection (associated with weak DOFs) is generated if the discontinuity is weak. If the discontinuity is strong, two enriched nodes are created per intersection (each associated with either weak or strong DOFs) unless the node is at the crack   front (associated with weak enriched DOFs only). These enriched nodes, together with those of the original mesh, are later used to create tetrahedral integration elements.
It is worth noting that the engine described here would not be efficient when dealing with a great number of discontinuities. Conducting intersection tests for every element-discontinuity pair leads to the so called allpairs weakness, which results in an algorithm of complexity O (D E) for E mesh elements and D discontinuities. The efficiency of the interaction can be greatly improved by performing heuristics conducted on bounding boxes. Furthermore, space partitioning techniques can speedup finding elements intersected by discontinuities. Among these one could use octrees, k-d trees, R * -trees, and grids [37]. Finally, these techniques could be supplemented with fast marching methods and graph traversal algorithms to look only for neighboring elements of already cut elements.
A crucial task of the geometric engine is to determine the location of crack fronts. Fig. 5 shows that a planar region, which represents a 3-D crack, intersects with a tetrahedral element partially. In order to determine the location of nodes along the crack front, intersections tests are performed between the crack and all faces of the tetrahedron. Numerical precision is of major concern when dealing with floating point arithmetic, especially for determining whether a point lies on an edge or on a face. In this work we adopt a tolerance that depends on the mesh size h for finding intersections.
Other operations performed by the geometric engine include creating groups of nodes and elements that are associated with discontinuities (the latter are used for assigning the correct material properties), generating integration elements, and updating the original mesh data structure.

Creation of integration elements
Although it is not strictly necessary, creating an integration tetrahedral mesh is convenient for four reasons: (i) Lagrange shape functions in integration tetrahedra are used to construct the enrichment functions. Storing integration elements in an ordered tree data structure further facilitates the construction of enrichment functions at each level of hierarchy-the resulting recursive quadrature algorithm is general and its computer implementation is straightforward; (ii) By creating integration elements at either side of the discontinuities we ensure the enrichment functions are smooth, and thus they can be integrated with the minimum number of integration points (parent  shape and enrichment functions are linear, and thus their corresponding derivatives are constant); (iii) Tetrahedral elements are the simplest 3-D elements to split, so they ease the handling of multiple discontinuities within a single parent mesh element; and (iv) Integration elements help in the post-processing stage so that results can be visualized correctly. It is worth noting that the creation of this integration mesh is a procedure bound locally to discontinuities, and as such it is completely different than performing remeshing.
Efficient algorithms for splitting a tetrahedron is one of our main concerns. There are four distinct configurations that result from completely splitting a tetrahedron with a planar region [34], as shown in Figs. 6(a)-6(d): (a) two tetrahedra; (b) a tetrahedron and a pyramid; (c) a tetrahedron and a prism; and (d) two prisms. A pyramid can be considered to be composed by two tetrahedra (Fig. 6(e)). Therefore, we can get three new tetrahedral elements in case (b). Normally, tetrahedralizing a prism ( Fig. 6(f)) can produce three tetrahedra under the premise that any two diagonals of three quadrilateral sides share the same vertex. However, in the case that the diagonals of the body are already prescribed (for example by the fact that neighboring elements have been tetrahedralized), it might be impossible to create new tetrahedra directly. This situation is alleviated by adding an extra vertex (called Steiner point) inside the prism (see Fig. 6(g)). Then, this prism can be split into eight tetrahedra that share this vertex.
While fully splitting a tetrahedron into several tetrahedra can be accomplished by using basic splitting rules [38], the complexity of the geometric operations increases rapidly when dealing with partial splits (see Fig. 7) or when a crack resides completely inside a tetrahedron. The latter situation could be alleviated by refining the mesh, but this in turn increases the computational cost. Hence, in this work we use constrained Delaunay tetrahedralization for creating new tetrahedral sub-domains [39].
Since each background element intersected by discontinuities is tetrahedralized in sequence, the constrained Delaunay algorithm is needed to avoid the creation of a non-conforming integration mesh, which as explained later has important consequences. Fig. 8 shows such a case, where two contiguous tetrahedral elements split by Mesh tetrahedral element e i intersected first by c1 (which creates the first level of integration elements) and then by c2 (which creates the second level). The generated children elements are stored in an ordered tree hierarchical data structure (right). a discontinuity could have sub-tetrahedra with non-matching faces. In other words, there exists the possibility of creating integration tetrahedra in element e m with triangular faces (1, 2, 4) and (2,3,4), whereas those tetrahedra in element e n could have faces (1, 2, 3) and (1,3,4). In order to avoid this situation, the constrained Delaunay function sets the faces of already created integration elements as the constraints for new elements sharing the same face. It should be emphasized that we are not conducting remeshing since the partition of unity (PoU) of the background mesh is kept intact, and constrained Delaunay is only used to locally generate integration elements; because these are used for numerical quadrature and for constructing enrichment functions, their requirements in terms of aspect ratios are less strict than those needed for Delaunay meshes.
The discussion above deals only with volumetric elements. It should be noted, however, that lower-dimensional elements (3-node triangles and 2-node lines) are also needed at times, e.g., to apply boundary conditions-consider a distributed pressure applied to triangular faces of tetrahedra split by a crack. Therefore, lower-dimensional elements are also extracted from the created volumetric integration elements.

Hierarchical bookkeeping
An ordered tree data structure is used to store the relationships between original (parent) and new (children) elements. This hierarchical structure later facilitates the evaluation of weak and strong enrichment functions during the numerical quadrature of local element arrays. Children integration elements can in turn become parents of new children elements deeper in the hierarchy created by subsequent discontinuities.
The procedure is illustrated in Fig. 9 for a tetrahedral element e i cut by two discontinuities c1 and c2 . This mesh element becomes the parent element for integration elements generated at the first hierarchical level by the first discontinuity c1 . These children elements in turn become the parents for newly created integration elements at the second hierarchical level by a second discontinuity c2 . The final hierarchical structure is also illustrated in Fig. 9.
This procedure can be used to create integration elements for mesh elements split by an arbitrary number of discontinuities. An element split by n discontinuities has n + 1 hierarchical levels in the tree, and this tree has the depth of the element split with the most number of discontinuities.

Enrichment functions
In order to calculate the element stiffness matrix k e and force vector f e in (7), we need to evaluate shape and enrichment functions, together with their derivatives. Lagrange shape functions N are obtained as in standard FEM, following an iso-parametric formulation. Here we focus on the evaluation of enrichment functions ψ and χ .

Single discontinuity
Consider a tetrahedral element e i intersected by a single crack c , as shown in Fig. 10. Because the element is split by a single discontinuity, all children elements created will be leaves (integration elements) in the ordered tree. The orientation of the crack c , given by its normal n, is used to define positive and negative regions. Three intersection points, labeled x i , x j , and x k , are determined by the geometric engine. Enriched nodes are then created at these locations and denoted with the same symbols henceforth for simplicity. 1 A positive (negative) sign superscript is used if the node is referred in the positive (negative) side of the crack. In addition, c subdivides e i into four children (integration) elements such that e (1) 1 lays on at the positive side, and e (1) 2 , e (1) 3 , and e (1) 4 are on the negative side. Enrichment functions will be constructed with the aid of standard Lagrange shape functions in these integration elements. Now consider the leaf integration element e (1) 1 on the positive side, which has three enriched nodes x i + , x j + , and x k + . The Lagrange shape functions in this element, i.e., N (1) , are used to create the weak and strong enrichment functions acting on this element. These functions can be expressed as and where w z and s z are weak and strong scaling (weight) factors for the zth enriched node, respectively, yet to be determined.
Similarly, for a children element on the negative side, such as e (1) 2 , weak and strong enrichment functions acting on the element are given by and where w z and s z are the zth enriched node's scaling factors for the negative side, and N (1) are the Lagrange shape functions of this integration element. In order to enforce C 0 -continuity of the weak enrichment function, the weak scaling factors in both sides are given by the same expression (i.e., w z + = w z − , z = {i, j, k}). For the weak scaling, we use a theoretically derived optimal function in order to preserve stability in the formulation [40]: where 0 ≤ ζ z ≤ 1 is the relative length between the enriched and original nodes measured along the corresponding element edge. For enriched node are the Euclidean distances between x i and x 4 , x 1 measured along the edge of the parent tetrahedron (see Fig. 10). The strong scaling factors are responsible for the field jump. Similarly to DE-FEM in 2-D [1], we choose the strong scaling factors so that the field jump has a unit magnitude, As a result, the associated strong DOFs β z physically represent the displacement jump at x z . For example, the strong scalings corresponding to the enriched node x i are obtained by setting s i + = −ζ i and s i − = 1 − ζ i . The weak and strong scaling factors for enriched nodes x j and x k are obtained analogously.
For enriched node x i , which is connected to integration elements e (1) 1 and e (1) 2 , the enrichment functions ψ i and χ i are: and respectively. Fig. 11 shows a schematic representation of weak and strong enrichment functions, which attain a maximum absolute value at the location of the enriched node and decrease linearly to 0 at other element nodes (both standard and enriched).
It is worth mentioning that in DE-FEM the condition [[χ z ]] = χ | x z + − χ | x z − = 1, z = {i, j, k} is fulfilled by using the partition of unity property of the background mesh element [1]. Specifically, for enriched node x i , (16) is equivalent to using one parent nodal shape function on one side of the discontinuity and the negative of another parent nodal shape function on the other side. To wit, where the superscript (0) is used to indicate that the function belongs to the parent mesh element.

Multiple discontinuities
We now consider the more general situation where a mesh element is cut by more than one discontinuity. For constructing enrichment functions in a hierarchical way, we traverse the tree starting from leaf elements.    Consider the case shown in Fig. 12, which shows a mesh element crossed by two non-intersected discontinuities.
The figure shows the hierarchical relationship between e i (mesh element), e (1) 4 (at hierarchical level 1), and e (2) 5 , e (2) 6 , and e (2) 7 (at hierarchical level 2). We take enriched node x q created after splitting children element e (1) 4 by a second discontinuity. This enriched node is associated with leaf integration elements e (2) 5 , e (2) 6 and e (2) 7 at the second hierarchical level. We follow the same strategy as that used for handling one discontinuity to construct the enrichment functions associated with x q . To wit, and where w q , s q ± are, as before, the weight factors, and N (2) ι (x), ι = { q + , q − } the Lagrange shape functions. It should be mentioned that an enrichment function at the second level of hierarchy-deepest level in this example-is nonzero only in elements connected to its corresponding enriched node (the enriched node's support). The function is exactly zero in integration elements created at shallower hierarchical levels outside the enriched node's support.
The same strategy is used to obtain enrichment functions for any enriched node at any hierarchical level.
Integration element e (1) 4 , that is created by the first discontinuity, is connected to mesh nodes x 1 , x 2 , x 3 and to enriched node x k . For the latter, the weak and strong enrichment functions ψ k and χ k are obtained, respectively, by and After evaluating enrichment functions (and their derivatives), we append them to the Lagrange shape functions of parent elements to calculate the vector ϕ (and the derivatives are appended to those of parent shape functions to get the strain-displacement matrix B). For instance, for one of the integration elements in the second level of hierarchy, such as e (2) 5 , the arrays ϕ and B are given by ϕ = [ where N (0) ι (x), ι = {1, 2, 3, 4} are the Lagrange shape functions of the parent mesh element e i . Elements not intersected by discontinuities are handled as in the standard FEM. On cut elements, the hierarchical structure adopted for handling an arbitrary number of intersections eases the computer implementation for conducting numerical quadrature. As an alternative to recursion, Algorithms 1 and 2 outline pseudo-code using an iterative loop for obtaining local arrays by Gauss integration. The algorithm shown is the same irrespective of the number of discontinuities.
In order to construct the local arrays for a (leaf) integration element e (k) at hierarchical level k, we input the modified mesh M consisting of nodes and elements. The node set is composed by the original mesh nodes N and the newly created enriched nodes N e . Similarly, the element set E is composed by the original uncut mesh elements and the hierarchical tree data structure H, which contain not only all children elements but also cut mesh elements. The function getLocalArrays proceeds as in standard FEM, where the stiffness matrix k and the force vector f are obtained by performing Gauss quadrature. Specifically, at each quadrature point ξ i , the contributions to the local arrays are obtained by evaluating arrays of shape and enrichment functions N (and their derivatives B), weighted by the product between the quadrature weight i and the Jacobian determinant j.

Algorithm 1 Local stiffness matrix and force vector for an integration element
-update nodal force vector return k, f end function The differences with respect to standard FEM are accounted for in functions. For a given quadrature point, the idea is to traverse the hierarchy from the leaf element e (k) towards its root parent mesh element e (0) , so that all enrichment functions (together with their derivatives) are calculated and stacked in their corresponding arrays along the way. Only then the partition of unity shape functions (and corresponding derivatives) of the parent mesh element are accounted for. Enrichment functions and their derivatives are obtained by first iterating over the enriched nodes belonging to the leaf element e (k) . Depending on the type of discontinuity, contributions will be made to the weak term and possibly to the strong term when dealing with a strong discontinuity. Lagrange shape functions N (k) and their derivatives N (k) ,ξ are evaluated at ξ (k) i . Shape functions are used to calculate the global coordinates x and their derivatives to compute the Jacobian matrix J (k) (which is used throughout to compute derivatives with respect to global coordinate x). Notice that we save the determinant j of the Jacobian matrix of the leaf integration element only. In the function enrich, N (k) and N (k) ,ξ , together with the scaling factors w j and s j obtained from the map W, are used to construct the enrichment functions ψ (k) , χ (k) and their derivatives ψ (k) ,x , χ (k) ,x ; these are then stacked to F and F ,x arrays, respectively. Then, we obtain the parent element (either mesh element or another integration element) from the ordered tree H and the master coordinate at that element is obtained by an inverse mapping. The iteration continues until we reach the original mesh element (condition fulfilled by k = 0 in the algorithm).
Once all enrichment functions and their derivatives are computed, we add the contribution of the parent element through N (0) and their derivatives N (0) ,ξ , which are also stacked in their corresponding arrays. To account for multiple DOFs per node in elasticity, the arrays are expanded and stored in N and B; in the end these are returned together with the saved Jacobian determinant j.
]} -stack functions and derivatives ← H(k, e (k) ) -get parent element in the hierarchy ξ

Numerical examples
In this section several numerical examples are provided to illustrate the performance of DE-FEM in 3-D. Unless explicitly stated, geometrical parameters, Young's moduli, and magnitudes of applied tractions are understood under any consistent unit system.

Discontinuous patch test
Following [1], we first verify the DE-FEM formulation through a discontinuous patch test. Fig. 13(a) shows a schematic, where a cube is split by a horizontal crack c (the red plane) and tractionst = t e 1 and 2t are applied to the front surface above and below the crack, respectively. The back surface is fixed and the top and bottom ones are constrained along the e 3 direction as shown. An elastic modulus E = 10 and Poisson's ratio ν = 0 are used. The original finite element discretization comprises 16 nodes and 42 linear tetrahedra (see Fig. 13(b)); the new discretization after introducing the crack is also shown in Fig. 13(c).
The resulting stress field is shown in Fig. 13(d) on the deformed configuration. As apparent from the figure, DE-FEM passes this discontinuous patch test since there are two different constant stress states at either side of the crack. This means that DE-FEM is capable of generating two independent kinematic configurations in 3-D.  As mentioned in Section 3, a non-conforming enriched mesh may be generated by the Delaunay algorithm during the creation of children integration elements. In order to display the effect of using a non-conforming mesh, consider two contiguous mesh elements, shown in Fig. 14(a) in different colors. We modified the connectivity of their corresponding children elements below the crack by hand. Fig. 14(b) shows that the stress is no longer constant. In this work we have chosen to use constrained Delaunay to avoid these situations altogether. The error with respect to the exact solution is then measured in the energy norm, i.e.,

Convergence
where u is the analytical solution, u h is the FE approximation by either DE-FEM or standard FEM, ε and ε h are their corresponding strains, and C is the elasticity tensor. Fig. 16 shows deformed configurations for modes I, II, and III, obtained by standard FEM 16(a)-16(c) and DE-FEM 16(d)-16(f). For the latter, the deformation field is accurately visualized after post-processing [1]. Convergence results are summarized in Fig. 17, where we see that convergence rates of DE-FEM are the same as those of standard FEM with the use of matching meshes. It is worth noting that convergence rates are not optimal for either method due to the lack of special techniques for handling the singularity along the crack front.
To investigate the influence of integration elements' aspect ratios, the element-wise error in the energy norm, given by  is evaluated for different locations of the crack on the non-matching mesh shown in Fig. 15(h). The crack, which moves from the middle of the background element to the position close to background mesh nodes, cuts original hexahedral units (composed of six tetrahedra) and creates integration elements (see Figs. 18(a) and 18(b)). The darkest integration elements shown in the figure are then considered for investigating the influence of aspect ratio Υ , which we define following [41] as where l avg is the average edge length, and V is the tetrahedron volume. The aspect ratios considered range from 1.35 to 894.76. It should be mentioned that the aspect ratio of an integration element is not always minimum when the crack cuts through the middle of the background mesh. Fig. 18(c) shows the squared local error in the energy norm defined by (25) for the two integration elements e I and e II shown in Figs. 18(a) and 18(b). The figure also shows the error measure normalized by the integration elements' volume. Because the local error involves integrating over element's volume, we normalize to single out the effect of aspect ratio. Indeed the solid curves show that increasing Υ also increases the local error, although the effect is more pronounced at low values of aspect ratio. The dashed curves show the combined volume/aspect ratio effect. The increase in local error caused by increasingly bad aspect ratios is counterbalanced by the element's volume, which decreases as the crack approaches the background mesh nodes (Fig. 18(b)). For the element in the corner, which has the constant aspect ratio Υ = 1.24, the local error without normalization shows the same tendency as the above two integration elements due to the decreasing element volume, and the normalized error ranges from 1 × 10 −5 to 1 × 10 −4 . Consequently, bad aspect ratios in integration elements have not much influence in the global measure defined in (24), and reducing the mesh size even further makes the contribution of the local error in integration elements even less significant.

Stress intensity factors
We now study DE-FEM's capability to extract stress intensity factors (SIFs) for the single edge-crack tension problem of Fig. 19(a), where for a crack width a = 2, we set a/w = 0.5, l/w = 1.5, and h/w = 1.75. Tractions t = ±σ 0 e 3 are applied at top and bottom surfaces. A local Cartesian coordinate system (e ′ 1 , e ′ 2 , e ′ 3 ) is located at the center of the block, where x ′ 2 is measured along the crack front. The material properties are E = 10, and ν = 0.3. The non-matching mesh used, shown in Fig. 19(b), is defined over a 39 × 39 × 39 Cartesian grid, resulting in 355 914 linear tetrahedra.
We use the domain integral method to calculate the energy release rate and the M 1 integral method to extract SIFs [42,43]. For a point along the crack front, which is parameterized by the coordinate s (see Fig. 20), the local interaction energy release rate is given by where superscripts (1) and (2) indicate the actual and auxiliary fields, respectively, and K (1,2) i , i ∈ {I, II, III} denote their corresponding SIFs. I (s) can be extracted from the total interaction energy release rateĪ (s) under the Fig. 18. (a,b) Integration elements created when crack cuts a hexahedral unit composed by six tetrahedra: (a) crack through the center; and (b) crack cutting close to background mesh nodes. The darkest integration elements are considered for investigating the influence of aspect ratio Υ in the local error. (c) Squared local error in the energy norm defined by (25) with and without normalizing by its volume, shown by solid and dashed curves, respectively. While the former curve shows the local error increases with increasing values of aspect ratio, the latter shows that the effect of high aspect ratios is counterbalanced by smaller volumes. assumption that it is constant within the integration domain V (s) of length 2η that encloses a crack front segment (see Fig. 20). I (s) andĪ (s) are given, respectively, by andĪ where q k is a continuous function with zero value outside V (s), and µ k (s ′ ) are the components of the unit vector generated from the junction of planes at point s ′ which are perpendicular to the crack front and tangential to the   (27), (28), and (29), an appropriate auxiliary equilibrium state is chosen. For instance, if mode I is set as the auxiliary field, then K (2) I = 1, K (2) II = K (2) III = 0, and the mode I SIF of the actual field is K (1) I = E I (s)/(2(1 − ν 2 )). A similar strategy is used to obtain SIFs for modes II and III. For the integration domain to compute SIFs, we set η = 12/39 and r = 16/39. Fig. 21 shows the normalized SIF K I /σ 0 √ πa as a function of the relative coordinate x ′ 2 /l along the crack front, and the corresponding values are given in Table 1. Under plane strain conditions, the analytical normalized SIF is 2.8284 [44]. The results obtained by DE-FEM show the same trend as other reference solutions; the SIF values approach the result obtained under the plane strain condition as x ′ 2 reaches the center. Moreover, in accordance with other references, DE-FEM also displays a decreasing K I as x ′ 2 reaches the boundary of the domain. For the study of Raju and Newman [45], singularity elements (where square-root terms are introduced to approximate the displacement field) were used to capture the singular stress field along the crack front. Generating the required mesh was, however, quite involved. For the study of Sukumar et al. [12], which used X/GFEM, special enrichment functions were used to capture the singular stress field. Yet, the required computer implementation is also more involved since, among other things, singular enrichments need special integration procedures [26]. There is always a trade-off between accuracy and complexity, and DE-FEM can obtain SIFs in a fully mesh-independent manner with a simpler computer implementation than that needed for X/GFEM, at the expense of some accuracy lost-1.4% below the value obtained by X/GEM at the center [12]. All these results were obtained on a coarse mesh, so for reference we also added SIFs computed by X/GFEM where the values were obtained using the mesh with local refinement along the crack front [46]. The SIF at the center, which is affected by the finite domain boundary, is higher than the one under the plane strain assumption.

Stability
In this section we study how the condition number of the stiffness matrix grows with mesh refinement for both DE-FEM and standard FEM. As mentioned in the introduction, enriched methods like X/GFEM suffer from stability issues when discontinuities get arbitrarily close to nodes of the background mesh, and discontinuityenriched methods like HIFEM/DE-FEM are not oblivious to this problem. For an integration element, the enriched components of the local stiffness matrix are inversely proportional to the element's volume (in fact to the determinant of its Jacobian matrix) [40]. Therefore, vanishingly small integration elements quickly deteriorate the overall condition number of the system matrix. For weak discontinuities alone, i.e., for IGFEM/HIFEM, an optimal function based on a 1-D analysis can be used to scale weak enrichments, resulting in a stable formulation [40]. However, it is also shown in that work that even if weak enrichments are not scaled, using a simple Jacobi-like preconditioner also results in a stable condition number. This is the approach we follow next. As it is commonly done in SGFEM, we study the condition number of a modified stiffness matrixK = D K D, where D i j = δ i j / √ K i j and thusK ii = 1. This modification is tantamount to applying a Jacobi-like preconditioner. The condition number (K ) is then obtained as the ratio between the maximum eigenvalue λ max and the minimum (non-zero) eigenvalue λ min . We compare this condition number to that of the standard FEM component (K uu ) (refer to (10)).
Two problems are used to study the stability. The first one consists of a solid cube containing a slanted crack with an orientation that favors the creation of arbitrarily shaped integration elements (see Fig. 22(a)). The planar crack's normal is n = 0.1617e 1 + 0.1819e 2 + 0.9699e 3 . For the second example, we introduce an extra spherical material interface with radius r = 0.59 that intersects the planar crack (see Fig. 22(b)). The structured background meshes considered here are the same as those used in the convergence study.   Fig. 22(a); and (b) the intersecting discontinuities shown. Fig. 22(b) show that, while the condition number of the raw system matrix (K ) is unstable, applying a simple preconditioner results in (K ), which grows roughly at the same rate as that of the standard FEM component (K uu ).
The results of both cases are reported in Fig. 23, where the condition number is shown as a function of the reciprocal of the mesh size. Our reference curve is that of the standard FEM component (K uu ), which scales . The figure shows that the condition number of the raw stiffness matrix (K ), without the use of a preconditioner, is unstable. Using the preconditioner results in the curve (K ), which has about an order of magnitude higher condition number than the standard part. However, the rate of growth is roughly the same, and therefore DE-FEM is stable if a simple preconditioner is applied.

Bone fracture
We now use DE-FEM to simulate the fracture mechanics of a femur whose geometry was obtained from [47], as shown in Fig. 24(a). The femur is subjected to a prescribed displacement fieldū = ±1 mm e 1 at its ends and a crack c that partially splits the bone is introduced.
Bone can be considered a composite material with cancellous (inner) and cortical (outer) phases with completely different material properties [48]. For simplicity, here we consider them as isotropic materials with elastic constants E = 10.4 GPa, ν = 0.12 (cancellous bone) and E = 18.6 GPa, ν = 0.3 (cortical bone) [49,50]. In addition, a level set function is introduced to describe their interface. Fig. 24(b) shows the finite element mesh used, composed of 7355 tetrahedral elements. Although the mesh represents the femur external boundary accurately, the mesh is non-conforming not only to the interface between materials, but also to the crack c .
The results of this simulation are shown in Fig. 25, where we see that new children (integration) elements are created properly at each side of the interface (Fig. 25(a)). Fig. 25(b) displays the resulting von Mises stress field, highlighting the stress concentration near the crack front and the higher stress in cortical bone. This example showcases DE-FEM's versatility not only in modeling both weak and strong discontinuities with a unified formulation, but also in resolving them even if they intersect the same element. This is possible due to the hierarchical implementation discussed in Section 2, which allows for the splitting of elements by several discontinuities seamlessly.

Summary and conclusions
In this paper we introduced DE-FEM for solving problems with both weak and strong discontinuities in 3-D. We discussed our implementation of a geometric engine to handle complex 3-D configurations. By means of constrained Delaunay, the geometric engine handles all possible scenarios that arise when splitting tetrahedral elements with  complex interfaces and/or cracks. In addition, a hierarchical data structure was discussed for resolving an arbitrary number of discontinuities within a single mesh element. This data structure is then used construct hierarchically weak and strong enrichment functions.
The integration mesh that is created in our DE-FEM implementation serves several purposes. Firstly, the tetrahedra that result from subdividing parent mesh elements are used to construct the enriched approximation. Indeed the enrichment functions are built with the aid of Lagrange shape functions of integration subdomains. Secondly, by using tetrahedra as subdomains we guarantee robustness when handling multiple discontinuities since tetrahedra are the simplest 3-D volumes to split. Thirdly, creating integration tetrahedra with sides aligned to discontinuities ensures that the discontinuous functions within these elements are smooth. Other integration schemes that do not require subdividing into integration elements, but which can otherwise recursively perform the numerical quadrature, are also possible for weak discontinuities. However, because of the non-smooth nature of the functions, the computational cost is higher because several levels of recursion are needed along discontinuities [51]. Lastly, the integration mesh is also used to visualize the post-processed displacement field.
Through a discontinuous patch test we showed that the DE-FEM formulation is able to generate two independent kinematic fields. We also showed that DE-FEM converges at the same rate as the standard FEM with matching meshes and that DE-FEM can be used to extract stress intensity factors. Finally, we simulated the challenging problem of resolving the interplay between weak and strong discontinuities within a complex-shaped femur model. In Table 2 we compare DE-FEM to both standard FEM and to X/GFEM. We see that DE-FEM also decouples the underlying mesh from discontinuities, which is X/GFEM's most attractive feature. However, contrary to X/GFEM, DE-FEM also retains the attractive properties of standard FEM. Because enrichment functions are local to cut elements by construction, there is no need to use the mesh's partition of unity shape functions to localize enrichments. Also, since enrichments are exactly zero at nodes of the background mesh, DOFs associated to these keep their physical meaning. Another consequence of enrichments vanishing at mesh nodes is that it is straightforward to impose Dirichlet boundary conditions. Non-homogeneous essential boundary conditions can be prescribed strongly also at enrichment node locations through a simple multiple point constraint. In fact, there is ongoing work on the use of DE-FEM as an immerse boundary or fictitious domain technique with strong enforcement of Dirichlet boundary conditions [27]. Moreover, the requirements for the quality of the integration mesh are less strict than those of a proper matching mesh. Regarding stability, we showed that the condition number of the whole stiffness matrix grows at roughly the same rate as that of standard FEM. Therefore DE-FEM is stable with the use of a simple preconditioner.
The simple formulation in DE-FEM leads to a simpler computer implementation when compared to X/GFEM. The current formulation cannot get optimal convergence rates for problems with singularities because the current formulation is devoid of singular enrichments. These could be added if higher accuracy is required-e.g., for extracting stress intensity factors. Nevertheless, adding singular enrichments to the formulation would complicate the computer implementation, not only because of the need for bookkeeping new enriched DOFs (and variable DOFs per mesh node), but also because singular enrichments need to be integrated accurately, for which many more quadrature points are required. Adding DOFs also leads to increasing the cost of solving the linear system of equations. Although we did not use any enrichments to capture singularities, the predictions of SIFs are not far from other more accurate values for the same mesh size h; our predicted result at the center was just 1.4% below the value obtained by X/GFEM.
The geometric engine is the most intricate part of the implementation because of the complexity of computational geometry operations involved and the need for robustness with respect to the placement of discontinuities. This is the crux not only of DE-FEM and other interface-enriched formulations such as IGFEM or HIFEM, but also for X/GFEM when dealing with discontinuities. For DE-FEM we showed the importance of creating a conforming enriched discretization for integration. Although the geometric engine described in this work suffices for the types of problems described, there is wide room for improvement and this will be the subject of subsequent work.

Acknowledgment
where µ is the shear modulus, κ and ν 1 are equal to 3 − 4ν and 0 under the plane strain, or (3 − ν)/(1 + ν) and ν under the plane stress, respectively, r and θ are polar coordinates centered at the point along the crack front. Stress fields around the crack front are given, respectively, by where ν 2 equals ν (0) under the plane strain (plane stress).