A stable interface‐enriched formulation for immersed domains with strong enforcement of essential boundary conditions

Generating matching meshes for finite element analysis is not always a convenient choice, for instance, in cases where the location of the boundary is not known a priori or when the boundary has a complex shape. In such cases, enriched finite element methods can be used to describe the geometric features independently from the mesh. The Discontinuity‐Enriched Finite Element Method (DE‐FEM) was recently proposed for solving problems with both weak and strong discontinuities within the computational domain. In this paper, we extend DE‐FEM to treat fictitious domain problems, where the mesh‐independent boundaries might either describe a discontinuity within the object, or the boundary of the object itself. These boundaries might be given by an explicit expression or an implicit level set. We demonstrate the main assets of DE‐FEM as an immersed method by means of a number of numerical examples; we show that the method is not only stable and optimally convergent but, most importantly, that essential boundary conditions can be prescribed strongly.


INTRODUCTION
Immersed boundary techniques eliminate the need for intricate meshing algorithms by decoupling the external boundary description from the discretization mesh. In this paper, we introduce a novel method for immersed problems, where enrichments to the finite element approximation are associated to locations on the immersed boundary. This facilitates the strong enforcement of essential boundary conditions on nonmatching edges, which was not possible in immersed methods until now.
The Finite Element Method (FEM) is a well-established numerical procedure for the analysis of a wide range of problems in physics and engineering. It requires a discretization mesh whose elements align with the domain's external boundaries and internal interfaces. However, for problems with complex geometries, creating a good-quality matching or geometry-conforming mesh is a demanding procedure, where admissible element aspect ratios have to be guaranteed. The burden of meshing increases even further in cases where the domain boundaries are changing in subsequent analyses, as is the case in optimization procedures. Local mesh-adaption techniques 1-4 such as Universal Meshes 5 and the Conforming to Interface Structured Adaptive Mesh Refinement (CISAMR) 6,7 have been proposed to account for changing boundaries. These methods ensure that elements with proper aspect ratios are created. Alternatively, immersed boundary methods alleviate the burden of remeshing by describing the external boundaries independently from a usually structured discretization mesh. 8,9 In doing so, part of the complexity is merely shifted from the mesh generation to the integration scheme and mesh interactions. This shift is reasonable in cases where mesh generation is cumbersome, such as for complex-shaped problem domains 10 and for moving interfaces during optimization. 11 A number of immersed boundary techniques have been proposed which deal with external boundaries independently from the mesh. In the unfitted FEM, first introduced by Barrett and Elliot 12 and later improved by Hansbo and Hansbo,8 the basis functions in the intersected elements are restricted to the respective integration domains, resulting in a doubling of the number of basis functions in intersected elements. This method was developed further into the cut FEM (CutFEM) [13][14][15] by adding ghost penalty terms to improve the condition number of the resulting matrices. CutFEM has been studied in the context of both solid 15,16 and fluid mechanics, 11 and has recently been applied to the topology optimization of flow problems. 11 The Finite Cell Method (FCM) 9,17 is a fictitious domain procedure based on the p-version of FEM (p-FEM) where the field variables extend smoothly outside the physical domain. FCM has also been used in topology optimization 18 and has been extended to handle geometries based on Non-uniform rational basis spline (NURBS). 19,20 While all of these methods undoubtedly allow greater flexibility by decoupling the discretization from the problem's geometric features, one of the core challenges remains prescribing essential (Dirichlet) boundary conditions. To the best of our knowledge, in existing immersed boundary methods, there is currently no way to strongly prescribe essential boundary conditions to a nonmatching element side. Instead, several methods have been proposed to weakly impose essential boundary conditions, such as utilizing Lagrange multipliers 13,21,22 or employing Nitsche's method. 14,23 Natural (Neumann) boundary conditions, on the other hand, do not generally pose much of an issue, as the only requirement is a means to integrate them accurately over the immersed boundary.
In parallel, the mesh-independent analysis of discontinuities that are internal to the domain has become established practice in enriched methods such as the eXtended/Generalized FEM (X/GFEM), [24][25][26][27][28][29] which originated from partition of unity methods. 30,31 They provide an elegant solution to handling strong or weak discontinuities, referring to whether the discontinuity is present in the field itself (eg, a crack) or in its gradient (eg, a material interface), respectively. When modeling discontinuities with X/GFEM, the standard finite element approximation is augmented by enrichment functions that incorporate the desired behavior otherwise missing by the use of a nonmatching mesh (in this case the jump in the field and/or its gradient). It has been shown that this method outperforms the standard FEM in situations where changing topologies are involved, such as solidification problems, 32 fluid-structure-contact interaction problems, 33 and topology optimization. 34,35 However, X/GFEM comes with its own set of challenges, such as the need for choosing appropriate enrichment functions that do not degrade accuracy 27 and the need for special formulations for prescribing interface conditions and essential boundary conditions. 36 More recently, research efforts have focused on obtaining stable formulations through stable generalized FEMs to address issues inherent to the formulation that result in ill-conditioned matrices. 37-39 X/GFEM has been studied in the context of immersed domain problems. However, as explained in the work of Cuba-Ramos et al, 40 it suffers from boundary locking as the Lagrange multiplier space overconstraints the problem.
The Interface-enriched Generalized FEM (IGFEM) 41 was introduced as a particular type of enriched FEM for the mesh-independent modeling of weak discontinuities. [42][43][44] It differs from X/GFEM in that it places the enrichments exclusively on nodes collocated along discontinuities, instead of associating them to nodes of the original mesh. This significantly simplifies implementation, as enrichment functions are straightforward to construct by using Lagrange shape functions of integration elements, ie, subdomains created for the purpose of integrating the element local stiffness matrix and force vector. IGFEM, which converges optimally with mesh refinement, 41,42 has been demonstrated successfully in the modeling of fibre-reinforced composites, 42 the multiscale damage evolution in heterogeneous adhesives, 43 and microvascular materials with active cooling. 41,42,45 Multiple interfaces crossing a single element can be resolved recursively by using a hierarchical implementation of IGFEM called the Hierarchical Interface-enriched FEM (HIFEM). 44 This makes it straightforward to analyze interfaces that are very close together, or even n-junctions, where n interfaces meet at a single point inside an element, allowing for the analysis of composite materials such as polycrystalline microstructures. IGFEM has also been used in combination with adaptive meshing, 46 with NURBS interfaces, 47 and in an optimization setting. 48,49 Cuba-Ramos et al 40 demonstrated IGFEM as an immersed boundary method, by using Lagrange multipliers to weakly impose Dirichlet boundary conditions on non-matching boundaries. More recently, Aragón and Simone 50 introduced the Discontinuity-Enriched FEM (DE-FEM) as a generalization of IGFEM to treat both weak and strong discontinuities with a unified formulation. In addition to the flexibility of dealing with both discontinuity types, DE-FEM inherits all of the virtues of IGFEM/HIFEM: • The construction of both weak and strong enrichment functions is straightforward, as they are based on the standard Lagrange shape functions of integration elements; • Because enrichment functions vanish at original mesh nodes, the Kronecker-delta property is maintained in standard nodes, allowing essential boundary conditions to be applied in the same way as in standard FEM; • With the use of a diagonal preconditioner, or a proper scaling factor for the enrichment functions, 51 the formulation used for treating weak discontinuities is stable, ie, the condition number increases at the same rate as that of standard FEM under mesh refinement; • A hierarchical implementation of DE-FEM can analyze multiple discontinuities and n-junctions within a single element. Cracks and interfaces are allowed to intersect each other as well. 55 In this paper, we place IGFEM, HIFEM, and DE-FEM in the context of immersed problems, using a strong enforcement of essential boundary conditions based on multiple point constraints. In the absence of strong discontinuities, the DE-FEM formulation simplifies to that of IGFEM/HIFEM. As external boundaries are weak discontinuities, the method for imposing Dirichlet boundary conditions holds for IGFEM, HIFEM, and DE-FEM alike. We demonstrate the versatility of DE-FEM by means of an immersed-discontinuous patch test that includes a material interface and a crack, within just one element. Convergence of immersed IGFEM is demonstrated by means of Eshelby's inclusion problem. An example with a slightly rotated background mesh is used to demonstrate the stability of the method. Lastly, it is shown that the method is capable of dealing with complex 3D geometries by means of immersing an intricate level set function within a volumetric mesh. Summarizing, the main novelties in this paper are: • The use of DE-FEM as an immersed boundary method, where a material is only assigned to the physical parts of the domain. We include cracks and material interfaces in the immersed problems as well; • To the best of the authors' knowledge, this is the first immersed method where Dirichlet boundary conditions can be applied on nonmatching boundaries in a strong manner, by means of multiple point constraints.

