Interpolation-based immersogeometric analysis methods for multi-material and multi-physics problems

Immersed boundary methods are high-order accurate computational tools used to model geometrically complex problems in computational mechanics. While traditional finite element methods require the construction of high-quality boundary-fitted meshes, immersed boundary methods instead embed the computational domain in a background grid. Interpolation-based immersed boundary methods augment existing finite element software to non-invasively implement immersed boundary capabilities through extraction. Extraction interpolates the background basis as a linear combination of Lagrange polynomials defined on a foreground mesh, creating an interpolated basis that can be easily integrated by existing methods. This work extends the interpolation-based immersed boundary method to multi-material and multi-physics problems. Beginning from level-set descriptions of domain geometries, Heaviside enrichment is implemented to accommodate discontinuities in state variable fields across material interfaces. Adaptive refinement with truncated hierarchical B-splines is used to both improve interface geometry representations and resolve large solution gradients near interfaces. Multi-physics problems typically involve coupled fields where each field has unique discretization requirements. This work presents a novel discretization method for coupled problems through the application of extraction, using a single foreground mesh for all fields. Numerical examples illustrate optimal convergence rates for this method in both 2D and 3D, for heat conduction, linear elasticity, and a coupled thermo-mechanical problem. The utility of this method is demonstrated through image-based analysis of a composite sample, where in addition to circumventing typical meshing difficulties, this method reduces the required degrees of freedom compared to classical boundary-fitted finite element methods.


Introduction
Computational modeling has become an integral part of engineering of all types, and developments in manufacturing and design have increased the complexity of the models required.Multimaterial problems are now ubiquitous, appearing in the design and analysis of composites [56], advanced additive manufacturing products [49], and multi-phase system analysis [72].These new design spaces present challenges to modeling methods, primarily in the discretization of intricate geometries and domain interfaces, and the unique discontinuities in solutions that result from material interactions.Finite element methods (FEMs) are among the most widely used computational tools for structural analysis.Classical FEM relies on sufficiently refined boundary-fitted meshes to discretize both the geometric domain of a material and its solution's function space.The accuracy of these methods is closely tied to mesh quality [13], especially when high-order methods are employed [21].Thus, considerable care must be taken in mesh construction before analysis can be performed [40].As the complexity of multi-material problems increases, generating these high-quality conforming meshes required by FEM becomes increasingly challenging, especially in three dimensions.Even for single material problems, mesh generation and refinement can consume up to 80% of the design-through-analysis time of engineers [10].More complicated PDEs involving multiple state variables also present unique modeling challenges.In such multi-physics problems, the various fields may have distinct discretization needs and the fields must be coupled.
Immersed boundary methods circumvent conforming mesh generation by embedding the geometric problem domain into a background grid constructed on a geometrically simple domain., Initially proposed to track fluid-structure interfaces in [54], similar classes of immersed methods were likewise developed by the solid mechanics community to accommodate discontinuities in solution fields without remeshing to create boundary-fitted meshes.The partition of unity method (PUM), introduced in [4], leverages the concept of enriching solution functions using a priori knowledge of the location of discontinuities.This was combined with classical FEM in [65] and [64] to introduce the generalized finite element method (GFEM).A similar enriched method known as the "eXtended" finite element method (XFEM) [8,9,47] was also introduced to model crack propagation and other discontinuous problem.As opposed to the a priori knowledge used to place enrichment functions used by GFEM, XFEM is characterized by adaptive enrichment schemes.Immersed boundary methods have also been extended to include high-order methods.The finite cell method, [53,57] utilizes p refinement of Lagrange polynomial basis functions to increase error convergence rates within an immersed framework.In the field of meshfree methods, the concepts of immersed or embedded methods have been used to model heterogenous materials [60], and to enhance solution accuracy and stability near material interfaces in fluid-structure interaction problems [30].
Isogeometric analysis (IGA) directly utilizes the geometric representation used in most computer-aided design (CAD) software in analysis [31] and was initially proposed to address the issue of generating high-quality boundary-fitted meshes.IGA has also been enhanced through the application of immersed boundary methods for the analysis of 'trimmed' CAD geometries [15].The Bspline basis functions used in IGA offer additional advantages over classical FEM including improved geometric representation, higher levels of continuity, and improved per-degree-of-freedom accuracy compared to more common nodal finite element basis functions [32,33].These attributes make B-splines attractive choices for certain applications, including Kirchoff-Love shell analysis [38].The combinations of IGA and immersed boundary methods are often called immersogeometric methods [37], and have been validated for use with hierarchically refined T-spline CAD models [58].The XFEM class of enriched methods was applied to IGA to create the eXtended isogeometric method (XIGA) in [51], exploiting level-set geometric descriptions and sophisticated integration algorithms to solve complicated PDEs involving multiple materials.XIGA was extended to use truncated hierarchically refined B-splines (THBsplines) for multi-physics problems in [61].
While offering elegant solutions to the problem of modeling complex geometries, existing immersed boundary software is currently limited to custom research codes.This is in part due to the complexity of generating custom quadrature rules on each cut background element [18].These custom quadrature rules require implementations that depart significantly from the highly optimized integration algorithms used in classical FEM codes.For example, an immersed boundary functionality called MultiMesh [34] was developed for the popular open-source FEM software FEn-iCS [1], however it proved too difficult to maintain and was not ported to the more recent FEniCSx [6].
Classes of immersed boundary methods have been developed to reduce the difficulty of integrating over cut cells, including approximate domain methods and specifically the shifted boundary method [41,42].The shifted boundary method maps the boundaries of a computational domain to a mesh conforming surrogate domain and has been extended to high-order methods in [3].Interpolation-based immersed boundary methods exist to address the same concerns.
This work utilizes approximate extraction to retrofit classical FEM codes to perform immersed analysis.Extraction was introduced to represent B-spline basis functions with Bézier polynomials for implementation of IGA [12].Bézier extraction was generalized to Lagrange extraction in [59].With Lagrange extraction, each spline in a B-spline function space is represented as a linear combination of Lagrange polynomials belonging to an interpolatory function basis {N i }, which satisfies the Kronecker delta property δ ij = N i (x j ), where x j are the nodal coordinates.Because of this interpolatory property, Lagrange extractionbased methods are also called interpolation-based methods.IGA and other partition of unity methods are challenging to implement as shape functions are not the same from element to element [46]¿ making interpolation-based methods attractive.Lagrange extraction has been utilized to implement IGA within existing finite element software in [35,36], which uses the software FEniCS, and in [68], which uses the software Code Aster.
In [24], approximate extraction was applied to immersed boundary methods to address the same integration challenges presented by IGA.In interpolation-based immersed boundary methods, a domain is embedded in a structured background mesh, which is used to define a background basis.A boundary-fitted foreground mesh is quickly generated by decomposing the cut elements of the background.The foreground mesh, which does not need to adhere to the usual mesh quality metrics [40], is used to define a boundary-fitted Lagrange polynomial foreground basis.The background basis is interpolated with the foreground basis, and the interpolated basis can then be used to solve PDEs with Galerkin-type methods.As the interpolated basis is represented as a linear combination of Lagrange polynomials defined upon a boundary-fitted mesh, the integration can be performed with classical FE software, as illustrated with FEniCS in [24].This work also employed Bspline background bases, combining immersed and isogeometric methods for an immersogeometric method.The interpolation-based immersogeometric method differs from previous extraction-based IGA methods in its use of approximate instead of exact extraction.With approximate extraction, the interpolated background basis is not everywhere equivalent to the actual background basis, resulting in additional interpolation error.The additional interpolation error is demonstrated to be bounded by the optimal method error, thus the method still yields optimal convergence rates.Easing exact interpolation constraints to permit approximate extraction allows for easier and more efficient implementation.
This work expands interpolation-based immersogeometric methods to multi-material and multiphysics problems.In the previous work [24], a single uniform background basis was used for all fields of a single material problem.Here, a framework is presented to interpolate multiple enriched hierarchically refined background bases and extract these bases to a single foreground mesh.Hierarchical refinement allows for local refinement around material interfaces, which are described by level set functions.Heaviside enrichment applied at material interfaces allow the background bases to approximate discontinuous state variable fields.Within this new framework, separate background bases may be used to approximate different fields in multi-physics applications, which are interpolated using a single foreground basis, allowing for easy coupling.This framework is implemented in the next generation open-source software library FEniCSx [6], with code available at [23].Multimaterials heat conduction and linear elasticity are modeled to demonstrate the error convergence rates of the proposed workflow.A coupled thermo-mechanical problem illustrates the combined multi-material and multi-physics capabilities of the interpolation-based immersed boundary framework.
The outline of this paper is as follows: Section 2 provides an overview of classic B-splines and the generation of the truncated hierarchically refined B-splines (THB-splines) used for local refinement.Section 3 describes this method's treatment of the geometric description of material interfaces and the Heaviside enrichment of solution spaces with discontinuities at material interfaces.Section 4 details the novel interpolation-based immersed boundary method's application to multi-material and multi-physics problems and its implementation workflow within existing FEM codes.Section 5 will provide numerical results validating and expanding upon this method, and finally Section 6 will draw conclusions and indicate future work to be done with this method.

