Boolean finite cell method for multi-material problems including local enrichment of the Ansatz space

The Finite Cell Method (FCM) allows for an efficient and accurate simulation of complex geometries by utilizing an unfitted discretization based on rectangular elements equipped with higher-order shape functions. Since the mesh is not aligned to the geometric features, cut elements arise that are intersected by domain boundaries or internal material interfaces. Hence, for an accurate simulation of multi-material problems, several challenges have to be solved to handle cut elements. On the one hand, special integration schemes have to be used for computing the discontinuous integrands and on the other hand, the weak discontinuity of the displacement field along the material interfaces has to be captured accurately. While for the first issue, a space-tree decomposition is often employed, the latter issue can be solved by utilizing a local enrichment approach, adopted from the extended finite element method. In our contribution, a novel integration scheme for multi-material problems is introduced that, based on the B-FCM formulation for porous media, originally proposed by Abedian and Düster (Comput Mech 59(5):877–886, 2017), extends the standard space-tree decomposition by Boolean operations yielding a significantly reduced computational effort. The proposed multi-material B-FCM approach is combined with the local enrichment technique and tested for several problems involving material interfaces in 2D and 3D. The results show that the number of integration points and the computational time can be reduced by a significant amount, while maintaining the same accuracy as the standard FCM.


Finite cell method
The finite cell method (FCM) is an embedded domain approach based on an unfitted discretization and high-order shape functions [1][2][3], enabling an efficient and accurate simulation of geometrically complex structures including hole regions [4,5] and inclusions of additional material [6,7]. As point of departure, a general multi-material problem is depicted in Fig. 1a. Here, = 1 ∪ 2 is the physical domain of interest, ∂ its boundary, and D together with N B Márton Petö marton.petoe@ovgu.de 1 Institute of Mechanics, Otto von Guericke University Magdeburg, Magdeburg, Germany 2 Institute for Mechanics, Technical University of Darmstadt, Darmstadt, Germany the Dirichlet and Neumann boundaries, respectively. Finally, 12 is the material interface between 1 and 2 . In the FCM framework, is embedded into a larger domain e (Fig. 1b), giving rise to the fictitious domain fict = e \ of theoretically zero stiffness. Due to its simple shape, a straightforward discretization of e is possible by Cartesian meshes (Fig. 1c). Thus, the often cumbersome and error-prone generation of a geometry-conforming discretization by means of distorted finite elements can be avoided.
Note that due to the unfitted discretization, cut cells arise, that require extra attention regarding the implementation of boundary conditions [8,9], the condition number of the global system of equations [10][11][12], and the numerical integration of cell matrices (Sect. 1.2). For the fictitious domains, the volume integrals are penalized by the indicator function α = 1 in 10 −q in fict (1) (c) Finite cell discretization in order to avoid a large energy contribution of fict to the total system (for a more detailed mathematical formulation of the FCM involving void regions, see Sect. 4.1). In contrast to void regions, the displacement field along 12 is generally C 0 -continuous. Such displacement fields are poorly approximated by a linear combination of smooth shape functions over the cut cells. Regardless of the polynomial degree, the insufficiently captured kink in the displacement field causes severe oscillations in the strain and stress fields, which lead to a severely deteriorated accuracy of the simulation and suboptimal convergences rates [6].
In the extended finite element method (XFEM) [13][14][15], an accurate approximation of the displacement field for multimaterial elements is achieved by a local enrichment of the Ansatz space. Here, the key idea is based on introducing additional shape functions that are constructed to be C 0continuous along 12 . In the context of the FCM, Joulaian and Düster proposed the hp-d-FCM and hp-d/PUM-FCM approaches for dealing with material interfaces [6,7]. In both cases, the discontinuity is dealt with on a superimposed overlay mesh. In the former approach, the local behavior around the interface is captured by a geometry-conforming mesh, while in the latter one, the local enrichment of the XFEM is adopted. Note however, that in the hp-d/PUM-FCM, the unfitted overlay mesh is not necessarily aligned to the base mesh. Thus, it is possible to define larger overlay patches that span multiple cells in the enrichment zone. The weak form corresponding to the multi-material FCM in conjunction with the local enrichment approach is given in Sect. 4.2.

Integration of discontinuous functions
For an accurate FCM simulation, a proper integration of the cut cell matrices is of crucial importance. Since for discontinuous integrals, the integration accuracy of Gaussian quadrature rules is severely deteriorated, more suitable numerical integration schemes are required. Approaches being available for this purpose can be classified as follows: 1. A main branch is based on partitioning the integration domain into several integration sub-cells using quadtree/octree-decomposition techniques (QTD/OTD).