FORMULATION
Consider a solid body, represented by an open domain Ω ∈ ℝ d referenced by a Cartesian coordinate system spanned by {e i } i=1..d , as shown in Figure 1. The body is composed by matrix and inclusion materials, denoted as Ω m and Ω i , respectively, such that the domain closure, denoted by an overbar, is defined as Ω = Ω m ∪ Ω i and Ω m ∩ Ω i = ∅. The boundary of the domain Ω ≡ Γ = Ω ⧵ Ω, with outward unit normal n, is composed by disjoint lower-dimensional manifolds Γ u and Γ t , such that Γ = Γ u ∪ Γ t and Γ u ∩ Γ t = ∅. We assume Γ u ≠ ∅ is the region where essential (Dirichlet) boundary conditions are prescribed. Similarly, Γ t is the region with prescribed natural (Neumann) boundary conditions. The figure also shows a traction-free crack Γ c ⊆ Γ t that interacts with both materials and might also coincide with the material interface Γ i . The crack is parameterized by the curve s, which also serves to define the orientation of the crack in space so that a unique normal n c can be defined. The normal vector field is also used to identify two regions in the immediate vicinity of the crack, shown in the figure with positive (+) and negative (−) symbols. Thus, two points that coincide in space but are situated at each side of the crack are separated after deformation by a vector , which denotes the displacement jump as a function of a local coordinate s through the crack centerline.
The setting just described corresponds to that of a solid that contains both weak and strong discontinuities. The weak discontinuity results from the mismatch in material properties between the matrix and the inclusion, while the strong discontinuity arises from the crack. In this work, we restrict ourselves to traction-free cracks, but extension to cohesive or pressure-loaded cracks is straightforward. For simplicity, the graphical representation of the problem shows a single inclusion and a single crack, but the generalization to multiple inclusions and cracks is straightforward as well. We are interested in solving the linear elastostatics problem on Ω. We denote the displacement field u in the matrix (Ω m ) and inclusion (Ω i ) by u m and u i , respectively, ie, u ≡ u| Ω , = m, i. Given the prescribed displacement with interface conditions and constitutive and continuity relations = tr ( ) + 2 , In Equation (1), · denotes the divergence operator and ∶ Ω → ℝ d × ℝ d is the second-order stress tensor that follows Hooke's law for isotropic linear elastic materials, which can be fully characterized by the Lamé parameters and .
Let (Ω) ≡ [(Ω)] d be a vector-valued function space on Ω, where each vector component of v ∈  belongs to the first-order Sobolev function space H 1 (Ω). Similarly, let  0 (Ω) ⊂  be the subset that satisfies homogeneous boundary conditions on Γ u . To deal with nonhomogeneous boundary conditions, we define the linear variety  ⋆ =ũ+ 0 as a translation of  0 by the vector-valued functionũ ∶ũ i ∈ H 1 (Ω),ũ| Γ u = u, ie, every element of  ⋆ satisfies the nonhomogeneous essential boundary condition.
The weak formulation is: Find u ∈  ⋆ such that or equivalently, find u ∈  0 such that B(u, v) = L(v) − B(ũ, v), ∀v ∈  0 . The linear and bilinear forms are given, respectively, by and For solving the problem, we choose a domain Δ ⊂ ℝ d that fully encloses our original problem domain (Δ ⊇ Ω), as illustrated in Figure 2A. This hold-all domain is discretized by finite elements so that (B) integration elements are stored in a tree structure; (C) new nodes and integration elements are added to the discretization ith element and e i ∩ e j = ∅, ∀i ≠ j. An interaction between the discretization Δ h and the problem's geometric features (boundary, interfaces, and/or cracks) then takes place, resulting in new nodes, created at the intersection of the discontinuities with the edges of elements in the discretization. Similarly, intersected, or cut, elements are divided into subdomains (called integration elements), which are added to a hierarchical data structure. For instance, Figure 2A shows an element e i that is traversed by the interface Γ i and the crack Γ c . First, this background element e i is split into three integration elements (e (1) i , i = 1, 2, 3) according to the crack Γ c . Integration element e (1) 1 is in turn split into three integration elements (e (2) i , i = 1, 2, 3). The resulting hierarchy is shown in Figure 2B. The new discretization containing new nodes and the element hierarchy will be denoted henceforth as Δ h H . Despite their name, integration elements are used for more than just integration. In fact, the purpose of these integration elements is fourfold: (i) they are used for integration of the element matrices; (ii) enrichment functions are constructed as linear combinations of their standard Lagrange shape functions; (iii) they ensure that the enrichment functions are smooth and can thus be integrated with the least number of Gauss points; and (iv) they are used to ensure the field can be displayed correctly after postprocessing.
The finite-dimensional form of Equation (8) is then solved on Δ h H by choosing our trial solution u h and our weight function w h from the discontinuity-enriched finite element space In Equation (11), h represents the index set of all nodes in Δ h H from the original discretization (shown in Figure 2C with and symbols for the degrees of freedom (DOFs) outside and inside the domain, respectively). Similarly, w and s denote index sets of weak and strong enriched nodes, respectively, which are the result of the aforementioned interactions with the background mesh (shown with and symbols in Figure 2C, respectively). The first term in (11) corresponds to the standard FEM part, where N i and U i are the Lagrange shape function and DOF vector, respectively, associated with the ith mesh node. The standard FEM space is augmented by enrichment functions. The result is an enriched space that is spanned by enrichment functions that contain the kinematic description of the discontinuities. In the approximation, i ( i ) is the enrichment function that reproduces the weak (strong) discontinuity with associated enriched DOF vector i where P 1 represents the space of first-order polynomials on e.
As previously stated, DE-FEM is a particularly versatile method which allows arbitrary configurations of interfaces, domain boundaries, and cracks. Any number of interfaces and domain boundaries are allowed to come arbitrarily close to one another, and intersections of domain boundaries with cracks are allowed. As a result, an element e i ∈ Δ h H can be intersected by an arbitrary number of discontinuities. The method has this advantage by virtue of the hierarchical construction of the shape functions in elements split by multiple interfaces, as illustrated earlier in Figure 2A and in Figure 3. Here, we follow the work of Soghrati 44 for the hierarchical construction of enrichment functions. At the element level, the approximate solution u h ∈  h can be written as where h ≡ ℤ + = {1, 2, … , D} is the index set of hierarchical levels resulting from D discontinuities that interact with the element, and h s ⊆ h represents the subset associated with strong discontinuities, since strong DOFs are not present in weak discontinuities. Note that an element intersected by D discontinuities will have D levels of hierarchy in the ordered tree. Because enrichment functions are constructed with the aid of Lagrange shape functions in integration elements, these functions are nonzero only in the domain of the cut element by construction; the functions attain their maximum absolute value at the location of the enriched node and ramp linearly to zero at nodes of the cut element. The procedure for constructing hierarchical shape functions for both weak and strong enrichments is outlined in Figure 3.
The numerical quadrature of every integration element e i ∈ Δ h H is conducted hierarchically in order to obtain the local stiffness matrix k i and force vector f i as where C is the constitutive matrix and d denotes the differential operator. In the Appendix, pseudocode is given for the computer implementation, where the traversal over the hierarchy is executed in a loop. Subsequently, following standard procedures, the discrete system of linear equations KU = F is obtained, where and denotes the standard finite element assembly operator. An important issue in immersed boundary methods is the ill-conditioning of the resulting system matrices, which can hinder obtaining a solution. Immersed methods suffer from conditioning issues due to nearly-zero contributions to the system matrix, originating from intersected elements with low volume fractions belonging to the physical domain.