Hierarchical B-splines
Immersogeometric analysis combines aspects of two classes of methods, immersed methods and isogeometric methods.As in isogeometric methods, this work employs splines to represent state variable fields.The following is an overview of spline fundamentals, starting with a review of multivariate B-splines, proceeding with a discussion of hierarchical refinement, and concluding with the process for generating truncated hierarchically refined B-splines (THB-splines).

A review of multivariate B-spline function spaces
Recall first the construction of a 1D univariate Bspline basis with polynomial order p: {B i,p (ξ)} n i=1 , where n is the number of basis functions.The domain is discretized with a knot vector Ξ = {ξ 1 , ξ 2 , ..., ξ n+p+1 } such that {ξ i } n+p+1 i=1 ⊂ R and ξ 1 ≤ ξ 2 ≤ ... ≤ ξ n+p+1 .The functions are then constructed recursively from the piecewise constant basis function using the Cox-de Boor recursion formula [11] If no interior knots are repeated, the basis is C p−1continuous at each knot in the interior domain, and C ∞ -continuous between the knots.The basis will form a partition of unity if the first and last knots are repeated p + 1 times.Higher dimension basis functions can be constructed by applying tensor-product operations to the univariate functions, such that where d p is the parametric space dimension, and there are d p knot vectors Ξ m = {ξ m 1 , ξ m 2 , ..., ξ m nm+pm+1 }, where n m is the number of basis functions and p m is the polynomial order of the m th parametric direction.Here i = {i 1 , ..., i dp } is a multi-index and p = {p 1 , ..., p dp } is the vector of polynomial degrees.The B-spline basis of the set of these functions is denoted by B p := {B i,p }.

