On the stability and interpolating properties of the Hierarchical Interface-enriched Finite Element Method

Abstract The Hierarchical Interface-enriched Finite Element Method (HIFEM) is a technique for solving problems containing discontinuities in the field gradient using finite element meshes that do not conform (match) the domain morphology. The method is suitable for analyzing problems with complex geometries or when such geometry is not known a priori. Contrary to the eXtended/Generalized Finite Element Method (X/GFEM), enrichments are placed on nodes created along interfaces, and a recursive enrichment strategy is used to resolve multiple discontinuities crossing single elements. In this manuscript we rigorously study the approximating properties and stability of HIFEM. A study on the enrichments’ polynomial order shows that the formulation does not pass the patch test as long as enrichments do not replicate the approximating properties of partition of unity shape functions. Regarding stability, we show that condition numbers of system matrices grow at the same rate as those of standard FEM—and without requiring a preconditioner. This intrinsic stability is accomplished by means of a proper construction of enrichment functions that are properly scaled as interfaces approach mesh nodes. We further demonstrate that, even without scaling, using a simple preconditioner recovers stability. The method’s stability is further demonstrated on the modeling of challenging thermal and mechanical problems with complex morphologies.


Introduction
While the finite element method (FEM) has become the standard numerical procedure for simulating a wide range of problems in solid mechanics, creating appropriate finite element (FE) meshes for problems with complex morphologies is still a laborious and time-consuming part of the modeling process. Moreover, because of the computational cost associated with creating FE meshes that match the domain geometry, this approach is particularly undesirable in iterative procedures where the geometry is not known a priori, such as in topology optimization or in problems with evolving discontinuities.
In order to avoid remeshing, several numerical techniques have been proposed over the years. In the Universal Meshes approach [1], the FE discretization is transformed so that finite elements conform to the description of discontinuities. This is accomplished by using an iterative smoothing procedure that modifies nodal locations while retaining the original connectivity table of the finite element mesh. In the Conforming to Interface Structured Adaptive Mesh Refinement (CISAMR) [2,3], a conforming discretization to discontinuities is created without the need of an iterative approach, but new elements are added to the mesh close to discontinuities to better characterize the field gradients. In the Conformal Decomposition FEM (CDFEM) [4,5], cut elements of the discretization are decomposed into conformal simplexes that conform to discontinuities. In the eXtended/Generalized FEM (X/GFEM) [6][7][8][9], a method that is based on the Partition of Unity Method (PUFEM) [10], discontinuities are resolved implicitly by means of an enriched formulation. In X/GFEM, the standard FE space is augmented by means of an enriched term that is able to reconstruct field (gradient) discontinuities in the approximate solution.
Enriched formulations provide an elegant solution to decoupling discontinuities from the underlying FE discretization. In this context, Soghrati et al. [19] introduced the Interface-enriched GFEM (IGFEM) to handle an arbitrary number of weak discontinuities-where the field gradient is discontinuous, as for example when analyzing problems with multiple material phases. Starting from the original (unmodified) finite element mesh, which will be henceforth referred to as the background mesh, this technique places enriched nodes at the intersection locations between weak discontinuities (e.g., material interfaces) and element edges. IGFEM yields accuracy and rates of convergence that are on par with standard FEM-but without the need for matching discretizations [19][20][21]. Therefore, IGFEM combines the best features of both standard FEM and X/GFEM; while decoupling geometry from discretization-X/GFEM most appealing feature-IGFEM also features a straightforward enforcement of essential (Dirichlet) boundary conditions (BCs) and mesh nodes retain the Kronecker delta property. IGFEM has been successfully used for modeling heat transfer in fiber-reinforced composites [20] and actively-cooled microvascular materials [19,20,22], and for modeling the multi-scale damage evolution in heterogeneous adhesives [23]. The mesh decoupling virtue of IGFEM while treating weak discontinuities within the problem domain was later extended to treat external boundaries [21], positioning the method as an immerse boundary or fictitious domain technique. Also, a generalization of IGFEM to handle strong discontinuities (e.g., cracks) was introduced recently by Aragón and Simone in the Discontinuity-Enriched Finite Element Method (DE-FEM) [24].
One of the main limitations of IGFEM is its inability in creating appropriate enrichment functions for elements cut by multiple interfaces. In the absence of an adaptive refinement scheme, the mesh size h has to be kept small enough to ensure every element is crossed at most by a single interface. This mandates for a priori information about the problem's discontinuities throughout the entire simulation, which is often not available as when interfaces evolve in time or when the final geometry is the result of an optimization procedure. Furthermore, when modeling problems where several material interfaces meet at a junction (e.g., polycrystalline materials), even refining the background mesh cannot avoid the presence of multiple intersecting interfaces within an element. To eliminate this restriction, Soghrati [25] introduced the Hierarchical Interface-enriched FEM (HIFEM), where a recursive algorithm is used for constructing enrichment functions in elements cut by an arbitrary number of interfaces. Another unique advantage of HIFEM is that material interfaces need not be aware of one another, which considerably facilitates the computer implementation of this method since enrichment functions associated with each interface can be computed independently. HIFEM was later extended to treat 3-D problems in [26] and also combined with higherorder elements in [27] for simulating the failure response of heterogeneous adhesives [28]. HIFEM yields the same accuracy and convergence rates of its predecessor [25,26], and because HIFEM is a generalization of IGFEM, the former acronym will refer to both methods henceforth.
Both HIFEM and its predecessor IGFEM suffer from the same stability issues as those suffered by X/GFEM, whereby global stiffness matrices become ill-conditioned when material interfaces pass arbitrarily close to nodes of the background finite element discretization. Because iterative solvers are faster than direct solvers as long as sparse systems of equations are well-conditioned, stability has been a matter of concern ever since our first publication on IGFEM [19]. To circumvent the issue of ill-conditioning, in [19] we first proposed the scaling of enrichment functions so their contribution becomes vanishingly small as their corresponding interfaces approach mesh nodes. The originally proposed scaling function failed to produce well-conditioned systems in some cases, so it was found later that another scaling function performed better [29]. Here we propose yet another scaling factor, which contrary to that found empirically by Soghrati et al. [29], is derived mathematically by minimizing the condition number Fig. 1. Solid composed of two phases Ω 1,2 separated by an interface Γ 12 . The boundary Γ , with outward unit normal n, is decomposed into Γ u and Γ t ; these are, respectively, regions with prescribed Dirichlet and Neumann boundary conditions. of a system matrix resulting from a 1-D analysis. We show this new scaling also works for higher-dimensional problems and we compare its performance to previous functions in [19] and [29].
In this article we thoroughly study the stability of HIFEM. We demonstrate that the condition number of the raw system matrix (immediately obtained after assembly) grows as O ( h −2 ) , which is the same rate as that of standard FEM with matching meshes. This work echoes ongoing research in pursuit of stable enriched formulations, for instance stable GFEM (SGFEM) [30][31][32]. Nevertheless, contrary to works on SGFEM, where stability is studied on a preconditioned system matrix, we show that a simple modification to HIFEM's enrichment functions immediately yields an intrinsically stable method. Further, we show that if enrichment functions are not scaledwhich results in an unstable formulation-applying the Jacobi-like preconditioner used in SGFEM recovers stability. Jacobi preconditioning has also been used in the context of CDFEM to remove the ill-conditioning associated with arbitrary element aspect ratios [4]. In this manuscript we also provide new insight on the interpolating properties of HIFEM. Although categorized as a technique that decouples discretization from discontinuities, it is still necessary to subdivide elements cut by material interfaces into smaller sub-domains to perform numerical integration and to evaluate enrichment functions. Finally, we investigate the approximation properties of HIFEM and show that the polynomial orders used in integration elements and their compatibility with those of the background mesh are crucial to the performance of HIFEM-i.e. to ensure passing simple patch tests and yielding optimal rates of convergence for problems that do not exhibit singularities.