Another branch is represented by the Moment Fitting
(MF) approach, where in each cut cell, a unique quadrature rule is derived [26,27]. By pre-defining the position of the integration points, the computational complexity of this method can be significantly reduced [28][29][30][31]. Furthermore, by appropriate combination with an adaptive space-decomposition technique [32], or by using a nonnegative least square solver [33,34], application of MF to non-linear problems is also possible. Finally Düster and Allix proposed and investigated the combination of moment fitting with the local enrichment approach [35]. 3. A different idea is seen in the equivalent polynomials [36][37][38] method, which is based on the replacement of the discontinuous integrand by a continuous one yielding the same integral value. Thus, similar to the MF, the computationally often expensive space-partitioning can be avoided. 4. Further methods use the Divergence Theorem for reducing higher dimensional integrals to surface and line integrals [39]. Here, radial basis functions [40,41] or pre-derivation can be used [42] for evaluating the antiderivatives required by this method.
by the community nowadays [46][47][48][49][50]. However, if a high integration accuracy is desired, it requires a large number of integration points, having a significant impact on the computational time. In our previous articles, we introduced an effective approach based on merging the integration sub-cells both in 2D [51,52] and 3D [53] for reduced computational time.
In this contribution, the Boolean FCM (B-FCM) [1], whose aim is also the reduction of integration points, is investigated and extended for multi-material problems. The B-FCM for single material problems was also investigated in Ref. [51], where promising results were obtained. The B-FCM extends the integration procedure by Boolean operations and can be easily combined with standard spacedecomposition techniques. Originally, it was developed for single-material domains with hole regions. The advantage of the method is best demonstrated by Fig. 2, where Fig. 2a, depicts a cut cell, over which the function α f should be integrated. The standard FCM approach follows the QTDbased integration scheme as depicted in Fig. 2b, requiring 10 sub-cells. In the context of the B-FCM, the same integration accuracy can be obtained for the given example by employing 2 sub-cells only. This is achieved by first, integrating over the entire domain, where α is set to a value of 1 and second, by subtracting the contribution of the fictitious domain, which means that α is set to 0 in the physical part of the sub-cell and takes the value of 1 in the fictitious part. It was shown in Refs. [1] and [51] that by using this approach a significant amount of integration points can be saved.

Motivation
As mentioned in Sect. 1.3, the B-FCM was originally developed for problems that (i) solely consist of and fict , where (ii) the discontinuity in the integrals is caused by the piecewise constant indicator function α. On the contrary, in case of multi-material problems, the physical domain is composed by n d material sub-domains and the discontinuity in the integrand is caused by an arbitrary function D resulting from the enriched Ansatz space. Thus, for extending the single-material B-FCM to a multi-material version, two main steps are needed: (i) a generalization of the Boolean space-partitioning technique for decomposing the multi-material integration domain and (ii) an adjustment of the integrands over the Boolean sub-cells, which depend on the chosen enrichment type.

Simplifying integrals with Boolean operations
In this section, the key idea of the Boolean integration scheme is presented for embedded problems that contain several subdomains with arbitrary functions defined over them.

Two domains
Let ∈ R 1 be a union of the disjoint domains 1 and 2 . Over , we consider the continuous and discontinuous functions f (ξ ) and D(ξ ), respectively, where the latter one is defined as Then, the integral of the product D f over is composed by integrating over the sub-domains 1 and 2 as Assuming that D 1 is a well defined and known function not only over 1 , but also over , Eq. (4) can be reformulated, such that where the integral of the term D f is computed by the following two steps: First, the continuous function D 1 f is integrated over the entire domain and second, the integral value of the first term is augmented by (D 2 − D 1 ) f over 2 . Example: An example visualizing Eq. (5) is given in Fig. 3, where For simplicity, f (ξ ) = 1 is chosen, 1 and D(ξ ) defined as D 1 = sin(8ξ) in 1 and D 2 = ξ 2 in 2 . The total integral value of the function D f is depicted by the gray area in Fig. 3a. In Fig. 3b, the first integral on the right hand side of Eq. (5) is visualized. Finally, Fig. 3c depicts the contribution of the second term, where the red (green) color indicates regions that need be subtracted from (added to) the gray area in Fig. 3b in order to obtain the correct integral value.  Of course, there is nothing special about choosing 1 and D 1 as the base term in Eq. (5). The formulation works equally well for 2 and D 2 . In this case, however, it must hold that D 2 is defined not only over 2 but also over

Multiple domains
The formulation given in Sect. 2.1 can be also extended to higher dimensional domains composed by n d > 2 disjoint sub-domains according to Eq. (2) as well, where the discontinuous function D(ξ ) is defined individually over each sub-domain, such that Then, for a chosen sub-domain i, Eq. (5) generalizes to where in the second term, j = i is excluded due to the integral value being zero. Note that each sub-domain i can be chosen as the basis of the formulation, assuming that D i can be extended from i to . This is fairly simple for functions that are explicitly given, such as, piece-wise continuous material properties. However, for locally enriched multi-material problems, more steps are required for identifying D i over (see Sect. 5). Example: In Fig. 4, a two-dimensional example is given for Eq. (10), where a rectangular domain ∈ R 2 consists of three sub-domains (n d = 3). In the current example, the discontinuous function D over is a piece-wise constant function with the values D 1 = 2, D 2 = 4 and D 3 = 1 over the sub-domains 1 , 2 , and 3 , respectively ( Fig. 4a) Also in this example, f = 1 is selected for simplicity. Choosing the sub-domain i = 1 as the basis for the formulation, the first and second terms in Eq. (10) are depicted in Figs. 4b and 4c, respectively. It is easy to see that the integral values based on the direct integration over the individual sub-domains (I ) and on the Boolean formulation (I B ) are indeed equal Dj − D1 f dξ