Hierarchically refined B-splines
Hierarchiacally refined B-splines (HB-splines) are constructed using nested sequences of spline spaces created by repeated knot insertion.Following the algorithms presented in [25], an HB-spline basis begins with the construction of a sequence of r tensor-product spline spaces V l , each of which has an accompanying B-spline basis B l , and tensor-product Cartesian mesh K l , where elements are denoted by K.A sequence of subdomains Ω l are chosen, such that Ω l−1 is a subregion of Ω l , and each Ω l can be discretized with mesh elements K ∈ K l .Here r is the depth of refinement.The HB-spline basis H := H r−1 is then constructed recursively by the algorithm: In essence, each subsequent level's basis H l+1 is formed from the union of the set of basis functions from the previous level whose support is not in the new level's subdomain Ω l+1 , and the set of functions from the new basis B l+1 whose support is within Ω l+1 .This is illustrated with a 1D mesh in Figure 1a.
The hierarchically refined basis H is associated with a hierarchically refined mesh K, defined as

Enforcing the partition of unity property through truncation
While a useful tool for applying adaptive refinement to IGA, HB-spline bases violate the partition of unity (PU) property.To regain the this property, the hierarchically refined bases are truncated as in [25] and [26].In addition to forming a partition of unity, truncation reduces the size of some basis functions' supports, thereby reducing the bandwidth of the resulting system of equations when compared to a non-truncated HB-spline basis.
A given multivariate basis function B l ∈ B l can be represented as a linear combination of the more refined functions of level B l+1 : where c B l+1 B l are coefficients relating the coarse basis function B l to the finer function B l+1 .The truncation of B l removes the contributions from B l+1 ∈ B l+1 with support contained within Ω l+1 , such that Using a similar algorithm to that given in Eq. ( 6), the truncated basis T := T r−1 can be constructed by as illustrated in Figure 1b.

Immersed material interfaces
Workflows to solve multi-material PDE require functionalities to both describe the geometry of material interfaces and to represent the associated discontinuities in the state variable fields.In this work, level set functions (LSFs) are utilized to implicitly describe the geometry of material interfaces, and a generalized Heaviside enrichment strategy in conjunction with a set of interface terms is employed to represent the required discontinuities at material interfaces.

Representing interface geometry through level set functions
The level set method, developed in [52], has been used to describe interfaces in the extended finite element method (XFEM) [8,47] and extended isogeometric analysis (XIGA) [51,61].Following these works, the domain geometry is implicitly represented using LSFs ϕ i (x).An iso-level ϕ t of the LSF describes the interface Γ ± between two subdomains Ω + and Ω − such that With n LSFs, this method can represent up to 2 n subdomains.Materials are then associated with these subdomains using a multi-phase level set model as in [69], where phases are identified by phase indices P. Phase indices P are assigned with characteristic functions f i , such that Phases are then mapped onto material subregions.
In this work LSFs are discretized using linear basis functions from a THB-spline basis B k ∈ T , where ϕ k i are the coefficients associated with LSF ϕ i .The LSF are linearly interpolated such that the coefficients are the nodal values ϕ i k = ϕ i (x k ).This discretization is used to construct material indicator functions and to enrich background basis functions.

Heaviside enrichment of basis functions at material interfaces
Heaviside enrichment have been widely used in PUFEM [4], GFEM [65], and XFEM [9] as a means to represent strong discontinuities within elements.The enrichment strategies presented in most existing literature, such as in [29,67], adds enriched basis functions for each material.This work adopts the enrichment strategy presented by [51] which instead considers the material connectivity in the individual basis functions' supports to alleviate artificial numerical stiffening arising around small geometric features.This stiffening is caused by interpolating the state variable field in locally disconnected domains of the same material by the same basis function.The high-order, higher-continuity B-spline basis functions with large supports employed in this work, alongside the complex material layouts presented in Subsection 5.3, would exacerbate local stiffening effects and lead to an increased solution error if a global enrichment strategy were instead chosen.Additionally, function wise enrichment results in fewer added degrees of freedom than global enrichment, which reduces overall system size.For a given function B k with support supp(B k ), the phase IDs, defined in Eq. ( 13), are used to identify the For L k distinct subregions, the basis function B k requires L k enrichment levels.This enrichment is then achieved through indicator functions ψ m k , such that the enriched basis functions can be expressed as The enriched basis functions constructed from a single non-enriched basis function, are shown in Figure 2 for a two-material configuration.The biquadratic B-spline B k depicted in Figure 2(a) is enriched assigning one material to the inside and one material to the outside of the ellipse shown in Figure 2 2(d).allowing for the representation of strong discontinuities at the material interface.Interface conditions are enforced weakly; for example C 0 continuity can be enforced at the interface using Nitsche's method [2].