HIFEM formulation
The domain of interest is defined as a region Ω ⊂ R d , as schematically shown in Fig. 1 for a 2-D problem. The domain boundary ∂Ω ≡ Γ = Ω \ Ω , with outward unit normal n, is composed by disjoint regions Γ u and Γ t , where Dirichlet and Neumann BCs are prescribed, respectively. The domain Ω is composed of disjoint subdomains-which represent different phases-such that Each subdomain Ω i has its own boundary Γ i with outward normal n i such that ∪ N i=1 Γ i \ Γ represent internal interfaces between subdomains. We use the notation Γ i j to denote the interface between the ith and jth subdomains, i.e.
Without further formalities, the weak (variational) form of a steady-state thermal boundary value problem is expressed as: Given the thermal conductivity κ, the heat source f , and the prescribed heat fluxq, find the temperature field u ∈ U such that where U and V are conveniently chosen function sets, where functions in the former satisfy the prescribed temperatureū on Γ u . Similarly, the weak form of the elastostatics boundary value problem reads: Given the body force b and the prescribed tractiont, find u ∈ U such that where σ = Cε is the Cauchy stress tensor, with C the constitutive tensor and ε (u) = 1 2 (∇u + ∇u ⊺ ) the linearized strain tensor. Function sets U, V are vector-valued, and as before u| Γ u =ū.
Henceforth, the symbol u will be used to denote the primal variable irrespective of the type of problem, but note that depending on the context (elastostatics) u may denote a vector. Thus, either of the problem statements above can be written in abstract form as: Find u ∈ U such that where the bilinear and linear forms can be inferred by inspecting, respectively, the left-hand and right-hand sides of Eqs. (1) and (2). As customary, we follow the Bubnov-Galerkin approach and (3) is discretized by selecting finite-dimensional function sets U h and V h , as discussed in the following section.