Combination of Boolean operations and space-tree decomposition techniques
Generally, using Eq. (10) over Eq. (4) does not have major advantages, when i has a simple shape (e.g., quadrilateral, triangular, hexahedral, tetrahedral, or even circular domains [1]), such that the integrals can be computed either analytically, or via readily available numerical quadrature rules. However, when i has a complex shape, Eq. (4) is often computed using the QTD and OTD space-tree decomposition techniques. As it will be shown later, if Eq. (10) is combined with a suitable space-decomposition strategy, the Boolean approach is more beneficial than the pure QTD/OTD, as it leads to a smaller number of integration points, while maintaining the integration accuracy.

Standard space-tree decomposition
The computation of Eq. (4) using the QTD/OTD approach is based on a recursive decomposition of the integration domain into a disjoint set of n sc quadratic/cubic sub-cells = n sc k=1 ω k , where ω k is the domain of the kth sub-cell. The partition procedure is applied only to the cut sub-cells, leading to a set of integration domains with decreased size and increased density in the vicinity of the discontinuity, as depicted in Fig. 5 for the domains defined by Eqs. (11)- (14). The tree-depth is controlled by the refinement level R. Finally, the total integral value is computed by integrating the product D f via Gaussian quadrature over the individual sub-cells Each sub-cell has its own local coordinate system η ∈ R d , where d is the dimensionality of the problem. Furthermore, the geometry mapping to the parent cell is established by the transformation η = Q η→ξ (ξ ) involving a uniform scaling Fig. 5 Partitioning of the the parent cell via QTD with R = 3 and geometry mapping of the sub-cell ω k and translation. Since integrals in Eq. (17) are computed in the local space of the individual sub-cells, the appropriate change of integration limits has to be taken into account by 1] denotes the sub-cell in its local coordinate system η. Furthermore, both the Jacobian matrix J η→ξ = grad η ( Q η→ξ ) and its determinant are constants due to the nature of the geometry mapping. For sake of compactness, wherever it is possible in this article, the integrals over the sub-cells will be formulated w.r.t. dξ rather than to dη.

Boolean space-tree decomposition
The proposed Boolean integration scheme extends the conventional QTD/OTD algorithm three-fold, such that the resulting sub-cells are still square/cube-shaped, but (i) they are typically overlapping, (ii) have special labels that contain instructions for computing the integrals over the individual and 7, respectively. The actual meaning of the labels denoted by L and the color-coding associated with them will be discussed in the following sub-sections. From Fig. 5 vs. Figure 6 and Fig. 7a vs. Figure 7b, the key idea and the potential of the Boolean approach for reducing the number of subcells n SC and integration points n IP is already visible. In fact, in Sect. 3.3 and 6, it will be shown that the much more efficient distribution of integration points via the B-QTD/B-OTD algorithm leads to reduction in computational time by approximately 70-80%, while maintaining the same accuracy.

Sub-cell labels
The Boolean approach sketched in Eq. (10) is based on integrating a discontinuous function in two steps: (i) a continuous function is integrated over a simple square/cube-shaped domain and (ii) a discontinuous function is integrated over a specific sub-domain, where the integrand results from the subtraction of two known functions. Although in the B-QTD/B-OTD, the integral value is not computed in two steps like in case of Eq. (10), but in n sc steps each sub-cell is directly connected to one of the two terms in Eq. (10). Thus, although multiple calculations are required, in essence, the same procedure as described before applies. The actual purpose of a given sub-cell and the form of the discon-tinuous function D in the integrand for that specific sub-cell, is explicitly encoded in the developed labeling system. In particular, two different label types have been developed, which are comprehensively discussed in the following paragraphs.

Label type 1:
The label of type 1 is composed by two domain IDs i and j in a specific order, and refers to the following integral over the given sub-cell The above integral essentially corresponds to the first term on the right hand side in Eq. (10), where the whole domain of interest is integrated and the discontinuity is ignored, regardless of whether the sub-cell is cut or not. However, instead of just using a single value of D i like in Eq. (10), the term D i − D j is inserted, such that an overlapping sub-cell structure can be exploited. Note that label type 1 is always assigned to sub-cells created at the levels R < R max . Additionally, we would like to express integrals, where only D i or −D j is used. Therefore, an additional zero value D 0 := 0 is declared, while keeping in mind that the actual domain IDs start from 1, not from 0, cf. Eq. (2). Thus, using i = 0 or j = 0, special cases of Eq. (20) are also possible, which are needed for the algorithm to work Label type 2: The second label type contains only a single domain ID j, and refers to the integral where, unlike for type 1 labels, the discontinuity is taken into account due to the presence of D. The above integral corresponds to the second term on the right hand side in Eq. (10). Here, the discontinuity is taken into account and hence, the integration error resulting from the label type 1 sub-cells is corrected. Note that for points lying in j , the integrand is vanishing due to Eq. (9) This fact constitutes a key feature of the B-QTD/B-OTD algorithm, since every sub-cell and integration point located in j can be directly discarded. Similar to Eqs. (21) and (22), if j = 0 is used, Eq. (23) turns into simply integrating the discontinuous function D f over ω B due to the applied definition D 0 = 0.

Labelling of new sub-cells
The labeling procedure is illustrated in Fig. 8 for the partitioning step of a quadratic sub-cell ω k ∈ R 2 with the type 2 label L[ j]. In principle, there are three cases how the B-QTD can proceed from this stage, which are elaborated below. Note that at the very beginning of the Boolean decomposition algorithm, j = 0 is set for the label.