FIGURE 3
Hierarchical splitting of an element crossed by a material interface (Γ e i ), defined by element subdomains Ω e i and Ω e m , and a crack Γ c (shown in red). The hierarchical construction of the corresponding enrichment functions ki and ki is illustrated. Observe that all enrichment functions can be written as linear combinations of Lagrange shape functions of integration elements. Note also that there is no strong enrichment function associated with the crack tip Several strategies have been proposed to prevent ill-conditioning in immersed methods. A possible solution is to add an artificial stiffness to the fictitious part of the domain. However, adding stiffness to void regions modifies the problem at hand, and it is not straightforward to strike a balance between stability and accuracy. Another option is to modify the FE function space; one could either simply remove all basis functions that deteriorate the condition number, or scale the FE basis functions. Lastly, a preconditioner could be employed to improve the condition number of system matrices. In this work, we choose to remove all DOFs belonging to nodes in void areas from the system matrices, while setting the material properties in the void regions to zero. We also use a diagonal preconditioner to avoid ill-conditioning. However, an optimal scaling for enrichment functions, which is a function of the geometric properties of the intersection and the material properties on each side of the interface, was proposed recently and could be used instead of the preconditioner. 51

Treatment of boundary conditions
An important issue in immersed methods relates to prescribing boundary conditions on the nonmatching boundaries Γ u and Γ t , which may coincide with the interface Γ i and crack Γ c . In immersed methods, in general, it is not straightforward to apply boundary conditions on these edges, as no DOFs are explicitly related to the boundaries. Consequently, boundary conditions have to be imposed in a weak manner. In DE-FEM, on the contrary, prescribing boundary conditions is straightforward due to the fact that enrichment functions vanish at nodes of the background mesh and that their associated enriched nodes are collocated exclusively along discontinuities. In fact, in DE-FEM, it is possible to impose Dirichlet boundary conditions in a strong manner by solving local problems, which is a distinctive asset that gives DE-FEM an advantage over other immersed domain methods.