Hierarchical interface-enriched finite element discrete formulation
We discretize Ω into M non-overlapping finite elements e such that Ω , and int (·) denotes set interior. HIFEM does not require creating elements that conform to material interfaces and thus a structured mesh is the simplest discretization choice. Trial solution and weight functions are taken from the following enriched space: with ι h and ι e denoting index sets for mesh and enriched nodes, respectively, where the latter are created in HIFEM along material interfaces. In the first term,ũ i and N i represent the standard degree of freedom (DOF) and Lagrange shape function, respectively, associated with the ith mesh node. Similarly, in the augmented part, α i and ψ i represent respectively, the generalized (enriched) DOF and enrichment function associated with the ith enriched node. The parameter s i is a scaling factor that is incorporated in the formulation to improve the condition number of the global system matrix; this parameter will be discussed thoroughly in Section 4. As noted previously,ũ i and α i have dimensions commensurate to the underlying physics problem. Noteworthy,ũ i physically represents the primal field value at the ith node; this is because enrichment functions in HIFEM vanish at mesh node locations by construction, and thus standard mesh nodes retain the Kronecker delta property of standard FEM. It is worth noting that the enriched space given by (4) is devoid of enrichment functions that resolve singularities. Consequently, the formulation does not yield optimal convergence rates in problems where such singularities are present, e.g., interface junctions in polycrystalline microstructures.
The main advantage of HIFEM over IGFEM is the ability to handle multiple interfaces interacting within a single element. This is made possible via a hierarchical enrichment strategy, where enrichment functions are built recursively over the hierarchy of interfaces. Therefore, at the element level, the trial solution u h ∈ S e 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. In Section 3 we provide the required algorithmic considerations for recursively building the hierarchical shape functions ψ ki in Eq. (5). Following standard procedures, by choosing both the trial solution and the weight function from S e we arrive at the discrete form of (3). The global system of equations K U = R is obtained by assembling the local contributions of arrays obtained at element level. In other words, the global stiffness matrix and the global force vector are obtained, respectively, by K = A i k i and R = A i r i , where A denotes the FE assembly operator. It is worth noting that local arrays have the following structure: where subscripts u and α refer to standard and enriched parts of the approximation, respectively. In Section 4 we discuss the scaling factor s ki in Eq. (5) and its impact on the condition number of the global stiffness matrix K , which is obtained as where λ max and λ min are the maximum and minimum (non-zero) eigenvalues of K . The hierarchical nature of the method can be understood by looking at (5). The term "hierarchical" in HIFEM therefore refers to the way the enriched finite element space is constructed, where enrichment functions are added to the partition of unity shape functions over hierarchical levels to resolve the discontinuous field gradient. This hierarchical nature is somewhat similar to the p-version of the finite element method ( p-FEM) [33], where the linear partition of unity is kept intact regardless of the interpolant order. Indeed in p-FEM the linear shape functions are augmented by higher-order shape functions to build the final interpolant. Notice that as in p-FEM, in HIFEM the partition of unity property of the final interpolant is lost within cut elements-although because enrichment functions vanish at mesh nodes, the Kronecker-δ property still holds at these locations.

Interpolating properties of HIFEM
IGFEM, and subsequently HIFEM, were developed to resolve accurately the primal field (i.e., displacement or temperature) and with a discretization that is fully decoupled from the discontinuities' number and locations. Because discontinuities can cut mesh elements in any way, one would immediately ask whether integration elements with high aspect ratios influence the approximating properties of the method. This is particularly true for IGFEM/HIFEM since enrichment functions are built from Lagrange shape functions in integration elements. Although the integration elements' aspect ratios do not affect numerical quadrature, they do affect the recovery of gradient fields (e.g., the stress field) [34][35][36][37]. In Soghrati et al. [34] and Nagarajan and Soghrati [35] it was reported that because of badly shaped elements, the stress values at interfaces are overestimated. A detailed study of this issue in the context of immerse boundary (fictitious domain) problems by van den Boom et al. [36] showed the problem is not significant at boundaries where Dirichlet boundary conditions are enforced exactly. In the context of fracture problems, where HIFEM was extended to treat also strong discontinuities via the Discontinuity-Enriched Finite Element Method [24], the effect of aspect ratios was studied thoroughly by Zhang et al. [37].
In this article we discuss another important aspect of the approximating properties of HIFEM. As explained below, we look at how enrichment functions are constructed and what happens when enrichment functions do not have the same polynomial interpolation order as the partition of unity shape functions of the background mesh.

Enrichment functions for triangular elements
Fig. 2 schematically illustrates the hierarchic subdivision of a 3-node triangular (T3) element for the purpose of constructing enrichment functions ψ ki . We first subdivide the background mesh (parent) element cut by a material interface into subdomains (children), resulting in the first level of hierarchy ( Fig. 2(b)). During this process, other material interfaces that possibly intersect the parent element are neglected. Once the subdivision is completed, the resulting children elements serve as parents for the next material interface that builds the second hierarchical level ( Fig. 2(c)). This process is repeated until all interfaces that intersect with the original parent (root) element are accounted for, and this procedure can be used regardless of their number and location. The corresponding children hierarchy is stored in an ordered tree data structure-or better, an ordered forest, where every mesh element is a root in the tree.
With the hierarchy of integration elements in place, constructing the enrichment function ψ ki associated with the ith generalized DOF created at the kth hierarchy level-independent of other levels-is done recursively as where N e is the Lagrange shape function of child e created at hierarchy level k that has a unit value at enriched node i. In other words, the shape function ψ ki is constructed considering all n hl children elements created at the kth hierarchy level that are connected to node i. The enrichment functions given by (8) are linear and thus compatible with the interpolating order of their parent element. Using these enrichments, the convergence properties of IGFEM/HIFEM have been thoroughly studied in 4 and e (1) previous works for problems without singularities. For IGFEM, optimal convergence rates in L 2 and H 1 norms were demonstrated for 2-D straight and curved interfaces in the original manuscript [19], and for 2-D immersed boundary problems in [21]. For HIFEM, when dealing with both straight and curved interfaces in 2-D, optimal convergence rates in L 2 and H 1 norms were reported in [25].