Case 1:
If the given sub-cell is not cut, it must completely belong to a single domain ω k ⊂ i . The sub-cell is excluded from further partitioning and its label modifies to L[i, j] (orange box in Fig. 8).
If the sub-cell is cut, four children are created and the following two additional steps are executed for the Boolean labeling process: (i) Find the material ID i for the sub-cell ω B k , for which k i = ω B k ∩ i is the largest. In the current implementation, this is achieved by testing a set of sample points. (ii) Test whether k i occupies at least 1 /4 th of ω B k in 2D and 1 /8 th in 3D. This condition is denoted by C in Fig. 8. Case 2: If C = false, the parent sub-cell is deleted and all of its children inherit its label (green box in Fig. 8). Essentially, this scenario corresponds to the general procedure of the standard QTD, with the difference, that the sub-cells are also labeled. Since C is violated, even the largest sub-domain has an area less then 25% of that of the parent sub-cell. Thus, none of the children can belong to a single sub-domain only. Note that C = false only applies in cases where the subcells have at least 5 sub-domains. The benefit of the labeling procedure in this case is demonstrated in Fig. 9, where the example domain is discretized by a single cell with R = 1. Since C = true for 2 , the parent cell is kept and the labels are assigned according to Fig. 9b, where the white text boxes are used for numbering the subcells. In Fig. 9c, the vanishing integrand property given in Eq. (24) is utilized: Firstly, since sub-cells No. 3 and 4 lie completely in 2 , these can be directly discarded. Secondly, while sub-cells No. 2 and 5 are kept, a further reduction of integration points is still possible due to the integrand beeing D 2 − D 2 in ω B 2 ∩ 2 and ω B 5 ∩ 2 (green domain in Fig. 9c). Both, discarding the entire sub-cells and the integration points are possible due to the profound interplay between the assigned labels.

Remark 1
Note that performing a standard QTD and then assigning labels to the sub-cells is not possible. The labeling process has to be carried out parallel to the QTD.
The procedure depicted in Fig. 8 is repeated for each cut sub-cell until the final refinement level R max is reached. Every time C = true, both a parent sub-cell and some of its children are kept. Using R max = 3 for Fig. 9a results in the set of overlapping sub-cells depicted in Fig. 6, where the Raxis indicates the refinement levels the individual sub-cells were generated on. A special color-coding is used to indicate the different sub-cell labels, whose corresponding integrands are depicted in the right side of Fig. 6.
Examining an arbitrary point P = [−0.1, −0.35] ∈ 2 in Fig. 6, the integrand at P should evaluate to I = D 2 (cf. Fig. 4a). By adding up the integrands represented by the Boolean sub-cells containing the given point, (a) Exemplary domain.
it is evident, that the equality indeed holds. Note that on R = 3, the type 2 label is used as defined in Eq. (23).
Finally, the effect of the labeling concept on the integration points is depicted in Fig. 10, where 3 × 3 integration points are distributed in each sub-cell of Fig. 6 and the same coloring scheme is used. Additionally, for each integration point, the corresponding integrand is given. The difference between the different label types is clear: Since sub-cell labels at R < R max are not taking the discontinuity into account (red and blue dashed lines), all integration points are used in the sub-cells. For the sub-cells at R = R max , label type 2 is used, where the discontinuity (red and blue solid lines) is taken into account during the integration. These cells generally contain integration points with vanishing integrand, such that they can be excluded from the integration procedure due to Eq. (24).

Performance
In this sub-section, the performance of the B-QTD and B-OTD algorithms is investigated in terms of the reduction of integration points and achieved integration accuracy. Due to the Boolean sub-cell labels, the B-QTD has a much more efficient distribution of the integration points than the standard QTD. This can be clearly recognized in Fig. 7 presented earlier, where the integration points for the B-QTD are based on = 81 for the current example, when the QTD is extended by the Boolean approach. Using these quantities, the reduction of integration points is introduced as where r IP > 0 % indicates a meaningful reduction and r IP = 0 % means no reduction. For the example given in Fig. 7, a significant saving by r IP = 67.86% can be obtained. Similar to r IP , a reduction rate of sub-cells r SC can also be defined.
Both of these quantities are depicted in Fig. 11a for various refinement levels, where a reduction rate up to r IP ≈ 80% is observed. Note that due to the additional savings of integration points in the leaf sub-cells (cf. Fig. 10d), r IP ≥ r SC .
Furthermore, let e I express the relative integration error where I is the integral value obtained by different numerical methods and I ref the reference value. For the current example I ref is given in Eq. (15). Figure. 11b reveals that despite the significant reduction of integration points, the B-QTD yields the same integration accuracy as the QTD. The same significant reduction rate and no loss of accuracy can be observed for 3D problems as well, as depicted in Fig. 12 for a cubic domain composed by the sub-domains 1 , 2 , 3 , defined below.

Fig. 10
Integration points in the overlapping sub-cells of Fig. 6 and their corresponding integrands

Multi-material finite cell method
In this section, the discussion in Sect. 1.1 is extended by the required mathematical formulation for multi-material problems. Let us consider a material domain that is composed by n d disjoint sub-domains = ∪ n d i=1 i , which are either void regions or other material domains.