Dirichlet boundary conditions
We first start our discussion on prescribing essential boundary conditions on enriched nodes placed along the immersed boundary, which is a weak discontinuity described by weak enriched DOFs alone. By denoting x , the spatial coordinate of an enriched node associated with DOFs kj , we simply solve for the latter using (12) Note that only the enrichment functions that are nonzero at x need to be taken into account. In practice, this means that only the higher levels of hierarchy (n < k) need to be considered. Equation (15) can be used to find the value kj for all prescribed nodes, as a function of U i , ni and ni . The resulting system of equations can subsequently be written in the form of a multiple point constraint (MPC): where the vector [Ũ̃̃] ⊺ contains only the remaining nonprescribed DOFs. The matrix T contains the contributions of other nodes, both original and enriched, to the prescribed DOFs. The vector g contains the prescribed values u and .
The modified system matrix and right-hand side can then be computed as In these reduced matrices, the displacements are strongly enforced.
To illustrate the procedure, consider the triangular element e in Figure 3, split by a boundary Γ i so that part of the element lies in the domain Ω m and part lies in the domain Ω i . In this example, Ω m is void. The element is also split by a crack Γ c so that multiple integration elements are created hierarchically. We now prescribe the primary field variable u over Γ i (ie, u| Γ i = u), which is discontinuous at the crack. To obtain the DOF values that should be prescribed on the first level of hierarchy, corresponding to coordinates x 1 and x 2 in the figure, we evaluate (15) for 11 and 12 , respectively. Note that the DOFs in the void area are set to zero and are therefore removed from the system. Furthermore, enrichment functions from deeper levels of the hierarchy vanish at these coordinates, so only the shape function N 2 corresponding to U 2 will have a contribution For the next hierarchy level, only the coordinate x 4 is located on Γ i , but here, both weak ( ) and strong ( ) DOFs are present. Because physically represents the crack opening displacement, 50 their values are readily available once the jump in the displacement is known, ie, 21 Then, solving for at x 4 follows the same procedure just described for 11 and 12 . Because the prescribed field u(x ) is discontinuous at this location, we prescribe here the average displacement u avg (x 4 ). The strong DOF 21 then vanishes from the equation. These DOFs can later be interpreted as the crack opening around the average displacements. The resulting expression for 22 is written as Note that 11