Enrichment functions for quadrangular elements
The overall algorithm used for constructing children integration elements and evaluating hierarchical enrichment functions in quadrilateral elements is similar to that explained in the previous section. A material interface can split a quadrangle in: (i) two quadrangles; (ii) one triangle and one pentagon; and (iii) two triangles. The first case scenario can simply be handled by creating 4-node bilinear children elements. However, there is a caveat when dealing with the other two cases to ensure that the HIFEM approximation passes simple patch tests and yields optimal precision and convergence rates. This is illustrated schematically in Fig. 3 for case (ii) above. It can be seen that starting from the Lagrange parent interpolation, using 3-node triangular integration elements for evaluating enrichment functions fails to pass a simple single-element patch test. This is because the enrichment functions built with linear triangular elements lack the bilinear term that is present in the parent mesh element. In other words, while the Lagrange basis functions interpolating the field in the parent element involve the bilinear monomial, this term is absent in the linear enrichment functions of triangular children elements, resulting in a quadratic approximate field rather than the target linear response. One way to alleviate this situation considers incrementing the polynomial order of enrichment functions, as shown in Fig. 3(c), where 6-node quadratic Lagrange triangles are shown as integration elements. In this case, additional enrichments are placed to the nodes located at edge centers to increase the interpolating order of the enriched formulation.
This hypothesis is tested numerically with the simple patch test of Fig. 4(a), which consists of a plate with a prescribed temperatureū = 0 at the bottom edge andū = 1 at the top edge. Left and right edges are adiabatic. The plate is composed of two materials separated by a straight interface Γ 12 at x 2 = 0. The ratio between material conductivities is κ 1 /κ 2 = 10. This problem is then discretized by linear triangular and bilinear quadrangular elements, as shown in Figs. 4(b)-4(c), respectively. The different colors in the latter denote different elements (4 triangles and 2 quadrangles in Figs. 4(b)-4(c), respectively). Dashed lines delineate the sides of integration elements as well. Fig. 4(d) shows the results of the triangular mesh, where the temperature field has been extruded and colors represent the magnitude of the heat flux. It is shown that the magnitude of the flux is constant throughout the domain, indicating that this case indeed passes the patch test. Rigorously, the linear Lagrange interpolation in the parent element is matched by the linear enrichments built from 3-node linear children triangles. In the quadrangular mesh case of Fig. 4(e), using 3-node triangular integration elements for evaluating the enrichment functions cannot replicate the 4-node quadrilateral parent mesh element, as explained earlier. The patch test is recovered by using 6-node quadratic integration elements.  Normally, enriched nodes are not added when an interface passes exactly through mesh nodes, as the gradient in the field can be readily accommodated by the C 0 -continuity property of the underlying FE mesh. However, the aforementioned case (iii) is another situation that requires special consideration. In this case the material interface cuts the parent element diagonally (passing through two opposite nodes). Even though the interface does not intersect with the element edges, generalized DOFs are still needed to capture the corresponding weak discontinuity. Although this special case scenario rarely occurs in practice, it cannot be ignored in a robust HIFEM implementation.

Stability of HIFEM
HIFEM-and its predecessor IGFEM-suffer from ill-conditioned system matrices when nodes of the background mesh are located arbitrarily close to material interfaces. In this section we aim to shed more light on the impact of the scaling factor in Eq. (4) on the condition number of the resulting global stiffness matrix. The discussion that follows is in the context of the heat equation, but similar considerations can be drafted for elastostatic problems. Consider two DOFs i, j associated with standard mesh nodes. The contribution of the parent (background) mesh element to the local (element-wise) stiffness matrix is evaluated as where integration is performed on the master domain □ parameterized by ξ ≡ (ξ, η), as it is traditionally done with isoparametric formulations. Also, J is the Jacobian matrix of this transformation, ∇ ξ N i ≡ [ N i,ξ N i,η ] ⊺ , and κ is the thermal conductivity, which is assumed constant within the element. Now consider two generalized DOFs m, n associated with enriched nodes in an integration subelement of the aforementioned parent mesh element. Assuming only one level of hierarchy, the contribution to the local stiffness matrix of the corresponding enrichment functions can be written as where △ is the master domain of the child element parameterized by the coordinate ξ c , J c is the Jacobian matrix of this new transformation, and s m and s n are the scaling factors assigned to each enriched node. If k mn ≫ k i j , the resulting local matrix would suffer from ill-conditioning. Comparing Eqs. (9) and (10) shows , as both are computed in the local coordinate system. Thus, what could cause a significant difference between components of the local matrix would be the Jacobian matrices associated with the parent and children elements, i.e., if cond ( . In order to minimize this effect, depending on the shape and aspect ratio of children elements, an appropriate scaling factor must counteract the effect of the Jacobian matrix to the matrix components associated with generalized DOFs such that To further clarify this point, consider for instance the case shown in Fig. 5, where a material interface cuts a 2-D element at a distance w ≪ h from a standard node. We look at the stiffness component k nn associated with enriched DOF α n , and without loss of generality, restrict the analysis to integration elements shown in the figure. Integration element e 3 has approximately the same area as the parent uncut element, and therefore J c 3 ≈ J. Referring back to Eqs. (9) and (10), not using any scaling factor would in any case yield k where k ii is any diagonal entry obtained from (9). For integration element e 1 , J −1 , which after substituting in (10) indicates that again not using any scaling yields k ii ). However, the behavior for integration element e 2 is quite different: Since J −1 i.e., its contribution to the stiffness matrix is much greater than that of the parent element as w → 0. Substituting s n = √ w in (10) would alleviate the problem and avoid an excessively large component in the stiffness matrix. The contribution of the three children elements shown in Fig. 5 is therefore k nn = k A similar reasoning can be followed to obtain the diagonal entry k mm associated with enriched DOF α m .
The issue of ill-conditioning, which is intrinsic to the method and also arises in other enriched formulations such as X/GFEM [31], has been a subject of discussion ever since the first article on IGFEM was published [19]. That original work proposed scaling the enrichment functions so that they vanish as interfaces approach mesh nodes. Referring to Fig. 6, the originally proposed scaling factor was  where 0 ≤ w i ≤ 1 is the relative distance from the enriched node to a mesh node (measured along their common edge). This factor did not seem to help much, so a modification was subsequently proposed in [29]: Here we propose a new factor to scale the enrichment functions: Similarly to those in (11) and (12), this new factor is also symmetric so the relative distance w i can be computed from either mesh node, which facilitates the computer implementation. In the following section we provide the rationale behind Eq. (13) and study its performance compared to previously proposed factors.