The interpolation-based immersed boundary method
One of the core challenges associated with classic immersed methods is the construction of custom quadrature rules for the various material regions within each intersected background element.Numerous solutions exist and have been used in custom research codes, such as octree refinement [19,37,57], interface reconstruction and tessellation [16,22,45], moment-fitting [48,66], and quadrature schemes using generalized Stokes theorem [27,28,63].However, the generation of such quadrature rules can generally not be implemented within existing finite element software without major changes to the software itself.
The interpolation-based immersed approach is instead designed to utilize the integration subroutines of existing FEM codes.To this end, a boundary-fitted foreground mesh is constructed with only minimal requirements on mesh quality.Element formation is then performed on the poor quality boundary-fitted mesh using existing standard finite element routines.Using Lagrange extraction operators [59] the resulting tangent matrix and force vector are projected into the enriched THB-spline space.This can be done either either on an elemental level during assembly, or globally afterwards.The resulting final problem uses an approximation of the enriched function space of the background mesh which is interpolated by the basis functions of the foreground mesh.
The following Subsection 4.1 will first introduce the thermo-elastic model problem.Subsection 4.2 will provide an overview of the interpolation-based approach specific to the multimaterial, multi-physics problems presented in this paper.For a more general and comprehensive introduction to the approach, we refer the reader to the authors' previous work on the topic [24].The generation of the boundary-fitted foreground mesh is discussed in Subsection 4.3.

Multi-material and multi-physics model problem
To illustrate the application of interpolationbased immersed boundary methods to multimaterial and multi-physics problems, a thermoelastic problem is introduced.This problem can be broken into a thermal subproblem and a structural subproblem.
Let a domain of interest Ω with closure denoted Ω be composed of n material subdomains A different thermal conductivity κ m may be associated with each material m for m ∈ M. A source term f : Ω → R, a boundary heat flux term q : ∂Ω → R on Γ q ⊂ ∂Ω, and Dirichlet boundary data T : ∂Ω → R on Γ T ⊂ ∂Ω are ascribed.The strong form for the thermal problem then reads as: where the intersections of the domain boundaries with the material subdomain boundaries.n denotes the surface normal.
The domain of interest Ω is embedded into a hierarchically refined background mesh K T , generated using the sequence of refined meshes K l and subdomains Ω l T . 1 Note that the largest subdomain Ω 0 T must be chosen such that the closure of the domain of interest is a subset, Ω ⊂ Ω 0 T .The mesh K T is associated with the enriched THBspline basis T T = {B T i }.The temperature field is discretized using the function space The discrete form can then be defined as: Find where R D T and R I T are Dirichlet and interface residual terms.The temperature Dirichlet residual is the result of Nitsche's method [50] enforcement of the Dirichlet boundary condition, Here the subscript (•) T refers to entities associated with the thermal subproblem.The subscript (•)u will refer to entities associated with the structural subproblem.T and u are used as superscripts when referring to entities that require subscripts for indices, such as B-spline basis functions (B T i and B u i ).
where β D T ≥ 0 is a user defined constant.The first integral of Eq. ( 22) will be negative for the symmetric version of Nitsche's method (which is employed in this work) or positive for the nonsymmetric version.h is taken as the characteristic element size on the foreground mesh, differing from the usual implementation where h is the element size on the background mesh.
The temperature interface conditions lines 2 and 3 of Eq. ( 19), are also enforced through a Nitsche-like method, resulting in the temperature interface residual where {•} = w i (•) i −w j (•) j is the weighted average of a given quantity.Motivated by the formulation in [2], these weights are defined as where h m is the characteristic size of the foreground element in domain Ω m bordering the interface facet, ω m is the characteristic material parameter, which for the thermal subproblem is κ m , and d p is the domain dimension.The penalty parameter γ ij T is defined as where β I T ≥ 0 is a user specified constant.The thermal subproblem can be stated compactly as the variational problem: Find where a T (T, θ) and L T (θ) can be computed from Eq. ( 21).The structural subproblem may utilize a differently refined background discretization.In this case a separate sequence of refined domains Ω l u may be selected.The structural subproblem's background mesh K u is constructed withthis sequence of subdomains and the same sequence of refined tensor-product Cartesian grids K l .The basis T u = {B u i } associated with the mesh is used for the components of the displacement field u.Each displacement component is discretized using the function space Using a similar derivation (given in full in Appendix A) as applied to the temperature subproblem the mechanical variation problem is compactly written as: Find The two subproblems are coupled through a constitutive model accounting for the thermal expansion by computing the mechanical strain ε m as where α is the thermal expansion coefficient and T 0 is the temperature in the reference configuration.I is the identity matrix.The fully coupled system is then given by the variational problem: where the form b(T, v) is a result of the coupling conditions and is given by

Interpolated basis functions
In traditional immersed boundary methods custom quadrature rules would be used to evaluate the integral in the weak form of the coupled problem, given in Eq. (30).The main idea of the interpolation-based immersed paradigm is to interpolate the background basis functions using a space of Lagrange functions defined on a foreground mesh, which can be integrated with classical quadrature methods.This workflow thus introduces an interpolated background function space for the thermal subproblem where the interpolated basis functions are defined as where is the Lagrange extraction operator.{N j } ν j=1 is the basis of a Lagrange FE space with nodal points x j such that N i (x j ) = δ ij .Here ν is the number of foreground basis functions.The same foreground space is used to interpolate the background bases for both the temperature and the displacement state variables.
The approximations of the temperature field becomes where {d T i } n T i=1 are the unknown coefficients and n T is the number of basis functions in the temperature field's interpolated background B-spline basis.The displacement is similarly discretized with vector valued function spaces as Details regarding this discretization are given in Appendix A.
The variational problem in Eq. ( 30) can be assembled using the interpolated bases in Eqs.(35) and (36) to form the linear system where are the matrix entries.
By consolidating the extraction operators for all components in a single matrix the linear system can be rewritten as where The quantities A θθ , A θv , A vv , b θ , and b v are evaluated and assembled on the boundaryfitted foreground mesh.As the foreground basis is thge typical conforming Lagrange polynomial basis, existing commercial or open-source FE software applying standard quadrature rules may be utilized.
The extraction operators can either be computed globally, as suggested in Eq. ( 45), or on an element level for more efficient implementation.However, the latter option requires a modification of the assembly maps in existing FE software.After solving the linear system in Eq. ( 44) with the interpolated basis, the solution is post-processed by projecting the solution for each field onto the foreground basis.
A visual example of the interpolation of an enriched B-spline basis function is given in Figure 2. The boundary-fitted foreground mesh in Figure 2(e) is used to define a discontinuous Lagrange polynomial foreground function space.The enriched background functions B 1 k and B 2 k are interpolated with this function space, as shown in Figure 2(f) and (g).