Neumann boundary conditions
As stated before, tractions on the boundary Γ t follow the same elementwise assembly procedure into f as in standard FEM. In order to apply a traction t on Γ t ⊂ e, one would simply compute where N ⊺ contains not only the element's Lagrange shape functions N i , but also enrichment functions ki and ki , as As a result, f e will add contributions to DOFs that are not located on the boundary Γ t . For example, if a traction were to be applied on Γ i in Figure 3 (still assuming Ω m to be void, and therefore removing U 1 and U 3 ), contributions would appear on DOFs U 2 , 11 , 12 , 22 , and 21 , as all of their corresponding functions are nonzero on the boundary.
In the case of Neumann boundary conditions, DE-FEM also benefits from placing nodes on the interface and creating integration elements, as the boundary Γ t is explicitly known and therefore can be used for integration. Note that, as the partition of unity is not retained in cut elements, the sum of all contributions will be higher than the total applied traction. However, the obtained solutions are correct, as demonstrated in Section 3.

NUMERICAL EXAMPLES
In the following examples, a consistent unit system is assumed. Furthermore, plane-strain conditions are used for the 2-D examples, and a quadrature rule that exactly integrates integration elements is adopted, ie, one integration point is used for the linear triangular (integration) elements.

The "ultimate" discontinuity patch test
In order to establish that DE-FEM can indeed recover constant states of stress, and to demonstrate the flexibility provided by DE-FEM, we devise an immersed patch test aiming at recovering multiple independent kinematic fields; the problem contains an interface and a crack that intersect each other within an immersed domain. The problem, schematically shown in Figure 4, consists of a square domain of area L × L. The material to the left (right) of the interface has a Young's modulus E 1 = 2 (E 2 = 20), and both materials have Poisson's ratio 1 = 2 = 0, to ensure that the vertical strain 22 is zero in both materials, and a constant analytical stress is obtained above and below the crack. Regarding boundary conditions, the plate is constrained in displacement as shown in the figure and subjected to a traction per unit length ||t 1 || = 1 and ||t 2 || = 2, respectively, below and above a crack Γ c at x 2 = 0.5L.
The analytical displacement field of this patch test is given by 0.5L 0.5L