Optimal scaling in 1-D
Consider the 1-D cylindrical bar of length L depicted in Fig. 7, which is clamped on the left end and subjected to a force P at the free end. A perfect material interface located at x Γ = wL divides the bar into two parts, each with its own axial stiffness k i = E i A i , i = 1, 2, where E i and A i refer to the Young modulus and area of the ith part, respectively. The displacement at the right end can be analytically evaluated based on contributions from both parts, as shown in the figure.
This problem can be solved straightforwardly using standard FEM by placing a node at x Γ . Therefore, for the standard FEM approximation with only two elements, the global stiffness matrix and force vector are given, Fig. 7. Clamped bar subjected to a load P and composed of two parts with axial stiffnesses k 1 and k 2 . The material interface is located at x Γ = wL. The figure also shows the exact displacement at the free end. respectively, by and In this equation B i is the strain-displacement matrix for the ith element, r is the reaction force (so the first row corresponds to the prescribed DOF at the left end). The condition number for the stiffness matrix is given explicitly by We now model the same problem using HIFEM with a single enriched finite element, where the enrichment function captures the jump in the displacement gradient at x Γ . Let the trial solution be defined as where the enrichment function is given by With N = [ N 1 N 2 sψ ] and B = dN/dx, the stiffness matrix and force vector for this 1-D problem can be evaluated as and respectively, where k ≡ k 2 − k 1 denotes the jump in stiffness. For this simple problem, an optimal scaling factor can be obtained analytically by minimizing (7): This optimal scaling function does not only depend on interface location w, but also on material properties at either side.
The standard FEM condition number in Eq. (16) is now compared to that of the stiffness matrix in (19) using the different scaling factors shown in Fig. 8, which are plotted as a function of interface location w. The plots are color-coded to ease matching with their corresponding condition number curves in Fig. 9 (the red curve corresponds to no scaling, i.e. s = 1). The plot also shows previously proposed scaling factors (11) (in purple) and (12) (in orange). In addition, the proposed scaling factor in (13) is shown in green. This scaling corresponds to the optimal scaling (21) for the specific case k 1 = k 2 . For all other cases, for which k 1 ̸ = k 2 (in blue), as the axial stiffness of the right component gets increasingly higher, the scaling approaches the straight line √ 2N 1 (w), with N 1 (w) = 1 − w. Figs. 9(a)-9(f) show the condition number of the stiffness matrix (K ) as a function of relative interface location w for several values of stiffness ratio k 2 /k 1 . Without loss of generality, curves in the figures were created considering k 2 ≥ k 1 . The condition number given by (16) for standard FEM is shown in Fig. 9(a), where it can be seen that, as one of the elements at either side of the interface gets increasingly smaller, the condition number grows unbounded. For the HIFEM approximation that uses one enriched element, results are shown in Figs. 9(b)-9(f). Not using any scaling factor (s = 1) results in Fig. 9(b), showing the same unbounded behavior. For the firstly proposed factor (11), whose results are shown in 9(c), there is no improvement from the previous two plots. A substantial improvement is obtained for scaling factor (12), for the condition number does not grow unbounded when approaching the mesh nodes.
The minimum condition number is attained when using scaling factor (21) (c f . Fig. 9(e)). In these curves, the condition number approaches the value (K ) = 1 as the interface approaches mesh nodes, although its magnitude increases with increasing mismatch between stiffness values. Further, when the two materials are the same, the condition number is exactly 1 regardless of interface location (a case that would result in a zero enriched DOF as it would not be needed to represent the jump in the gradient). Because (21) is more complex than previously proposed factors-and thus harder to implement-we also studied the behavior of Eq. (21) simplified for no stiffness mismatch, i.e., k 1 = k 2 . The simplified expression, given by Eq. (13), produces a constant condition number regardless of interface location, which also increases in magnitude with increasing stiffness mismatch (see Fig. 9(f)).