Foreground mesh generation
Both state variable fields are interpolated using the same Lagrange FE space {N j } ν j=1 which is constructed on a boundary-fitted foreground mesh.To generate this foreground mesh, the elements of a background mesh intersected by interfaces are decomposed into triangles and tetrahedrons whose facets approximately reconstruct the interfaces.
To ensure the foreground mesh is sufficiently refined, a background mesh for decomposition K D is constructed from the series of hierarchically refined meshes K l , where Ω l D ⊇ Ω l T ∪ Ω l u is a series of subdomains containing the union of the subdomains used for each subproblem's discretization.This ensures that K D is at least as refined as the meshes used for each subproblem's discretization and allows for additional refinement of the foreground to improve geometric resolution.The construction of a decomposition mesh K D is illustrated in Figure 3.
The decomposition mesh K D is triangulated by first applying a pre-defined triangulation to the intersected background elements, as shown in Figure 4 (a-b).This pre-defined triangulation forms 4 triangular elements in 2D domains or 24 hexahedral elements in 3D.Through root finding along elemental edges, the location of the isocontour is found, indicated by the black dots in Figure 4 (b).A subdivision template is then applied to each intersected triangle/tetrahedron, as shown in Figure 4 (c), to further subdivide the triangles/tetrahedrons into a set of triangles and tetrahedrons whose facets follow the LS isocontour.The last step is repeated recursively for each individual LSF ϕ h i which enables sharp geometric corners and edges to be captured where multiple interfaces meet.This approach has been used previously by, e.g., [62].The resulting approximation of the interface is piecewise linear and depends on the resolution of the decomposition mesh K D .
Note that the sliver elements and poor aspect ratios in meshes produced by this method are still suitable for the interpolation of the background basis functions which is not bound by the typical mesh quality constraints of traditional FEM [39].The resulting foreground mesh is of mixed element type (triangles and quadrilaterals in 2D, or hexahedrons and tetrahedrons in 3D) and contains hanging nodes.To accommodate the hanging nodes on the foreground mesh and to adequately interpolate the discontinuous enriched background basis functions described in Subsection 3.2, this method employs discontinuous Galerkin type elements for foreground function spaces.As the original THB-spline background basis maintains at minimum C p−1 continuity within each material domain, the continuity of the interpolated basis is likewise C p−1 -continuous where interpolation is exact [59].This work expands upon previous results using approximate extraction methods [24], where the constraints placed on the foreground Lagrange basis are lessened while numerical accuracy is maintained.
The sliver elements resulting from this procedure do not themselves present problems with the interpolation-based workflow.However, due to the arbitrary location of material interfaces with respect to cell boundaries, it is possible ⊆ Ω l u , shown in the third row, are defined for each state variable.For the decomposition mesh the series of subdomains Ω l D , shown in the last row, is defined such that the domain on each level l contains the union of the l th level domains for both each variable field.The hierarchically refined meshes K T , K u , and K D , shown in the rightmost column, are constructed with the mesh series K l in the top row and their respective subdomain series.for cells to be cut such that only a small portion of a basis function's support resides inside a given material domain.These small cell cuts result in sparsely supported basis functions, which can present issues with stability and linear conditioning.
Numerous strategies exist to mitigate these issues and were recently reviewed in [55].Strategies include basis function removal [20], ghost stabilization [14], and basis function agglomeration [5] or extension [15], and have the potential to be implemented within the presented interpolation-based framework.The benchmark problems presented in this work do not require special treatment of sparsely support basis functions and the authors leave the implementation of these stabilization strategies to future work.

Numerical Results
The accuracy of the proposed method is demonstrated through the study of several benchmark problems.Problems were defined in the Pythonbased open-source FE code FEniCSx, using foreground meshes and extraction operators generated by the open-source XIGA code MORIS, available at github.com/kkmaute/moris[44].The source code with which the results were generated is also available at github.com/jefromm/EXHUMEdolfinX [23].