FIGURE 4
A plate composed of two materials is split in two by a crack Γ c , denoted by the red line. Tractions t 1 and t 2 are applied to the right edge, and the displacement is prescribed as shown This problem is then discretized within a single background element. This is possible because, in this case, the linear displacement field, combined with these material properties, will induce a constant state of stress. It has to be noted that, in physical problems containing multiple discontinuities, the interaction between these is likely to be more complex. A finer mesh would then be preferred to accurately capture the physics in that region. Even then, it is important to be able to describe multiple discontinuities within a single element for several reasons: • A crossing between two nonmatching discontinuities will typically occur in a single element, even after refinement.
This may occur, for example, in composite materials, crack branching (or merging), and in crack propagation in composite materials or immersed domains. • The method is more robust when discontinuities can come arbitrarily close to one another. While a mesh refinement step might be required, a linear interpolation between two interfaces might give a reasonable first estimate of the solution.
The hierarchical immersion of the problem in a single three-node triangular element is illustrated in Figure 5, together with the integration elements that result from the interaction with the mesh at each level of the hierarchy. Figure 6A shows the numerical displacement field that matches the analytical one exactly. Note that, although only the physical domain is colored, the integration elements outside the domain clearly undergo deformation as well. However, because Young's modulus in that region is exactly zero, these deformations do not lead to any stress. As shown in Figure 6B, DE-FEM is able to recover the analytical displacements and stresses exactly, despite the fact that this problem contains an interface and a crack and is immersed in a single, nonmatching, background element.

Immersed Eshelby's inclusion problem
The accuracy and convergence of the formulation are tested on the classical Eshelby's inclusion problem, as illustrated in Figure 7A. A circular inclusion of a compliant material (E 1 = 1, 1 = 0.25) and radius r i = 0.9 is immersed in a stiffer material (E 2 = 10, 2 = 0.3) with a radius of r o = 2. The circular boundary is subjected to a prescribed displacement along the outer radius 1.

2.
(B)  (23), is retrieved at either side of the cracked domain. The deformations of the integration elements that lie outside the physical domain do not induce any stress The analytical solution for the displacement field is where is a function of the Lamé parameters of the materials 1 , 2 and 1 , 2 : and the Lamé parameters are defined as = E ∕ (1+ )(1−2 ) and = E ∕ 2(1+ ) . The analytical displacement field shown in Equation (25) is nonlinear in terms of spatial coordinates and therefore cannot be recovered exactly by the linear interpolation space of the discretized model.
This circular domain is immersed in a square structured mesh as shown in Figure 7B, creating nonmatching edges for both the internal material interface Γ i as well as the external domain boundary Γ o . On the boundary, the exact solution is prescribed in the form of a Dirichlet boundary condition. Contrary to the weak enforcement of nonhomogeneous Dirichlet boundary conditions in the work of Cuba-Ramos et al, 40 here we prescribe it using multiple point constraints, as described in Section 2.1.1. The numerically obtained displacement field is shown in Figure 8. Furthermore, the tractions on the Dirichlet boundary are recovered by mapping the nodal forces (obtained as KU = F) back to d−1 elements. The recovered traction profile is smooth, as illustrated in Figure 8. Figure 9 shows an elementwise evaluation of the error fields in both the energy norm and the L 2 -norm, which we define for an element e as and where u is the analytical solution, u h is the numerical solution, (u) is the strain due to displacement field, and C the elasticity tensor. As this error evaluation is done on the level of integration elements, all elements are either completely inside or completely outside the physical domain. As we are not interested in the latter, we simply omit them from the error analysis. It is clear that the error is concentrated mostly along the material interface, as the circular interface is discretized as a piecewise linear discontinuity. Particularly, the error in the energy norm may be large in integration elements with bad aspect ratios, due to inaccurate computations of the derivatives of the field. 6,7 Note, however, that the error on the domain boundary is relatively low, as the enriched DOFs are prescribed to represent exactly the imposed displacement there.
To further investigate the influence of aspect ratio on the error, we analyze a series of Eshelby's inclusion problems on a single mesh that is matching to the material interface, but nonmatching to the outer boundary. Figure 10 shows that the error on the matching interface is of a similar order as the error in IGFEM. Changing the radius of the outer boundary while keeping the background mesh fixed changes the aspect ratio of the integration elements. The error in an integration element close to the immersed boundary, indicated in Figure 10 with a red outline, is plotted as a function of the aspect ratio in Figure 11. Both the absolute error ||u − u h || 2 (e) and the relative error ||u − u h || 2 (e) ∕||u|| 2 (e) are shown. From the former, it can be seen that the absolute error in the element decreases with aspect ratio as the volume of the elements decreases. From the relative error, it becomes clear that the aspect ratio indeed has an influence, but the error remains bounded. Therefore, elements with bad aspect ratios will not have much influence on the global error measure under mesh refinement, as is investigated next.  In order to analyze convergence under mesh refinement, we use the standard L 2 -and energy norms of the relative error, given respectively by and The convergence behavior under mesh refinement of the method is illustrated in Figure 12. Immersed DE-FEM, where both the domain boundary and the material interface are nonmatching, as described before, is compared with standard FEM, where a mesh that conforms to both circular boundaries is used. In Figure 12A, the error in both norms is plotted against the number of DOFs, whereas in Figure 12B, the error is plotted against the mesh size h. Dotted lines with slopes corresponding to optimal convergence rates are provided for reference. It is clear from the results that immersed DE-FEM has an optimal rate of convergence. Furthermore, it performs comparable with the standard FEM with respect to accuracy.
For h-convergence, DE-FEM has a slightly better accuracy, as details within a background element can be captured by DE-FEM, without a change in h. In summary, this example demonstrates that the method is optimally convergent for cases without singularities such as cracks. For the same number of DOFs, a similar error is obtained as with standard FEM, without the burden of creating a matching mesh. Furthermore, this example shows that the MPC method for prescribing Dirichlet boundary conditions in a strong manner is also able to correctly prescribe nonhomogeneous displacements on a nonmatching boundary.