Simple tests in 2-D
In this section we investigate the effect of previously studied scaling factors in a 2-D setting. An interface Γ 12 splits the problem domain into regions Ω 1 and Ω 2 , with corresponding Young's moduli E 1 and E 2 , respectively. We consider both a slight and a substantial mismatch between Young's moduli, a Poisson ratio ν = 0, and plane strain conditions. The condition numbers are reported as a function of the interface angle φ with respect to the horizontal direction given by e 1 . Two tests are then designed so that the interface traverses the computational domain in every possible way, as detailed below.

Rotating barycentric interface
This first test is designed so that the interface gets arbitrarily close to nodes of a single T3 element while keeping similar areas at either side of the interface. HIFEM adds four DOFs to the original six DOFs of the T3 element. This element has nodal coordinates x i = 1/3 3/4 e 1 ± 1/3 1/4 e 2 (i = 1, 2) and x 3 = −2/3 3/4 e 1 , a unit area, and its barycenter is located at x = 0. The results of this test are reported in Fig. 10, showing the condition number (K ) as a function of interface angle φ. Because of the problem's symmetry, only results in the range φ = [0, 2π/3] are presented. As reported earlier for the 1-D case, the results show that not using any scaling produces an   ill-conditioned stiffness matrix as the interface approaches nodes of the triangle. Both scaling factorss ands ⋆ , given respectively by Eqs. (12) and (13), recover stability-although the behavior of the newly proposed factor is slightly better. In addition, as shown previously in 1-D, an increased mismatch in Young's moduli shifts the curves (in this case two orders of magnitude).

Swiping interfaces
This second set of tests looks at vanishingly small areas on element assemblies. We consider an interface that swipes a computational domain in a way that the area of integration elements on one side of the interface tends to zero. The interface lies along the line that passes through a point q = −4/3 3/4 e 1 and has a varying slope that is given so that it swipes the domain from φ = −π/6 (for which µ (Ω 1 ) = 0, with µ (·) denoting measure) to φ = π/6 (with µ (Ω 2 ) = 0). Two sets of results are reported in Figs. 11(a) and 11(b), for assemblies of three and two elements, respectively.
For the three-element assembly, test results show the condition number remains bounded as long as either scaling parameters ors ⋆ is used. It is worth noting that the condition number remains bounded even when the areas of integration elements approach a zero value, for which φ = [−π/6, π/6]. In Fig. 11(b) we see that the condition Fig. 12. Model problem used to study the condition number of the stiffness matrix as a function of mesh refinement. A square of unit area is composed by two phases Ω 1 and Ω 2 , separated by either a straight (a) or a circular (b) interface Γ 12 .
number of the two-element assembly grows unbounded as φ → π/6, i.e., as the interface approaches the top node of the assembly. This is a pathological case that is encountered only when the interface slices the corner of the computational domain and is caused by the lack of support for the corresponding node. This is, nevertheless, a case that is rarely encountered in practice. Given their negligible impact on the domain morphology and corresponding HIFEM solution, such cases (interfaces merely chipping off domain corners) can be discarded to avoid ill-conditioning of the system stiffness matrix.

Effect of mesh refinement on structured meshes
In this section we study the stability of HIFEM on structured meshes as h → 0. Following the work of Kergrene et al. in the context of SGFEM [31], consider a square plate of unit area with two types of embedded interfaces (see Fig. 12). The first problem considers a straight interface Γ 12 passing through the point q = e 1 + e 2 at an angle −π/6 with respect to the plate's top edge. The second problem considers a circular inclusion with center point located at x c = 1/ √ 5e 1 + 1/ √ 3e 2 and radius r = 1/ √ 10. In both problems we approximate the steady-state conductive heat transfer, for which the weak formulation was given by (1). The ratio of thermal conductivity of Ω 2 with respect to Ω 1 is κ 2 /κ 1 = 10.

Raw system matrix
We first study the condition number of the system matrix just after the assembly, as done previously in Sections 4.1 and 4.2. This matrix can be written as where we made explicit standard and enriched components-together with their coupling terms. Our reference curve, denoted by (K uu ), is obtained by computing the condition number on the standard FEM component K uu . This curve, which scales as O ( h −2 ) with mesh size h, is close to what would be obtained by using standard FEM on matching meshes with the same mesh size.
Results of this study are summarized in Fig. 13, where the condition number is plotted as a function of the reciprocal of the mesh size h for both straight and circular interfaces (Figs. 13(a) and 13(b), respectively). While solid lines corresopnd to the condition number of the global system matrix K , dashed lines show the condition number of the enriched component K αα . The figures show that the scaling factor used in the formulation plays a fundamental role in the stability of HIFEM. While using no scaling (s = 1) results in a system whose condition number is not stable, using the proposed scalings ⋆ given by (13) recovers stability. As a result, the condition number of HIFEM grows as O ( h −2 ) . The figures also show that the newly proposed scaling factors ⋆ performs better than s, particularly for the circular interface, as results for the latter become increasingly unstable as meshes get finer.

Preconditioned system matrix
Results of preconditioning the system matrix after assembly are also presented. We study the condition number of a scaled system matrix [31]:  where the diagonal matrix D is defined such that D ii = 1/ √ K ii and thusK has unit values in the diagonal. As apparent from the results presented in Fig. 14, the scaling used is no longer important when applying this preconditioner, as the resulting condition number is the same. This has important implications, as a preconditioner could be used in iterative methods even when the formulation uses no scaling at all. In fact, as demonstrated by the condition number of the enriched portion of the stiffness matrix K αα , the curves for both scaling options overlap.