Finite cell method
The original FCM formulation is concerned with porous media, i.e., the given problem consists of a single material and multiple void regions [3][4][5]43]. It has been shown, that exponential convergence rates can be obtained [10]. Without derivation, the weak formulation of linear elasticity in the context of the FCM reads Here, u is the displacement field, v denotes the test function, and S h together with V h represent the appropriate finite dimensional function spaces, which in this article, are based on the spectral element concept [54][55][56]. Furthermore, B e and F e are bilinear and linear functionals, respectively, defined over the extended domain e In the above equations, α is the indicator function defined in Eq. (1) for distinguishing between the physical and fictitious domains, while ρ and C are the material density and stiffness, respectively. Furthermore,ü is the acceleration field, b stands for the body load vector, andt represents the prescribed tractions along the Neumann boundary N . Considering a two-dimensional setting for the sake of simplicity, in each finite cell, the displacement field 2 ] T is approximated using smooth shape functions where N u ∈ R 2×2n N is a matrix containing the shape functions associated with the n N nodes of the given cell and U (c) contains the unknown nodal displacements. The cell matrices as well as the global system matrices of the

Local enrichment
In the following, the focus is kept on multi-material problems without any void regions. Thus, α = 1 applies to all cases. Nonetheless, the material properties are still discontinuous due to C(x) = C i (x) and ρ(x) = ρ i (x) ∀x ∈ i . Furthermore, at the material interfaces, the displacements are generally C 0 -continuous while stress and strains obey certain jump conditions. For locally enriched problems, the weak formulation reads where u is composed of the smooth u u and enriched displacements u a , such that u = u u + u a . Furthermore, v u and v a are the test functions associated with the base V h u and enriched function spaces V h a , respectively. Following from Eq. (38), for the approximation of u (c) , Eq. (36) turns into where the first term is identical to the standard approximation in Eq. (36) and the second term includes the matrix of enrichment shape functions and enrichement DOFs A. There are multiple ways how the enrichment function ψ in Eq. (40) can be defined. In this contribution, we follow the modified-abs enrichment proposed by Moës et al. [14] ψ = N|ϕ| − |Nϕ|, which is based on a level-set function ϕ i of the embedded geometries, where the domain and boundary of the ith geometric entity is defined by In Eq. (41), N and ϕ are vectors containing the shape functions and level-set values associated with the individual nodes, respectively Furthermore, the last term in Eq. (41) ensures the C 0continuity of ψ along the interface, and subtracted from the first term, it yields a ψ that vanishes on the element edges that are not intersected by the boundary. Thus, parasitic terms and partially enriched transition elements can be avoided [57]. Note that in the general case, N u , N a and N can be based on different Ansatz spaces and meshes [6]. In this contribution, for all of these quantities, Lagrangian shape functions are used based on the Gauss-Legendre-Lobatto (GLL) nodal distribution. While for N u , a user defined polynomial order p is used, the polynomial degree p ψ of N depends on the complexity of intersecting interface [6]. Finally, following from Eq. (40), the polynomial degree of N a is p a = p + p ψ .

Discretization
Discretization of the coupled weak form leads to cell-specific stiffness matrices including terms that are only related to (i) the smooth displacements K c uu , (ii) the enriched displacements K c aa , and (iii) appropriate coupling terms K c ua = (K c au ) T . Following the QTD/OTD-based integration scheme often used in the FCM (see Sect. 3.1), K c is computed by integrating over the n sc sub-cells In the above equation, B u and B a are the strain-gradient operators associated with the standard and enriched displacements, respectively. Similar to B u , B a is constructed as where B (l) a is associated with the lth enrichment shape function, and for two-dimensional linear elastic problems it is of the form The differentiation of the individual entries w.r.t. x m is realized using the product rule with m = 1, 2. Based on Eq. (41), the derivative of the enrichment term ψ is Note that here the chain rule is utilized again and that the sgn-function results from differentiating the abs-function. Similar to Eq. (46), the cell-specific mass matrix has the structure Due to the enrichment, the highest polynomial degree in the integrands of K c and M c is 2 p a = 2( p + p ψ ). Thus, for an accurate numerical integration 2 , instead of the standard p + 1 integration points per direction, p + p ψ + 1 integration points are required. Finally, after assembly, the coupled global equation system without damping reads

B-FCM for multi-material problems
In this section, the Boolean integration approach given in Sect. 2 is combined with the local enrichment of the displacement field over cut cell. For simplicity, the focus is kept on cases where cut cells are intersected by nothing else but a single material interface. Thus, the cells contain two material domains 1 and 2 , without the presence of any void regions. However, more complex scenarios are possible and constitute a straightforward extension of the proposed approach.