Stability
Here we study the influence of refinement on the stability of the proposed methodology by looking at the condition number of the system matrices. A rectangular domain is evaluated on a background mesh that was rotated by 3 • , as illustrated in Figure 13A, to ensure that integration elements of all shapes and sizes are present. The background mesh is then refined in multiple steps, creating different irregular-shaped elements in each level of refinement. The stability is then evaluated. The problem is not subjected to any boundary conditions, but rigid body motions are accounted for by discarding the lowest three eigenvalues in the system, which are zero within numerical precision. As stated earlier, in this method, we remove all DOFs belonging to nodes in the void area of the system. The remaining stiffness matrix can in general be written as where K uu is the portion that corresponds purely to the standard FEM, K corresponds to weak enrichments, and K corresponds to strong enrichments. The off-diagonal matrices K u , K u , K contain coupling terms. Note that in this particular example, no strong discontinuities are present, and therefore the matrices K u , K and K are absent. Following the work of Kergrene et al, 39 we study the condition number of K and a modified matrixK = DKD, where D is a diagonal matrix, defined such that D ii = 1∕ √ K ii and thusK has unit values on the diagonal. The condition number is then obtained as where max and min are the maximum and minimum (nonzero) eigenvalues of either K orK. Ideally, the condition number of a DE-FEM matrix would scale under mesh refinement at the same rate as that of standard FEM, which scales as (h −2 ). In Figure 13B, the following results are shown: we compare the condition number of the full system without any preconditioner (K) and the preconditioned system (K) with the part of the matrix (K uu ) that corresponds to the standard part of the approximation. From these results, it is clear that the method indeed suffers from ill-conditioning if no measures are taken to prevent it. This happens because when an enriched node is placed close to an original node, integration elements with small areas are created. The condition number grows with 1∕j, where j is the Jacobian of these small integration elements, as was shown in the work of Aragón et al. 51 In DE-FEM without a preconditioner, the condition number is not a straight line as a function of h −1 , suggesting that the condition number is highly dependent on the intersection geometry in each particular background mesh. However, the results also show that a simple diagonal preconditioner suffices in reducing the condition number of the full matrix to that of standard FEM. Alternatively, a local scaling of the enrichment functions has been proposed in the aforementioned work 51 that results in the same reduction of the condition number. Consequently, the method is stable, and complex conditioning schemes are unnecessary in DE-FEM.

An immersed popcorn
As a final example, we demonstrate the method in 3D. To this end, we immerse a level set description of a popcorn shape 15,52-54 ( Figure 14A) into a structured cubic mesh, as illustrated in Figure 14B. This geometry is described by the function Following the work of Burman et al, 15 we choose r = 0.6, = 0.2, and A = 4. The material inside the domain is assigned a Young's modulus E = 2, a Poisson's ratio = 0, and a thermal expansion coefficient V = 0.01. A uniform temperature is applied on the surface Γ i of the popcorn, described by the zero level set. A thermal solve leads to a constant temperature throughout the domain, such that the material will try to expand uniformly. However, a homogeneous Dirichlet boundary condition is prescribed for the displacement field on the entire surface Γ i , preventing the popcorn from expanding.
The analytical result for this problem yields a constant state of stress throughout the entire physical domain with stress magnitude || || = 0.0346. In Figure 15, the numerical results are illustrated. It is clear that the constant temperature field with a value of T = 1, and the state of stress, with magnitude || || = 0.0346 are indeed retrieved.