Numerical examples
Four numerical examples are presented in this section to demonstrate the application of HIFEM for modeling realworld thermal and structural problems with complex geometries. The condition numbers of the system matrix in each problem are also reported for each example with and without using the scaling factor in the HIFEM formulation.

Intervertebral disc mechanical response
In this first example problem, we use HIFEM to simulate the mechanical response of the 2-D model of a human intervertebral disc (IVD) subjected to a compressive load. IVD injury is one of the most common reasons of lower back pain, and therefore understanding its mechanical behavior under different loading conditions is crucial to design appropriate preventive strategies in work-related environments. Fig. 15 illustrates the 36 mm × 6 mm geometrical model of IVD studied in this example. As shown in this figure, IVD is composed of three main regions: a central gelatinous region known as the nucleus pulposus, the annulus fibrosis (outer region with higher stiffness), and endplates that connect IVD to adjacent vertebrae. Material properties for each phase are given in Table 1.
The HIFEM approximation of the normal stress field in the y-direction (σ y ) in the IVD model subject to a uniform compressive traction vector applied on the upper vertebra bone in depicted in Fig. 15. The simulation is conducted using an axisymmetric elastic constitutive model, as IVD geometry is axisymmetric with respect to any plane cutting in through-the-thickness direction (one of such slices is shown in Fig. 15(a)). The inset of Fig. 15(a) illustrates a small portion of the 200 × 100 structured background mesh (obtained by fitting a bounding box to the IVD geometrical model) used for creating this HIFEM model, together with children subelements created along intersecting material interfaces in the domain. Regarding stability, it is worth noting that condition numbers of the stiffness matrix for this problem without and with the proposed scaling factor are 7.82 × 10 9 and 3.45 × 10 3 , respectively. This comparison shows the effectiveness of such scaling in building a well-conditioned stiffness matrix.

Particulate composite micromechanical behavior
In this example, HIFEM is employed to approximate the magnitude of the microscopic displacement field   u 1   in six repeating unit cells (RUCs) of an Al/Al 2 O 3 particulate composite with different volume fractions, as shown in Fig. 16. The micromechanical model used for describing the constitutive behavior of the composite is described in [38]. The periodic RUCs studied in this example are virtually created using the reconstruction algorithm presented in [39]. Each RUC has a length of l = 70 µm and periodic boundary conditions along all its edges, and is subjected to a macroscopic strain ε 0 in the vertical direction. The Al matrix has an elastic modulus of E m = 300 GPa and a Poisson's ratio of ν m = 0.33, while mechanical properties of embedded Al 2 O 3 particles are E p = 68.9 GPa and ν p = 0.21. A small portion of the 200 × 200 background structured mesh used for creating the HIFEM model for all six RUCs, together with the resulting children elements along material interfaces, is illustrated in the inset of Fig. 16(d). This figure also illustrates the magnitude of the microscopic displacement field   u 1   in each RUC, showing the ability of HIFEM to capture weak discontinuities along material interfaces. Once again, for the six RUCs studied in this example, incorporating the proposed scaling factor in the HIFEM formulation leads to 6-7 orders of magnitude reduction in the condition number of the stiffness matrix-e.g., from 4.29 × 10 10 to 5.72 × 10 3 for the RUC with the highest particle volume fraction shown in Fig. 16(f).

Damage propagation in a heterogeneous adhesive
In this problem HIFEM is used to simulate the damage evolution in the 3-D microstructural model of a heterogeneous adhesive with an epoxy matrix and embedded silica particles. Material properties chosen for this simulation are E m = 2.0 GPa, ν m = 0.4 for the epoxy matrix and E p = 71.7 GPa, ν p = 0.17 for the silica particles. The microstructure of a 0.8 mm × 0.8 mm × 0.4 mm RUC for this system is illustrated in Fig. 17(a). The inset of this figure illustrates a small portion of the 160 × 160 × 80 background mesh used for discretizing the domain, together with the children elements created along the interface of one of the embedded particles. The HIFEM model created for this damage simulation has more than 2.5 million DOFs.
To simulate the damage response, periodic boundary conditions are imposed along lateral faces of RUC, while a macroscopic tensile displacement jump of u M = 0.4 mm is applied to its top and bottom faces. A multiscale cohesive model relying on an isotropic continuum brittle damage model, as described in [23,28,38], is employed to simulate the failure response of the adhesive layer. In this model, a damage parameter δ, ranging from 0 (intact) to 1 (fully damaged), is employed to reflect the impact of damage on the potential energy Φ m as whereΦ m is the potential energy in the undamaged material and ε m is the microscopic strain tensor. The evolution of damage in the epoxy matrix can be evaluated as where χ t is the internal state variable describing the magnitude of damage. Also, the damage surface G(Φ m ) is characterized using a Weibull distribution function given by In (26), Φ in m is the energy threshold for the initiation of damage, while the scalar parameters p 1 and p 2 determine the shape of the damage surface. For this simulation, the parameters for simulating damage were chosen as Φ in m = 152 J/m, p 1 = 25, and p 2 = 1.3. The HIFEM simulation of the damage pattern in the adhesive layer at the last time step (after failure) is depicted in Fig. 17(b), showing damaged zones (cracks) in the epoxy matrix perpendicular to the loading direction. Before the initiation of damage, the condition number of the stiffness matrix without applying any scaling factor was (K ) = 5.6 × 10 11 . With the proposed scaling, the condition number was reduced to 6.1 × 10 4 , showing once again the ability of the proposed scaling parameter to construct a well-conditioned system. Noteworthy, after damage takes place and especially near the failure point, the stiffness matrix becomes ill-conditioned even when using the proposed scaling. Nevertheless, this ill-conditioning is purely due to significant discrepancy in material properties of damaged and undamaged elements, which can only be improved (but not completely alleviated) using an element-deletion approach-i.e., by not including fully damaged elements during the assembly of the stiffness matrix.