Resolving discontinuities in solution fields through Heaviside enrichment
In this Subsection, a multi-material beam undergoing a spatially varying heat load presents a weakly discontinuous temperature solution field.In this work, weak discontinuities refer to discontinuities in the gradients of solution fields.The solution is approximated with an interpolated Heaviside-enriched background basis, using the interpolation-based immersed boundary workflow.The error convergence results from this study demonstrate the accuracy of this method's enrichment scheme for multi-material problems.Beams in both 2D and 3D domains are considered.The beam is initially defined with corner coordinates (0, 0) and (L, H) in 2D, and (0, 0, 0) and (L, H, H) in 3D, with L = 5 and H = 1.The beam is divided into 3 sections, with interfaces at x = L/4 and x = 3L/4.Each section of the beam is assigned a thermal conductivity κ 1 = 1.0, κ 2 = 0.1, and κ 3 = 1.0.To ensure the nonconformity of the material interfaces with respect to the background mesh facets, the 2D beam is rotated about the origin by angle ϕ = 20 • , while the 3D beam is rotated about y-and z-axes by angles ϕ y = −5 • and ϕ z = 5 • , respectively.The 2D beam is embedded into rectangle with corner coordinates (-1.0, -0.5) and (5,3) and the 3D beam is embedded into a rectangular prism with corner coordinates (-0.5, -0.25,-0.25)and (5.5,1.75,1.75).The geometric configurations are shown in Figure 5.Each edge (in 2D) or plane (in 3D) of the beam is defined by a level set function which extends beyond the beam domain in the mesh.The functions are extended through the entire mesh to fully resolve the corners (in 2D) or edges (in 3D) of the beam.
Thermal diffusion is governed by the Poisson equation.The strong and weak forms of this problem are given by Eqs. ( 19) and ( 21) in Subsection 4.1.The source term is constructed from the exact solution In 2D the beam-aligned coordinate is expressed in global coordinates as and in 3D where x = [x, y, z] are the mesh coordinates.The exact solution is imposed as Dirichlet boundary data on the ends of the bar, x ′ = 0 and x ′ = L.
For the 2D domains a suite of meshes with average background element sizes h ∈ •[1, 0.5, 0.25, 0.125, 0.0625] was used.In 3D, the background element sizes h ∈ 1.16 • [1, 0.5, 0.25, 0.125, 0.0625, ] were used.The foreground meshes were constructed with the workflow described in Subsection 4.3.Results from the Fig. 5: The geometric configurations for the 2D (left) and 3D (right) domains.A three-material beam is embedded in a structured background grid and rotated such that material interfaces do not align with element edges (in 2D) or facets (in 3D) .Elements intersected by the level set functions defining the beam geometry are triangulated to form a boundary-fitted foreground mesh.Fig. 6: The temperature solution field plotted for both domain geometries, using a bi(tri)-quadratic B-spline basis interpolated with a bi(tri)-quadratic foreground Lagrange space.The solution is weakly discontinuous at material interfaces.
2D mesh with background element size h = 0.5 and from the 3D mesh with background element size h = 0.25 using bi(tri)-quadratic Lagrange foreground bases to interpolate a bi(tri)-quadratic enriched B-spline bases are shown in Figure 6.Error convergence rates are plotted for both bi(tri)-linear and quadratic enriched background spline spaces in Figure 7, interpolated with equalorder foreground bases.Ideal error convergence rates validate the interpolation-based immersed boundary workflow for multi-material problems.

Approximating curved geometries through local foreground refinement
A major challenge in the modeling of multimaterial problems is the discretization of material interfaces.In this work, local refinement of the foreground mesh is performed to increase geometric resolution without affecting the number of degrees of freedom in the solution space.In this example the linear elastic behavior of an infinite plate with an embedded circular inclusion of radius R = 0.5 is modeled, with local refinement employed to improve the approximation of the inclusion geometry.The inclusion is comprised of Material 1 with Lamé constants λ 1 = 497.16and µ 1 = 390.63,while the exterior plate is made of Material 2 with λ 2 = 656.79,and µ 2 = 338.35.A uniform isotropic eigenstrain of ε 0 = 0.1 is imposed on the inclusion.The weakly discontinuous analytic solution for the radial displacement is given in [70] as where Exploiting symmetry, only the upper right quadrant of the plate is modeled as shown in Figure 8. Symmetry conditions are enforced on the left and bottom edges of the domain, and the exact displacement is prescribed on the right and top edges.The solution domain is a 5 × 5 square, with a quarter circle at the lower left corner.
The strong and discrete forms of the multimaterial linear elasticity PDE are detailed in Appendix A. For this example the mechanical Fig. 8: Approximated radial displacement of the eigenstrain problem.The background mesh is shown in white, while the foreground integration mesh is overload in black.The background mesh remains the same as local refinement is applied to the foreground for improved geometric resolution.The images correspond to the coarsest refinement level with background element size h = 0.625.Fig. 9: Error convergence data for the eigenstrain problem, illustrating the efficacy of foreground refinement.With no foreground refinement, the L 2 error convergence rate is limited by the geometric error.With sufficient foreground refinement (3x LR), the convergence rate approaches the ideal of 3.
The foreground meshes used here are locally refined about the material interfaces, as seen in Figure 8.The background basis functions remain constant, meaning that there is no increase to the number of solution degrees of freedom with foreground refinement.
The convergence plots in Figure 9 show the expected error convergence rates for the bi-linear B-spline function space, while the L 2 error convergence of of the bi-quadratic function space is limited to the second order by the geometric error.This is confirmed by the reduction in error magnitude with local refined of the foreground mesh around the curved interface.With sufficient local foreground refinement the convergence rates approach the ideal rate.
The foreground meshes generated using this octree local refinement strategy include hanging nodes, which are difficult for most commercial or open-source finite element software to handle.To accommodate these hanging nodes this interpolation-based workflow employs discontinuous foreground Lagrange polynomial basis functions.Classical discontinuous Galerkin methods require augmentation of the variational form to enforce higher levels of continuity at cell interfaces [17], but this augmentation is not done within this workflow.The results in this work indicate that an interpolated basis of lower continuity than required of a basis used for conforming finite element methods can be employed in this method, provided the background basis meets the continuity requirements.This attribute of interpolated-based methods was first exploited in [24] to approximate fourth-order PDEs with a background basis of quadratic B-spline bases, which are C 1 continuous, interpolated with a foreground basis of quadratic Lagrange polynomials, which are only C 0 continuous.