SUMMARY AND CONCLUSIONS
In this work, we extended DE-FEM to treat immersed domain problems. The formulation, combined with a hierarchical implementation, results in a versatile mesh-independent method for solving problems containing both weak and strong discontinuities with a unified formulation. We have demonstrated, by means of a complex patch test, that a domain containing interfaces and cracks can be analyzed by enclosing it by a discretization that disregards completely all discontinuities. With the immersed Eshelby's inclusion problem, we showed that nonhomogeneous essential boundary conditions on a nonmatching boundary can be prescribed in a strong manner. Furthermore, this example was used to show that optimal convergence is attained. The stability of the method was investigated by means of a rotated mesh example, and it was found that the method remains stable under mesh refinement with the use of a diagonal preconditioner. Lastly, the extension to complex 3D structures was demonstrated. Discontinuities were represented both implicitly (via level sets as in the Eshelby and popcorn examples) and explicitly through line segments (in the immersed patch test and in the stability example).
As any other computational method, DE-FEM has advantages and disadvantages: Pros The method can truly decouple the mesh from discontinuities: different types of weak and strong discontinuities are allowed to lie arbitrarily close to one another, or even intersect within a single element. This makes the method robust against any placement of discontinuities and suitable for the analysis of complex geometries, such as fracture analysis in immersed domains. • The method is stable, ie, the condition number increases at the same rate as that of standard FEM, which scales with mesh size h as (h −2 ).

Cons
• Compared to existing immersed boundary methods, for example FCM, the geometric operations are more complex.
Enriched nodes need to be placed along the boundary, and integration elements need to be created. These geometric operations can be done efficiently but require specialized code: a geometric engine. 55 • The hierarchical construction of enrichment functions calls for a dedicated data structure, such as an ordered tree, to store the hierarchy. However, this small overhead enables the hierarchical implementation, which is a major asset. • For problems with material interfaces, it has been shown that IGFEM/HIFEM may overestimate stresses in integration elements with bad aspect ratios. 6,7 This issue, however, is not significant near Dirichlet boundaries, where essential boundary conditions are enforced exactly.
Considering DE-FEM's many interesting properties, the advantages over other immersed methods strongly depend on the application. As a final remark, for the same functionality, the complexity of the geometric operations is similar to X/GFEM. Any increased complexity in the geometric operations of DE-FEM immediately opens up a wide range of new capabilities, such as the hierarchical implementation, interface with sharp corners and crack tips inside an element, and immersed domains inside a single element. This trade-off between generality of the method and simplicity of implementation is to be expected, and in cases where this flexibility is not needed, simplified geometric operations might be executed. This makes DE-FEM suitable as an all-purpose immersed domain method and especially suitable for problems with complex configurations of discontinuities.
How to cite this article: van

APPENDIX PSEUDOCODE FOR HIERARCHICAL IMPLEMENTATION
For completeness, pseudocode for the hierarchical implementation of DE-FEM is provided. This algorithm assembles a leaf integration element by looping over the hierarchical tree and assembling all contributions. This pseudocode works in all dimensions and is based on the pseudocode given in the work of Zhang et al. 55 In Algorithm 1, the general assembly loop for integration elements is described. This procedure is the same as for standard FEM, except that a check is done to confirm that the integration element lies within the physical domain. To assemble the contributions of integration element e (k) inside the physical domain at hierarchy level k, local arrays are initialized. Then, in a loop over integration points, the Jacobian determinant j, the shape/enrichment function matrix N, and its derivative B are computed as described next.
Algorithm 2 describes the computation of N and B at every integration point. As a first step, the nodal coordinates of the integration elements are obtained. Then, in a loop over the enriched nodes of the integration element, the contribution is found for each level in the hierarchy. Starting from the leaf element, we obtain the global coordinate, the Jacobian matrix, the Lagrange shape functions, and their derivatives. If the current hierarchy level is a leaf, the Jacobian determinant is also computed and stored. Furthermore, the enrichment function and its derivatives, corresponding to the current enriched node, are appended to an array. Finally, the level of hierarchy is decreased, the parent of the current element is found, and an inverse mapping is conducted to find the local coordinate in the parent element. Once all hierarchical levels and enriched nodes have been visited, the shape and enrichment functions and their derivatives are expanded, according the number of DOFs per node, into the N and B matrices