Intuitive implementation of the cell matrices
In order to apply the features presented in Sect. 3.2, a clear separation of the discontinuous (D) and continuous terms ( f ) is required. In the case of local enrichment, the discontinuity in Eq. (46) is not only caused by the strongly discontinuous material properties related to C, but also by B a , that is derived from the weakly discontinuous enrichment function ψ. Thus, when integrating the enriched stiffness matrix over a given Boolean sub-cell ω B k , the integration labels defined in Eqs. (20) and (23) lead to and Note that in the above equations, the terms are separated such that all weakly/strongly discontinuous terms are within the brackets. As a next step, the extension of these terms to the entire integration domain is required. On the one hand, this is not an issue for C 1 and C 2 , which are often piecewise constant functions. On the other hand, B a1 and B a2 are polynomial functions and form a piece-wise polynomial B a , such that Since the discontinuity in B a is caused by the enrichment function ψ, the computation of ψ i and its spatial derivatives are required for calculating B ai . Based on Eqs. (48) and (49), for the sub-matrix of B ai associated with the lth enriched shape function, When using the modified abs-enrichment approach to construct ψ, (cf. Eq. (41)), the weak discontinuity is introduced by the term |Nϕ|. Thus, the enrichment function is divided along the isoline Nϕ = 0 into ψ 1 and ψ 2 , such that Since |Nϕ| = Nϕ for Nϕ > 0 and |Nϕ| = −Nϕ for Nϕ < 0, ψ 1 and ψ 2 read An example for this is given in Fig. 13, where Fig. 13a depicts a specific ψ whose parts over k 1 and k 2 are color-coded by green and orange. The extensions of ψ 1  Finally, since in Eqs. (58) and (59), the vectors N and ϕ contain polynomial and constant entries, respectively, the required partial derivatives of ψ 1 and ψ 2 in Eq. (56) are computed as The formulation presented in this section sketches an intuitive concept of formulating the integrands over the Boolean subcells. Although it utilizes a reduced set of integration points, Eqs. (53) and (54) contain additional matrix operations. An example for this can be seen in the computation of K c aa by means of Eq. (53), where instead of B T a C B a , B T ai C i B ai − B T a j C j B a j is computed, i.e., the number of matrix operations is increased by a factor of two. For avoiding this feature, the next sub-section presents a more compact formulation, reducing the number of unnecessary matrix operations.

Compact implementation of the cell matrices
In the current sub-section, the notion of Sect. 5.1 is extended by an improved separation of the continuous and discontinuous terms, such that unnecessary or redundant matrix operations can be avoided and the computational time of the Boolean approach is further reduced. More details concerning the reduction in computational effort are given at the end of this section.

More efficient separation of terms
As a first step, Eq. (48) is rewritten, such that where is a 3 × 2 matrix containing the spacial derivatives of ψ. The same approach also applies when computing the entire matrix Note that in this case B a is reproduced by using the terms N u and B u , which are generated only once for a given sub-cell. Using Eq. (63), let us rewrite the purely enriched term K c aa in Eq. (46): The above equation can be also expressed in a matrix form as where D aa ∈ R 5×5 contains all discontinuous properties associated with K c aa in a compact form. Since C is symmetric, (ψ C ) T = ψ T C, and the matrix D aa is also symmetric Using the same idea, the following two equations redefine the mixed terms associated with K c ua and K c au : Finally, the definition is also applied to K c uu , where due to C being the only discontinuous term, D uu = C applies Collecting all matrices containing the discontinuous terms in Eqs. (67)-(70), the set D is defined, which is, in the remainder of this section, associated with the enriched K c of Eq. (46):

Boolean operations
Let us now combine the formulation from the previous subsection with the Boolean approach. All discontinuous terms in D are composed by individually continuous parts, such that where The above expression only requires evaluating the matrices C i and the scalars ψ i and ∂ψ i /∂ x m , for i = 1, 2 and m = 1, 2, where the latter ones constitute i . Using this formulation, Eqs. (53) and (54) can be reformulated as and The application of the same approach to the enriched mass matrix M c and body load vector F c is presented in Appendices A.2 and A.3, respectively.
Recall that computing K c aa in Eq. (53) involves evaluating the matrices B a1 and B a2 , which are in turn used to determine the terms B T a1 C 1 B a1 and B T a2 C 2 B a2 as well as their difference (cf. end of Sect. 5.1). In the more compact version presented by Eq. (74), the matrices B a1 and B a2 are not computed explicitly and the Boolean operation is performed on significantly smaller matrices using D aa 1 and D aa 2 . Thus, while the intuitive approach is more straightforward in its implementation, the more compact approach presented in this section generally leads to a decreased numerical effort.