Image-based thermo-mechanical analysis of composite materials, utilizing multiple levels of local refinement
The capability of this method to tackle distinct discretization requirements of state variable fields within a multi-physics problem is illustrated here in a coupled thermo-elastic problem.This problem is posed on an alumina-epoxy composite sample undergoing simultaneous heating and loading conditions.Separate background discretizations are interpolated with a single foreground discretization for the two fields.
The geometry of the composite sample is taken from real micro-CT images converted to an implicit level-set description.The micro-CT image, taken from [71] and shown in Figure 10a, is made up of 200 × 200 pixels, with a pixel size  The process described in Subsection 4.3 was used to construct a foreground integration mesh, beginning with a uniform 80 element by 80 element axis aligned decomposition mesh.Two levels of local refinement were applied about the material interfaces, and the refined decomposition mesh was triangulate.The resulting foreground mesh, shown in 11a.contains 81,809 cells.
For comparison purposes, a similar workflow was used to generate a boundary-fitted mesh for use in classical FEM.Classical FEM with FEn-iCSx requires a single element type mesh without hanging nodes.To avoid hanging nodes and to sufficiently resolve the gradients of the state variable fields the decomposition mesh was uniformly refined three times forming a 640 element by 640 element square mesh.The cut cells were triangulated to create a boundary-fitted mesh, and then the remaining quadrilateral elements were triangulated.The mesh was then modified with the software package Coreform Cubit 2023.11 to improve mesh quality metrics.The resulting mesh, shown in Figure 11b, contains 1,675,860 cells.The method described here was used to ensure the level set descriptions of the material interfaces matched between this mesh and the foreground mesh used with the interpolation-based immersed boundary method.
With this sample, a heated compression-shear test was simulated.The top and bottom displacements were imposed as u top = [−0.01,−0.01]mm and u bottom = [0, 0]mm.The temperature at the top and bottom edges were specified as T top = 0 o C and T bottom = 100 o C. The environmental temperature was T 0 = 0 o C. The sides were left both traction and heat flux free.
The governing equations for heat conduction and linear elasticity were coupled via the constitutive law adding a thermal expansion component added to the total strain as in Eq. ( 29).Dirichlet boundary conditions for each field were imposed using Nitsche's method.Nitsche's terms enforcing continuity in T and u were also imposed upon the allumina-epoxy interfaces.
Bi-linear B-splines are used for the temperature field and bi-quadratic B-splines are used for each component of the displacement field, and both fields are interpolated using bi-quadratic Lagrange polynomials.For the boundary-fitted FEM comparison, Dirichlet boundary conditions were likewise enforced using Nitsches method, a bi-linear Lagrange basis was used for the temperature field, and a bi-quadratic Lagrange basis was used for the displacement.The results for the temperature field, the temperature gradient magnitude, the displacement magnitude, and the strain magnitude are shown in Figures 12, 13, 14, and 15, respectively.Within the figures, the foreground meshes are drawn in black while the background meshes are overlaid in white.The same foreground mesh, with two levels of local refinement to resolve the composite geometries, is used for each discretization.
Progressive levels of local refinement were applied to the background B-spline function spaces, as seen in the rows of Figures 12,13,14,and 15.This local refinement was implemented with truncated hierarchically refined B-splines, as described in Subsection 2.3.The regions of refinement with high-order splines are larger than those for the low-order to adequately support the truncated basis.This can be seen by comparing the white background meshes used for the temperature field depicted in Figures 12 and 13 with those used for the displacement field in Figures 14 and  15.The degrees of freedom associated with the various background discretizations are shown in Table 1.
With two levels of local refinement, the results were in qualitative agreement with the ones of the uniformly refined classical FEM example with far fewer degrees of freedom.As previously noted, the uniformly refined meshes were initially generated with the same workflow used to create the boundary conforming foreground meshes, a 'hands-off' method requiring only a bitmap image file.A more efficient mesh with fewer elements could be created for use in classical FEM but only with either significant user intervention or with sophisticated meshing software.Additionally, the workflow used here can easily be extended to 3D image stacks, which are not easily processed for classical FEM.