Actively cooled aluminum plate
This example, which is completely detached from the stability study, showcases the application of HIFEM for simulating the steady-state thermal response of an actively-cooled aluminum plate with embedded microchannels, as illustrated in Fig. 18. This figure shows the lower half of these 50 mm × 50 mm × 2 mm 3-D plates, showing the increasing temperature of the coolant (water) pumped through microchannels, each with a diameter of ϕ = 750 µm, as well the temperature field in the surrounding aluminum. A uniform heat flux of q = 55 kW/m 2 is applied on the bottom face of each plate, while all other faces have a convective boundary condition with an ambient temperature u ∞ = 20 • C. The thermal conductivity of aluminum is κ = 205 W/(mK) and a heat transfer coefficient h = 100 W/(mK) is chosen, while the thermal conductivity of water is κ = 0.59 W/(mK).
The governing equations used for approximating the steady-state heat transfer of the actively-cooled plates, under the assumption of laminar Poiseuille flow with a fully-developed velocity profile in microchannels, can be found in [40]. Note that the stabilized Streamline Upwind Petrov-Galerkin (SUPG) method is utilized to approximate the hyperbolic convection-diffusion equations governing the heat transfer in the coolant phase. In the absence of active cooling, the maximum temperature due to the prescribed heat flux could theoretically reach u max ≈ 5.700 • C, which is well above the melting point of aluminum. This heat flux corresponds to an intense heat source (q = 55 kW/m 2 ) which can be generated by an actual computer processor.
The simulation results, which describe the HIFEM approximation of the temperature field obtained by discretizing the plate with a 150 × 150 × 8 structured mesh composed of tetrahedral elements, are summarized in Fig. 18. For the results shown in the top row, a flow rate Q = 10 ml/min is pumped through each of the four microchannels. In the bottom row, the sum of these flow rates (Q = 40 ml/min) is collected into a single embedded microchannel. Given the fact that the pressure drop in a microchannel has a direct relationship with its length and the flow rate, the total energy needed for the active-cooling of plates in the former case scenario is lower due to the shorter length of each channel. However, the results shown in Fig. 18 indicate that collecting the entire flow in a single channel is more effective in reducing the maximum temperature u max in the plate, with the rightmost sinusoidal design yielding the lowest value, i.e., u max = 57 • C.

Summary and conclusions
In this article we have demonstrated the application of the Hierarchical Interface-enriched Finite Element Method for modeling problems with weak discontinuities. The mesh decoupling feature of HIFEM gives unique flexibility to the method, as it allows the treatment of problems with complex geometries using structured finite element meshes. This is due to the hierarchical enrichment strategy used in HIFEM, which enables handling multiple material interfaces that are in close proximity or even intersecting with one another. HIFEM was showcased with a set of 2-D and 3-D problems with complex morphologies, where all approximate solutions were obtained on structured meshes. New formulation aspects presented in this manuscript showed that special care has to be taken when dealing with quadrilateral elements to ensure proper reproduction of the underlying mesh interpolant. That discussion extrapolates straightforwardly to the treatment of brick elements in 3-D (trilinear hexahedral elements).
In Section 4 we showed that an integration element's contribution to the stiffness matrix is inversely proportional to its Jacobian's determinant. Further, because the latter is proportional to the integration element's measure, entries to the stiffness matrix can grow unbounded if no scaling is used as an interface gets arbitrarily close to a mesh node. The new scaling factor proposed to stabilize HIFEM, which was derived analytically from a 1-D problem, not only solves the issue of ill-conditioning in 1-D but also shows remarkable performance in higher dimensions. Although a better performing scaling function could potentially be derived for 2-D/3-D, the proposed scaling yields stable systems because it only needs to scale appropriately with the integration elements' measure (area/volume in 2-D/3-D) to prevent ill-conditioning. This observation is in agreement with the ad hoc scaling parameter proposed earlier by Soghrati et al. in [29], which also scales as √ w for w → 0 and thus yields stable systems. With a comprehensive discussion and several tests we conclude that HIFEM is intrinsically stable, i.e., the condition number of the stiffness matrix grows at the same rate as that of standard FEM. In order to obtain well-conditioned system matrices, a scaling factor is used when constructing HIFEM enrichment functions, effectively acting as a preconditioner to the method. As a result, ill-conditioning is completely overcome when interfaces are placed arbitrarily close to mesh nodes and no preconditioner is needed when using an iterative solver. It is worth noting, however, that even when such scaling parameter is not used-in which case the resulting matrix becomes ill-conditioned-HIFEM is still stable by applying a simple Jacobi-like preconditioner.