Circular plate with inclusion
In this section, a benchmark problem with a circular inclusion 1 embedded into a circular plate 2 is considered [58,59]. Along the outer boundary of the plate, radial displacement boundary conditions are applied where a and b are the radii of the inclusion and embedding plate, respectively. The parameter β is determined by a and b, together with the Lamé constants of the inclusion {r 1 , μ 1 } and the matrix {r 2 , μ 2 } Instead of simulating the entire domain, a quadratic domain with a side length of c is considered with appropriate non-homogeneous Dirichlet boundary conditions (Fig. 14a).
Regarding the geometry and material properties, the following parameters are used: a = 3 mm, c = 8 mm, b = 15 mm, E 1 = 10,000 MPa, ν 1 = 0.3, E 2 = 0.1E 1 and ν 2 = 0.27. Furthermore, a plane strain state is assumed. The analytical reference strain energy for the simulated domain is For simplicity, the quadratic domain is discretized by a single finite cell and spectral shape functions with polynomial degrees from p = 1 to p = 8 are tested. For ψ in Eq. (41), quadratic shape functions are used, which are sufficient for approximating the level-set function of the circle. For each step of the p-refinement, the QTD is performed with R = p + 2 and ( p + 10) 2 integration points per sub-cell are used 3 .  When using the B-FCM approach, the simulation accuracy is maintained when compared to the standard FCM, resulting in basically identical convergence of the relative error in the energy norm as depicted in Fig. 15a. Here, the B-FCM versions 1 and 2 refer to the different versions of the method introduced in Sects. 5.1 and 5.2, respectively. Additionally, Fig. 15a also depicts the severely deteriorated convergence curve of the case when no enrichment is used. Furthermore, the accuracy of the simulation is also demonstrated in Fig. 15b, where local values are evaluated for p = 8 along the cut line A-A', indicated in Fig 14b. The red curve depicts the radial displacements, in which the sharp kink at the material interface is well captured, while the blue curve depicts the non-oscillatory strain energy density ψ = 1/2 σ : ε. Thus, it can be stated that the B-FCM achieves very accurate results, while requiring significantly less integration points and computational time compared to the standard FCM. Figure 16a compares the wall-clock times spent on the cut cells in an absolute, while Fig. 16b in a relative manner Additionally, the reduction of integration points is also depicted, which is fairly constant for the given example within a typical refinement range of R = 3...10: r IP = 70−76%. The reduction of computational time r t starts at about 30%, however with increasing p and R, it reaches values of r t ≈ 60% for the B-FCM 1st version and r t ≈ 76% for the B-FCM 2nd version. The reason for the larger difference between r IP and r t is due to the fact, that for low values of p and R, the integration time of K c accounts for a relatively smaller chunk of the overall computational time spent on the cut cell. For higher values of p and R, the numerical integration is more costly, and the effect of reducing the integration points is more prominent. The additional time reduction by the blue curve is due to the more efficient formulation of the matrix operations in case of the 2nd version when compared to the 1st one. Note however, that despite these differences, both versions of the B-FCM (i) operate on the same number of integration points, (ii) reduce the computational time by a significant amount, while (iii) yielding the same accuracy as the standard FCM.

Cube with spherical inclusion
Next, a 3D test is conducted involving a cube ( 2 ) of side length a = 4 m and a spherical inclusion ( 1 ) with a radius of R = 2 m, as depicted in Fig. 17a. The problem is formulated using the method of manufactured solutions (MoMS) [60,61], i.e., a displacement field u * is given, for which the strong form of linear elastostatics 4 yields the corresponding body loads b = −div(σ ). These body loads, together with the appropriate boundary conditions are then given as input for the FCM simulation and the resulting u is compared to u * . Using this approach, the code functionality can be tested based on an analytical reference solution, and no overkill FEM solution is required. For the given example, a radial displacement field u * = [u * r , u * ϕ , u * θ ] is chosen in spherical coordinates    where Note that for u * r (r = R) = 0 and for any c = 1, the manufactured displacement field exhibits a kink at the interface ( Fig. 18b). In the current example, c = 0.01 is chosen and the deformed shape of the simulation domain is depicted in Fig. 17b. For the material properties, E 1 = 0.1 Pa, E 2 = 10 Pa and ν 1 = ν 2 = 0.3 are chosen. The body loads,  expressed in a spherical coordinate system, read for the inclusion and matrix, respectively. Due to our specific choice of c, E 1 = cE 2 , and thus, b r1 = b r2 . Consequently, tractions are continuous along the interface and jump conditions are naturally fulfilled (see Fig. 18d). In this special case, when following the MoMS approach, no weak boundary conditions along the interface are required and the problem complexity is kept at the necessary level. 5 The body load vectors for the B-FCM are given in Eqs. (100) and (101). 5 It is noted that in a forthcoming publication, the application of MoMS to embedded domain methods is studied in detail. To this end, special considerations needed in a embedded domain framework are comprehensively discussed and algorithms to make MoMS available for verification purposes are presented.
Finally, the analytically evaluated reference strain energy for the manufactured solution reads For the current example, the simulation domain is discretized by a single finite cell only and the polynomial degrees p = 1...6 are tested. Furthermore, a constant R = 8 is chosen for all polynomial degrees and ( p + 10) 3 integration points are used per sub-cell. The global simulation accuracy measured by the relative error in the energy norm, is depicted Fig. 18a, where similar to the 2D case, the B-FCM approach yields the same convergence curve for the current 3D problem as the FCM. Furthermore, for the B-FCM simulation with p = 6, Fig. 18b-d depict local quantities of the solution, such as the radial displacement, strain and stress fields along the cut line B-B' indicated in Fig. 17a. For all cases, the analytical reference solutions are also shown, demonstrating the high accuracy of B-FCM, while using a single finite cell only. Analogous to Fig. 16 of the previous example, Fig. 19 depicts the absolute computational times as well as the reduc-  tion of those (r t ) and of the integration points (r IP ). Due to the constant refinement level, the reduction of integration points is basically constant, r IP ≈ 80%. However, since (i) generating K c in 3D is computationally more intensive, and additionally, (ii) F c body also has to be computed, the difference between the presented B-FCM versions is more pronounced. For the 2nd version (blue curve), r t is about 70% for all values of p, while for the 1st version (red curve), the unnecessary matrix operations in K c lead to additional computational overhead for higher values of p, rendering the method less effective. Nonetheless, even in this case, 50-70% of the computational time is saved.