Conclusions
This work presents a new interpolation-based immersed boundary method for multi-material and multi-physics problems.This method employs enriched truncated hierarchically refined B-spline background spaces and discontinuous hierarchically refined Lagrange integration spaces.Domain geometry and material interfaces are represented by level set functions, which can be generated by, for example, geometric primitives or from 2D or 3D images.The domain is embedded in a grid with an associated B-spline basis and the level set function is discretized and used to compute Heaviside indicator functions.Basis functions are inspected and individually enriched if their support is intersected by domain boundaries.This function-wise enrichment offers large savings in system degrees of freedom when compared to previous published on global enrichment strategies.
Subdivision is used to generate a sequence of refined B-spline function spaces and associated     tensor-product meshes from the original function space and mesh.In this work, the meshes are refined along domain interfaces.A hierarchically refined B-spline function space is defined recursively, along with its associated hierarchically refined mesh.This hierarchical function space is then truncated to form a partition of unity.
The refined basis functions can then be enriched using the Heaviside enrichment functions when intersecting material interfaces.The enriched and refined background bases are interpolated by a discontinuous foreground basis.This foreground basis requires a boundaryfitted mesh, but is not subject to usual mesh conditioning constraints.The foreground meshes are thus constructed by the discretized level set function domain boundary descriptions and the hierarchically refined tensor-product meshes.Cells in the background mesh intersected by the domain and material interface geometry are triangulated to form a mixed-element type foreground mesh, which can be used by classical finite element codes to define a foreground basis.The level of refinement used to create the foreground mesh may exceed the level of refinement used for the background basis, allowing for greater geometric resolution without increasing the system's degrees of freedom associated with the background basis.
The background basis is interpolated using extraction operators.Extraction operators are constructed by evaluating the background basis at the locations of the foreground basis nodes.Existing classical finite element codes assemble the linear systems using the Lagrange foreground basis.The linear system is then projected onto the interpolated basis with the extraction operators and the system is solved with the interpolated basis.Interpolation allows the enriched and refined B-splines basis to be utilized in traditional finite element codes without the addition of complex integration subroutines, broadening the applicability of this method.
Several benchmark problems validated this method.Interpolated enriched bases were used in both 2D and 3D for a multi-material thermal diffusion problem with exact geometric representation.The numerical approximations were compared with analytic solutions and errors were computed.The L 2 and H 1 errors converged at ideal rates with mesh refinement.Foreground only refinement was implemented for a multi-material linear elasticity PDE involving a circular domain.When compared to the analytic solution, errors from a non-refined foreground mesh converged ideally with a bi-linear basis, but the L 2 error convergence rate of the bi-quadratic basis was limited to 2 nd order due to geometric error in the discretization of the domain boundary.With sufficient foreground refinement, the ideal 3 rd order convergence for the bi-quadratic basis was observed.Lastly, micro-CT images were used to generate a geometric discretization of an alumina-epoxy composite sample.Thermo-elastity was simulated, coupling separately discretized temperature and displacement fields through a thermal expansion component.Unique truncated hierarchically refined B-spline bases were used for each field, bi-linear for temperature, and bi-quadratic for displacement.The results were in qualitative agreement with the ones of a uniformly refined traditional boundary conforming finite element simulation.
In this work, the open-source code FEniCSx is used to demonstrate the efficacy of interpolationbased immersed boundary methods, and future work will expand implementation to other existing finite element codes.Within FEniCSx, additions to the interpolation-based immersed boundary workflow will be the implementation of stabilization techniques to address issues of linear conditioning.also acknowledge and thank Dr. David Kamensky for his assistance in developing the interpolation workflow.
(a) A 1D illustration of a sequence of subregions Ω l and quadratic B-spline bases B l .(b) The hierarchically refined mesh K, the HBspline basis H, the truncated basis T , and the verification of the partition of unity.

Fig. 2 :
Fig. 2: A B-spline basis is defined on a structured background mesh and an example function B k is depicted in (a).Using the geometry description of the material subdomains in (b), Heaviside enrichment is applied to form the discontinuous functions ψ 1 k B k and ψ 2 k B k , depicted in (c) and (d) respectively.A Lagrange foreground function space is defined on the boundary-fitted mesh in (e).The function space is used to interpolate the enriched background functions ψ1 k Bk and ψ2 k Bk , depicted in (f) and (g), respectively.
with k ∈ M and k ̸ = m, are the material interfaces, and [[•]] = (•) k − (•) m is the jump of a given quantity over an interface Γ km .The material fields are defined T m = T (x), x ∈ Ω m , and

Fig. 3 :
Fig. 3: A series of uniformly refined meshes K l are shown in the top row.A series of nested subdomains, Ω l+1 T ⊆ Ω l T , shown in the second row, and Ω l+1 u

Fig. 4 :
Fig. 4: Foreground integration meshes are formed by triangulating cut elements of the decomposition mesh K D .(a) A cell is intersected by the isocontour of the discretized level set function ϕ h = ϕ t defining a material interface.(b) The cell is subdivided into triangular cells and isocontour-edge intersections (indicated by black cicles) computed.(c) Using the intersection points as nodal points, the cell is further subdivided.

Fig. 7 :
Fig. 7: Ideal convergence rates are seen for both bi(tri)-linear and quadratic B-spline basis functions, which were interpolated with equal order foreground Lagrange function spaces.The convergence rates are ideal for both the 2D domain (top) and 3D domain (bottom).

Fig. 10 :
Fig. 10: Top:Micro-CT image of alumina-epoxy composite, where the white sections signify alumina particles and the grey is the surrounding epoxy.Bottom: Smoothed image used to generate the LSF geometric description

Fig. 11 :
Fig. 11: The whole domain is shown in the images on the left with the box indicating the region shown in the zoomed in view on the left.

Fig. 12 :
Fig. 12: The temperature field results are compared for four different mesh configurations.The left shows the entire domain with the box indicating the region shown in the zoomed in view on the right.

Fig. 13 :
Fig. 13: The the temperature gradient magnitude field is compared for four different mesh configurations.The left shows the entire domain with the box indicating the region shown in the zoomed in view on the right.

Fig. 14 :
Fig. 14: The displacement magnitude field results are compared for four different mesh configurations.The left shows the entire domain with the box indicating the region shown in the zoomed in view on the right.

Fig. 15 :
Fig. 15: The mechanical strain magnitude field results are compared for four different mesh configurations.The left shows the entire domain with the box indicating the region shown in the zoomed in view on the right.

Table 1 :
Degrees of freedom (DOFs) associated with the background function spaces at various levels of local refinement and with the uniformly refined FEM function spaces.