Acoustic meta-material
Finally, the proposed multi-material B-FCM approach is applied to the analysis of an acoustic meta-material placed on a metal plate, where the meta-material consists of a foam matrix and several inclusions (Fig. 20). The simulation of the surrounding acoustic field is omitted at this point for simplicity, and the focus is kept on the vibrating structure. For testing the proposed integration scheme, a harmonic analysis is conducted on the structure using the material properties given in Table 1.
The simulation domain is discretized by 5×30 finite cells, among which 28 are cut by the embedded material interfaces (yellow cells). The plate is subjected to a harmonic excitation wherep is the pressure amplitude, i the imaginary unit and = 2π f the angular excitation frequency. The system answerÛ is computed in the frequency domain by where a stiffness proportional damping is used Note that since the damping factor η is also yet another discontinuous material property (Table 1), Eq. (90), cannot be applied directly to generate the global stiffness matrix K . Instead, it is realized on the sub-cell level where K c k,i is the contribution to the stiffness matrix K c over k i ⊂ ω B k , i.e., the ith material domain in the kth subcell, and η i is the damping factor in that material region. Furthermore, n k d is the number of sub-domains in ω B k . Note that Eq. (91) holds both for the FCM and B-FCM approaches, i.e., regardless of how K c is computed. Table 2 compares the FCM and B-FCM approaches in terms of the number of sub-cells (SC) and integration points (IP) generated in the cut cells, as well as the computation times required for creating the local integration mesh (LIM), the cut cell procedure, assembly, and solver. In both cases, p = 5 and R = 5 is used with ( p +10) 2 integration points in each sub-cell. In the last row, the reductions of the different values are computed. Similar to the previous examples, a high reduction of integration points is obtained, r IP = 74.10%. While creating the LIM requires additional time in the B-FCM, its cost is negligible to the time spent on numerical integration. Thus, the overall time spent on the cut cells, that includes other miscellaneous tasks besides creating the LIM and integration, is reduced by 65%. Since the cut cells represent the computationally most intensive part of the assembly procedure, the assembly time is reduced by a similar amount. The deformed shape of the vibrating meta-material computed with the B-FCM is depicted in Fig. 21.

Conclusion
In this contribution, the B-FCM formulation for porous media (Sect. 1.3) is extended to multi-material problems for efficient simulation of complex structures with material interfaces in the FCM framework. The numerical integration in the proposed B-FCM approach is based on a quadtree-/octreedecomposition (QTD/OTD) of the cut cells (Sect. 3.1), as it is often done in the FCM. However, in the B-FCM implementation, sub-cells are equipped with special labels (Sect. 3.2), whose purpose is twofold: 1. They steer the decomposition, such that it generally results in an overlapping structure of integration domains (Figs. 6 and 10), leading to significantly less sub-cells and integration points. 2. The labels hold specific information regarding the form of the integrand over the individual sub-cells, such that Boolean operations can be used to successively compute the integral value over the overlapping sub-cells.
In our implementation, the C 0 -continuous displacement field along the material interfaces is captured by a local enrichment approach. The multi-material B-FCM requires certain changes to the integrand such that the Boolean nature of the sub-cells can be exploited. While these changes are straightforward for discontinuities in the material properties, further steps are needed for discontinuities caused by the chosen enrichment function, as derived in Eqs. (58)- (61). Numerical examples are conducted both in 2D and 3D, and the multi-material B-FCM concept is applied to the computation of stiffness and mass matrices as well as to the body load vector. The results show that compared to the conventional QTD/OTD-based FCM, its Boolean version requires up to 80% less integration points for reasonable refinement levels. Furthermore, despite a significant reduction in the computational effort, the global and local accuracy of the simulation is not altered. These statements apply to both presented B-FCM formulations in this paper. While the 1st version enables a more intuitive understanding of the integration over the Boolean sub-cells (Sect. 5.1), the 2nd version (Sect. 5.2) is based on a more efficient separation of continuous and discontinuous terms and thus, reduces the number of redundant matrix multiplications. For simpler problems, the difference in the two approaches is less prominent, as they reduce the time spent on the cut cells by ∼55-65% and ∼55-75%, respectively (Fig. 16). However, for problems with larger cell matrices, e.g., in case of 3D examples or when using high polynomial orders, the 2nd version is superior, leading to ∼70% time reduction, while the 1st version achieves roughly 50%.
In conclusion, while the implementation of the B-FCM requires slight modification of the decomposition scheme and integrands over the sub-cells, it reduces the computational effort of the QTD/OTD-based integration schemes by a significant amount, while enjoying the same level of robustness in combination with no loss in accuracy. If the simplest implementation of the B-FCM is required (e.g., for testing purposes), we recommend its 1st version, and if the maximum efficiency is desired, its 2nd version should be exploited.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A.2 Mass matrix
With the relation for N a given in Eq. (40), the enriched mass matrix over the sub-cell ω B k in the cell c is given by Using the proposed integration approach, M c over the Boolean sub-cells is computed by for the different label types. In the brackets, the discontinuity is taken into account, which in case of bi-material problems, requires the evaluation of the scalar terms ρ 1 , ρ 2 , ψ 1 , ψ 2 , ψ 2 1 and ψ 2 2 . Note that the product N T u N u must be computed only once for the given sub-cell.

A.3 Body load vector
The enriched body load vector over the sub-cell ω B k is computed by Using Boolean sub-cells for the numerical integration and the relation for N a given in Eq. (40), F c body over ω B k reads for the two different label types