Error-estimate-based Adaptive Integration For Immersed Isogeometric Analysis

The Finite Cell Method (FCM) together with Isogeometric analysis (IGA) has been applied successfully in various problems in solid mechanics, in image-based analysis, fluid-structure interaction and in many other applications. A challenging aspect of the isogeometric finite cell method is the integration of cut cells. In particular in three-dimensional simulations the computational effort associated with integration can be the critical component of a simulation. A myriad of integration strategies has been proposed over the past years to ameliorate the difficulties associated with integration, but a general optimal integration framework that suits a broad class of engineering problems is not yet available. In this contribution we provide a thorough investigation of the accuracy and computational effort of the octree integration scheme. We quantify the contribution of the integration error using the theoretical basis provided by Strang's first lemma. Based on this study we propose an error-estimate-based adaptive integration procedure for immersed isogeometric analysis. Additionally, we present a detailed numerical investigation of the proposed optimal integration algorithm and its application to immersed isogeometric analysis using two- and three-dimensional linear elasticity problems.


Introduction
Immersed finite element methods -such as, e.g., the finite cell method (FCM) [1], CutFEM [2] and immersogeometric analysis [3][4][5] -have been demonstrated to be suitable for computational problems for which the performance of mesh-fitting finite element methods is impeded by complications in the meshing procedure. In recent years, immersed finite element methods have been successfully combined with isogeometric analysis (IGA) [6], a spline-based higherorder finite element framework targeting the integration of finite element analysis (FEA) and computer aided design (CAD). Immersed methods provide the possibility to conduct IGA on volumetric domains based on a CAD boundary surface representations [5] or voxelized geometry data [7]. Moreover, immersed methods provide a natural framework for the consideration of trimming curves in IGA [4,[8][9][10][11][12][13][14][15]. Immersed IGA has been applied successfully to various problems in, amongst others, structural and solid mechanics [16,17], fluid-structure interactions [18][19][20] and scan-based analysis [7,21]. In comparison to mesh-fitting finite element methods, immersed methods require special treatment of various computational aspects. A prominent computational challenge that is inherent to immersed finite element methods is the integration over cut-cells, a problem that closely relates to the special treatment of discontinuous integrands in enriched finite element methods such as XFEM and GFEM. Since the geometry of the computational problem is captured by the integration procedure rather than by the ambient mesh in which the domain is immersed, cut-cell integration techniques must be capable of adequately capturing a wide range of cut-cell configurations. A myriad of dedicated integration procedures with this capability has been developed over the years in the context of immersed FEM (see [16,22] for reviews/comparisons) and enriched FEM (see [23]), which can be categorized as: • Octree subdivision: The general idea of octree (or quadtree in 2D) integration is to capture the geometry of a cut-cell by recursively bisecting sub-cells that intersect with the boundary of the domain, as illustrated in Fig. 1. At every level of this recursion, sub-cells that are completely inside the domain are preserved, while sub-cells that are completely outside of the domain are discarded. This cut-cell subdivision strategy was initially proposed in the context of the finite cell method in Ref. [24] and is generally appraised for its simplicity and robustness with respect to cut-cell configurations. Octree integration has been widely adopted in immersed FEM, see, e.g., Refs. [7,[16][17][18]. Various generalizations and improvements to the original octree procedure have been proposed, of which the consideration of tetrahedral cells [25,26], the reconstruction of the immersed boundary by tessellation of the lowest level of bisectioning [7], and the consideration of variable integration schemes for the sub-cells [27], are particularly noteworthy. Despite the various improvements to the original octree strategy, a downside of the technique remains the number of integration sub-cells that results from the procedure, especially in three-dimensional cases.
• Cut-cell reparametrization: Accurate cut-cell integration schemes can be obtained by modifying the geometry parametrization of cut-cells in such a way that the immersed boundary is fitted. This strategy was original developed in the context of XFEM by decomposing cut elements into various sub-cells containing cut-cells with only one curved side and then to alter the geometry mapping related to the curved sub-cell to obtain a higher-order accurate integration scheme [28]. This concept has also been considered in the context of implicitly defined geometries (level sets) [29], the NURBS-enhanced finite element method (NEFEM) [30,31] and the Cartesian grid finite element method (cgFEM) [32]. In the context of the finite cell method the idea of cut element reparametrization has been adopted as part of the smart octree integration strategy [33,34], where a boundary fitting procedure is employed at the lowest level of octree bisectioning in order to obtain higher-order integration schemes for cut-cells with curved boundaries. Reparametrization procedures have the potential to yield accurate integration schemes at a significantly lower computational cost than octree procedures, but generally compromise in terms of robustness with respect to cut-cell configurations.
• Polynomial integration: Provided that one can accurately evaluate integrals over cut-cells (for example using octree integration), it is possible to construct computationally efficient integration rules for specific classes of integrands. In the context of immersed finite element methods it is of particular interest to derive efficient cut-cell integration rules for polynomial functions. The two most prominent methods to integrate polynomial functions over cut-cells are moment fitting techniques [34][35][36][37][38], in which integration point weights and (possibly) positions are determined in order to yield exact quadrature rules, and equivalent polynomial methods [39,40], in which a non-polynomial (e.g., discontinuous) integrand is represented by an equivalent polynomial which can then be treated using standard integration procedures. Such methods have been demonstrated to yield efficient quadrature rules for a range of scenarios. A downside of such techniques is the need for the evaluation of the exact integrals (using an adequate cut-cell integration procedure) in order to determine the optimized integration rules. This can make the construction of such quadrature rules computationally expensive, which makes that they are of particular interest mainly in the context of time-dependent and non-linear problems, for which the construction of the integration rule is only considered as a pre-processing operation (for each cut-cell) and the optimized integration rule can then be used throughout the simulation.
• Boundary integral representation: Depending on the problem under consideration, it can be possible to reformulate volumetric integrals over cut-cells by equivalent boundary integrals. This reformulation, which has been proposed in the context of XFEM in Ref. [41] and in the immersed FEM setting in Ref. [42], is advantageous from a computational effort point of view, as the equivalent boundary integrals are generally less costly to evaluate. Moreover, provided that an adequate description of the boundary surface is available, higher-order accurate integration evaluations are obtained. A downside of integration techniques of this kind is that they require reformulation of the volume integrals, which makes them less general than standard quadrature rules.
In the selection of an appropriate cut-cell integration scheme one balances integration between robustness, accuracy and expense with respect to cut-cell configurations. If one requires a method that automatically treats a wide range of cut-cell configurations and is willing to pay the price in terms of (higher-order) accuracy and computational effort, octree integration is the compelling option. If constraints are imposed from an accuracy and computational expense point of view and one has some control over the range of configurations to be considered, alternative techniques such as cut-cell reparametrization are attractive. The need to balance between robustness, accuracy and expense has driven the development of hybrid integration schemes, such as smart octree integration [33] and adaptive moment-fitting [37], which allow one to attain an integration procedure with the desired properties.
Balancing accuracy with computational effort does not necessarily require the consideration of hybrid integration procedures, but can also be achieved by controlling the parameters of the integration procedures (listed above). This is particularly the case for octree integration, the accuracy of which can be controlled by the bisectioning depth and the integration orders used on the sub-cells. Optimization of these parameters to reduce the computational expense of octree integration without compromising its robustness with respect to configurations was proposed in Ref. [27], where an algorithm is proposed to select the integration order on the different levels of sub-cells. It is demonstrated that reducing the integration order with increasing bisection depth significantly reduces the number of integration points, without unacceptably compromising the accuracy.
The idea of Ref. [27] to reduce integration orders on certain levels of the octree subdivision is in agreement with the theory of finite elements. Strang's first lemma [43][44][45] provides a framework to incorporate integration effects in the error analysis of finite element methods. This lemma indicates that integration does not need to be exact in order to attain (optimal) convergence [43]. 1 This lemma has been considered for the analysis of CutFEM, see, e.g., Refs. [46,47]. In fact, the idea of reduced integration in finite element methods has been studied extensively over the last decades, with applications in the analysis of plates and shells [48] and mixed finite element methods [49,50] being prominent examples.
In this manuscript we propose an error-estimation-based adaptive algorithm to optimize the parameters of an octree integration scheme for cut-cells. Although the optimization procedure is here tailored toward the octree integration procedure, the presented theoretical framework is applicable to a range of cut-cell integration techniques. When applied to the octree integration procedure, the presented framework provides the capability to evaluate sub-cell integration errors in accordance with Strang's first lemma [43][44][45]. When the sub-cell errors are combined with a computationalcost estimate based on the number of integration cells, a refinement indicator that optimally balances computational effort and accuracy is obtained. The proposed algorithm deviates from the one proposed in Ref. [27] in that it is directly based on integration error evaluations. The effectiveness of the adaptive integration procedure is demonstrated in the context of the isogeometric finite-cell framework. Besides the ability to obtain optimal integration rules for cut-cells, the developed adaptive integration procedure provides a tool to postulate rules of thumb for the selection of integration orders at different levels of bisectioning, and to assess the quality of integration rules derived using alternative algorithms such as that proposed in Ref. [27,51].
The paper outline is as follows. In Section 2 we introduce the finite cell method, with a focus on the application of Strang's lemma to estimate integration errors. In Section 3 we study the influence of the octree integration scheme on the computational cost, and we define and evaluate the integration error of a sub-cell in a cut element. Based on the evaluated integration error, we present an optimization algorithm. The optimal integration procedure is investigated for two-and three-dimensional numerical test cases with various cut element configurations in Section 4. The developed algorithm is then applied to a two-and three-dimensional immersed isogeometric analysis problem in Section 5. Finally, concluding remarks are presented in Section 6.

The finite cell method
Because the integration procedure proposed in this work applies to a wide range of finite-cell simulations, in Section 2.1 we first introduce the finite cell method in abstract form. Based on this abstract problem setting, Section 2.2 presents an error analysis that incorporates integration errors. In Section 2.3 the Poisson problem is considered to exemplify the abstract derivations.

Formulation
We consider a domain Ω ⊂ R d (d ∈ {2, 3}) with boundary ∂Ω as shown in Fig. 2. The boundary is split into a part on which Dirichlet conditions are imposed, denoted by ∂Ω D , and a complementary part to which Neumann conditions apply, denoted by ∂Ω N . We consider a problem described by a field variable u. Let W ∋ u denote a suitable ambient space for the solution, equipped with a Hilbert structure corresponding to the inner product (·, ·) W and the association norm ∥ · ∥ W . Similarly, V is a Hilbert space with inner product (·, ·) V and norm ∥ · ∥ V , which encompasses the test space for the weak formulation of the problem under consideration. To accommodate the Dirichlet boundary conditions, let the spaces W 0 ⊂ W and V 0 ⊂ V be composed of functions that vanish at the Dirichlet boundary in a suitable manner. Moreover, let ℓ g ∈ W denote a lift of the Dirichlet data. We consider a weak formulation of the generic form: We assume that the bilinear form a : W × V → R is continuous on W × V and weakly coercive on W 0 × V 0 , and the linear form b : V → R is continuous. The weak formulation (1) is then well-posed; see, e.g., [45].
By trimming the elements in the background mesh, a mesh for the interior of the domain Ω is obtained: Similarly, meshes for the Dirichlet and Neumann boundaries are defined as: Note that elements in T h ∂Ω D and T h ∂Ω N are manifolds of co-dimension 1.
We consider a B-spline basis of degree k and regularity α constructed over the ambient domain using the Cox-De Boor recursion formula [52]. The span of this B-spline basis is denoted as with P k (K ) the collection of d-variate polynomials on K . The approximation spaces in the finite cell method are obtained by restricting the B-splines in S k α (A) to the domain Ω: A basis for W h , V h follows immediately from the restriction of the B-spline basis for S k α (A).
It is generally infeasible to impose restrictions on W h , V h according to (6) to form subspaces of W 0 , V 0 that retain suitable approximation properties. Hence, the finite cell method generally relies on weak imposition of the Dirichlet boundary conditions, typically by means of Nitsche's method [53]. Accordingly, the bilinear and linear forms for the approximation problem are adapted to weakly incorporate the Dirichlet boundary conditions, giving rise to a Galerkin formulation of the form: where the superscript h on the bilinear form a h explicit dependence on the mesh parameter. In general, the bilinear and linear forms in the finite cell method comprise contributions from both the interior and the boundaries. Corresponding to this typical setting, in the remainder of this work we assume the operators to be of the form where for each (interior or boundary) element K , the set {(x l K , ω l K )} l K l=1 represents a quadrature rule, i.e., a suitable collection of pairs of integration points in K and corresponding weights. We denote by Q the complete integration scheme, i.e., the collection of all the (interior and boundary) element-wise quadrature rules. The Galerkin problem corresponding to the approximate (bi-)linear forms a h Q and b h Q then writes: We assume that the integration scheme Q transfers the coercivity and boundedness properties of a h , b h to a h Q , b h Q , so that (10) is well posed. These assumptions generally pertain to the required order of the integration scheme in relation to the discretized problem (differential equation, basis function orders, etc.) under consideration. The reader is referred to, for example, Strang and Fix [43,Ch. 4] or Ern and Guermond [45,Ch. 8] for a detailed discussion on this topic.

Finite cell error analysis
The error of the approximated finite cell solution, u h Q , computed using the Galerkin problem (10) is composed of two parts, viz.: (i) the discretization error, defined as the difference between the exact solution to (1), u, and the approximate solution to (7) in the absence of integration errors, u h ; and (ii) the inconsistency error related to the integration procedure, which is defined as the difference between the approximate solution in the absence of integration errors, u h , and the approximate solution to (10) with integration errors, u h Q . In this section we present the error analysis in an abstract setting. A concrete example is provided in Section 2.3. To provide a setting for the error analysis, we denote by W (h) = span {u} ⊕ W h the linear space containing the actual solution to (1) and the approximation space , the solution of (1) is continuously embedded in W (h); cf. [45,Ch. 2]. We assume that the solution u of (1) is suitably regular, so that the bilinear form a h admits a continuous extension to W (h), i.e., there exists a constant C > 0 such that: In this setting, an upper bound for the error u − u h Q is provided by Strang's first lemma [44], which, following [45,Lemma 2.27], is written as where I h u represents the best approximation of u in W h , i.e., and α h denotes the inf-sup constant of the bilinear form a h Proof of (12) is standard; see, for instance, [45,Lemma 2.27]. The first term in (12) represents an upper bound to the discretization error. The second term in (12) bounds the consistency error, i.e., the integration error. Eq. (8) conveys that, in principle, integration errors emerge for both the volume integrals and the surface integrals. Since the focus of this work is on the optimization of the integration orders over octree sub-cells, in the remainder of this work we will restrict our considerations to the integration errors associated with the volume integrals, thereby implicitly assuming that the boundary integrals are evaluated exactly. When neglecting the integration errors associated with the boundary integrals, the error term in (12) associated with the inexact integration of the linear form is bounded by Note that the supremum is always considered over a function space with the zero function excluded, e.g., v h ∈ V h \{0}. For notational brevity we omit the zero exclusion in the remainder of this manuscript. Under the (non-restrictive) assumption that B h Ω is a local operator in the sense that e.g., if B h Ω corresponds to a differential operator, the summands in the ultimate expression in (15) are bounded by the element-integration-error indicators The element-integration-error indicators e b K can be computed element-wise and provide a bound for the integration error in the linear form: Similarly, the integration error for the bilinear form in (12) is bounded as with the element-integration-error indicators defined as Substitution of the bounds (18) and (19) into the error estimate (12) then yields which conveys that the element-integration-error indicators (17) and (20) yield control over the inconsistency error. This motivates the development of an algorithm to minimize the element-integration-error indicators; see Section 3.2.

Application to the Poisson problem
We apply the general theory presented in Sections 2.1-2.2 to the particular case of the Dirichlet-Poisson problem: with f : Ω → R and g : ∂Ω → R exogenous data. A suitable ambient space for the weak formulation of (22) is provided by W = H 1 (Ω), equipped with the usual H 1 -norm and inner product. Denoting by ℓ g ∈ H 1 (Ω) a lift of the Dirichlet data such that ℓ g | ∂Ω = g and by W 0 = H 1 0 (Ω) the subspace of functions in H 1 (Ω) that vanish at the boundary in the trace sense, the weak formulation of (22) assumes the form (1) with We denote by V h a finite-dimensional subspace of H 1 (Ω) according to the construction in Section 2.1. We equip V h with the mesh-dependent norm The symmetric Galerkin formulation in V h with weak enforcement of the Dirichlet conditions by means of Nitsche's method conforms to (7) with a h For suitably chosen constants c 0 and c 1 , the bilinear form is bounded and coercive on V h × V h and the linear form b h is bounded on V h (see Refs. [54,55] (14) apply. To elucidate the notation for the quadrature-dependent bilinear and linear forms for the Dirichlet-Poisson problem, we note that in this case: One easily verifies that both A h Ω and B h Ω satisfy the locality condition (16). Using these expressions, the element-errorindicators for the linear form (17) and the bilinear form (20) follow as: Substitution of these error indicators into (21) finally yields the H 1 -error estimate of the Poisson problem.

Optimized octree cut-cell integration
Based on the error analysis discussed in Section 2.2, in this section we develop an algorithm to optimize the distribution of integration points over octree-subdivided cut-cells. In Section 3.1 we first introduce the considered octree procedure and study its complexity with respect to the most relevant parameters. In Section 3.2 we then express the integration-error estimate in an operator-independent form, making it applicable to a class of operators. The per-cut-element evaluation of the operator-independent integration error is also discussed in this section. Finally, the developed optimization procedure in the form of Algorithm 1 is presented in Section 3.3.

Octree partitioning
To construct an explicit parametrization of the geometry (including its boundary) we herein consider the octree procedure proposed in Ref. [7], which is illustrated in Fig. 3 for a two-dimensional cut-cell. In this procedure an element in the background mesh that intersects the boundary of the domain is bisected into 2 d sub-cells. If a sub-cell is completely inside the domain, it is preserved in the partitioning of the cut-cell, whereas a sub-cell is discarded if it is completely outside of the domain. This bisectioning procedure is recursively applied to the sub-cells that intersect the boundary, until ϱ max -times bisected sub-cells are obtained. At the lowest level of bisectioning, i.e., for the ϱ max -times bisected sub-cells, a tessellation procedure is applied to construct a partitioning of the sub-cells that intersect with the domain boundary, resulting in an additional level, i.e., ϱ max +1. This tessellation procedure provides an O(h 2 /2 2ϱmax ) accurate parametrization of the interior volume [7]. Additionally, it provides a parametrization for the trimmed surface. See Appendix for details regarding the employed tessellation procedure.
Since the octree procedure provides a parametrization of the cut-cell and integration rules are generally available for all sub-cells [27,51], cut-cell integration can be performed by agglomeration of all sub-cell quadrature points. The accuracy of the cut-cell integration scheme can then be controlled through the selection of the quadrature rules on the sub-cells. In particular, polynomials can be integrated exactly over the cut-cell when Gauss quadrature of the appropriate order is used on all sub-cells. An example of a cut-cell integration scheme based on Gauss quadrature for third order polynomials is illustrated in Fig. 3, where the gray squares represent the Gauss points, and the relative size of the squares is representative for the integration weights.
Selection of Gauss points of optimal order on all sub-cells is evidently very attractive, in the sense that an adequate cut-cell integration scheme can be constructed by specification of the bisectioning depth ϱ max and the order of polynomials to be integrated exactly. However, as observed from Fig. 3, the obtained cut-cell integration rule is generally not computationally efficient, in the sense that the majority of integration points is formed on small sub-cells, as a consequence of which their relative contribution to the overall integral (observed from the size of the squares) is generally limited. As studied in Ref. [27], the computational efficiency of the cut-cell integration scheme can be improved by reducing the number of integration points on lower bisectioning levels. Although this reduction decreases the accuracy of the integration scheme, the obtained improvement in computational effort associated with the lower number of integration points outweighs this disadvantage. The algorithm proposed in Ref. [27] targets the optimization of this balance between accuracy and computational effort.
Since the balance between accuracy and computational effort also forms the basis of the error-estimation-based optimization procedure in this work it is important to understand the distribution of the number of octree sub-cells over the levels of bisectioning. To establish a scaling relation for the number of sub-cells we consider a random cut-cell of size h with trimmed surface area S, as illustrated in Fig. 4(a). The octree subdivision approximates the surface area at the maximum octree depth as where m 0 ϱmax denotes the average number of sub-cells at octree depth ϱ max that is intersected by the trimmed boundary, and wheres is the average surface area of a randomly cut unit cube in d dimensions. This surface approximation is valid under the condition that with surface fraction ς = S sh d−1 . For example, for the cut-cell illustrated in Fig. 4, ϱ max ≥ ϱ min = 1. From (28) the number of intersected elements at level ϱ max follows to be: Denoting the number of sub-cells at level ϱ that are completely inside of the domain by m + ϱ (see Fig. 4), and noting that the average volume of a randomly cut unit cube is 1 2 , the cut-cell volume can be approximated as which is divided by the volume of the background cell to obtain the cut-cell volume fraction: Assuming that ϱ max adequately resolves the volume fraction of the cut-cell, it holds that η ϱmax ≈ η ϱmax−1 with ϱ max ≥ ϱ min + 1 and Substitution of (32b) into the above equation yields from which the number of sub-cells at bisectioning level ϱ max can be obtained as Since it follows from (30) that the interface localizes in a single sub-cell if ϱ ≤ ϱ min (m 0 the number of preserved sub-cells across all levels is given by where the volume fraction η dictates which side of the interface corresponds to the interior domain. For example, for the cut element displayed in Fig. 4 the volume fraction is larger than a half, which implies that all but one of the sub-cells at level ϱ = 1 are preserved. If the complement of the element would be considered, none of the level ϱ = 1 sub-cells would be included in the partitioning.
If we denote the average number of integration points per octree sub-cell byq ϱ , the total number of integration points per level follows as wheret d is the number of sub-cells in which the lowest bisectioning level ϱ = ϱ max + 1 is tessellated (see Appendix). The total number of integration points for a level ϱ max octree partitioning then follows as This approximation reflects the strong dependence of the number of integration points on the surface fraction ς . When ς is very small, ϱ min is very large, as a result of which a significant number of integration points for non-intersected sub-cells can be present (depending on the volume fraction). When ς is large, or, more generally when ϱ max ≫ ϱ min , the number of non-intersected sub-cells is negligible. For the case that the Gauss order is selected equal on all levels, i.e.,q ϱ = q d line for ϱ ≤ ϱ max (with q line the number of integration points along a one-dimensional line segment) andq ϱmax+1 = q d tes for ϱ max + 1 (with q d tes the number of integration points in a tessellated sub-cell), the number of integration points (38) reduces to: This estimate conveys that the number of integration points scales linearly with the surface fraction, and exponentially with the octree depth. This exponential scaling depends on the number of spatial dimensions. In two dimensions the number of integration points doubles while increasing the octree depth by one, while in three dimensions it quadruples. Hence, in comparison to the two-dimensional setting, the number of integration points in three dimensions scales much more dramatically with the octree depth. A numerical study of the number of integration points in relation to the scaling relations presented in this section is provided in Section 4.

Cut-cell integration errors
In Section 2.2 the application of Strang's first lemma to the quantification of integration errors in the finite-cell method was discussed. Eqs. (18) and (19) convey that the integration error can be estimated by the summation of the operatordependent element-integration-error indicators e b K and e a K . This section discusses the evaluation of these errors and their localization to the integration sub-cells. In order to evaluate the cell-wise integration error, in Section 3.2.1 the integration errors are expressed as the product of an operator-independent integration error and a scaling term depending on the operator. Subsequently, in Section 3.2.2 the numerical evaluation of the operator-independent cut-cell integration error is discussed.

Integration error definition
The element-integration-error indicators e b K and e a K , expressed by Eqs. (17) and (20) respectively, depend on the operations A h Ω and B h Ω . It is generally desirable to apply a single integration scheme for all terms and, hence, to have uniform control over e a K and e b K . To enable a uniform treatment of both e a K and e b K , we first note that the integrals in the numerators of (17) and (20) constitute linear functionals on V h K . By the Riesz-representation theorem, there exist elements We now proceed under the assumption that the error introduced in the numerical integration is equivalent for the original functionals and their Riesz-representation, i.e., the difference in applying the integral rule to the left-and right-hand-side members of (40) and (41) is negligible. It then holds that for (20) and the following similar expression applies to (17) For T a , v h K in the polynomial space V h K , the product T a v h K resides in the double-degree polynomial space P K . By virtue of the equivalence of norms in finite dimensional spaces, P K can be equipped with a norm ∥ · ∥ P K such that Substitution of the operation-independent element integration errors (45) in the global error bound (21) yields The error bound (46) conveys that the integration errors can be controlled through the operator-independent elemental integration errors e p K .

Integration error evaluation
The element integration error (45) can be approximated by discretization of the space P K . We select the finite dimensional space P K as the tensor-product space of univariate polynomials of order k over the element K , i.e., with Φ(x K ) ∈ R np the basis of monomials of size n p = (k + 1) d and a ∈ R np (and, correspondingly, γ i,j ) the corresponding coefficients. The space P K is equipped with the Sobolev norm ∥ · ∥ H r with r ≥ 0. Using the polynomial basis (47), the element integration error (45) can be expressed as where λ max is the largest eigenvalue of the generalized eigenvalue problem 2 with eigenvalues λ i and eigenvectors v i , and with the vectors ξ andξ defined as The matrix G in (49) is defined as the Gramian matrix for the basis functions Φ associated with the H r (K ) Sobolev space.
Since the left-hand-side matrix in the eigenvalue problem (49) is the dyadic product of two vectors its rank is equal to one. As a consequence, the eigenvalue problem (49) only has one non-zero eigenvalue, with the corresponding eigenvector being in the direction G −1 (ξ −ξ). Since the Gramian matrix is positive definite, this single non-zero eigenvalue is positive and therefore equal to the maximum eigenvalue. This maximum eigenvalue and its corresponding eigenvector are equal to Substitution of λ max in Eq. (48) expresses the element integration error in terms of the vector ξ −ξ as: The polynomial corresponding to this error, i.e., the function in P K yielding the largest integration error follows directly from the maximum eigenvector v max as Note that the scaling of this function with respect to the norm (52) is chosen such that 2 Using a = ∑ np i=1ã i v i = Vã, with v i the eigenvectors in (49) normalized w.r.t. the Gramian matrix G, the expression in (48c) can be written as:

The cut element quadrature optimization algorithm
The ability to evaluate the maximal cut-cell integration error and the corresponding integrand that leads to this error serves as the basis for the quadrature-optimization procedure developed in this work. The error (56) is, in general, controlled by the octree depth, ϱ max , and the quadrature rules selected on all sub-cells. In the remainder of this work we assume the octree depth to be fixed at a level where the partitioning error is negligible.
For the space of polynomials considered above, i.e., those of Eq. (47), the vector ξ of basis function integrals can be evaluated exactly on the exact or approximated geometry and by appropriate selection of the Gauss integration order. As discussed in Section 3.1, the usage of exact Gauss integration on all sub-cells is impractical from a computational effort point of view in large scale finite-cell simulations. The conceptual idea behind the quadrature optimization procedure developed here is to determine integration rules with a significantly smaller number of integration points, and to use the evaluable integration error (56) to find the optimal distribution of these points over the sub-cells.
The developed optimization procedure is intended as a per-element pre-processing operation, which results in optimized quadrature rules for all cut elements in a finite cell simulation. In order to conduct this pre-processing operation, for each cut element the exact shape function integrals ξ and the Gramian G must be determined. These computations are relatively expensive, but especially in time-dependent or non-linear problems -where the optimized quadrature scheme will be re-used many times -the computational gain in the simulations outweighs the effort in this pre-processing operation. The computational effort will be studied in Section 4.3.1.
In Section 3.3.1 we first formalize the integration error minimization problem that serves as the basis for the optimization algorithm. In Section 3.3.2 the algorithm is then introduced and various algorithmic aspects are discussed.

The minimization problem
We consider an octree partitioning of the element K as described in Section 3.1 and denote this partitioning by ,ϱ the sub-cells corresponding to the octree level ϱ. On each sub-cell ℘ ∈ P a sequence of quadrature rules is defined: We assume that the quadrature rules Q ı ℘ are nested, in the sense that the set of polynomials on ℘ that are integrated exactly by Q ı+1 ℘ includes all polynomials that are integrated exactly by Q ı ℘ . The quadrature rule for the complete partitioning is defined as with ı = {ı ℘ | ℘ ∈ P} a list of length m = #P with sub-cell quadrature rule indices. For a given type of integration rule, e.g., Gauss integration, the index list ı ∈ N m fully determines the quadrature rule for the partitioned cut-cell K .
Both the total number of integration points in the cut element, defined by #Q ı ℘ , and the integration error (56) can be expressed as functions of the quadrature index list ı: with p K ,max according to (54) for the integration rule Q ı P . Using these functions the intended quadrature optimization procedure can be formulated as the constrained minimization problem minimize ı∈N m e p,ı K subject to #Q ı with q ⋆ the specified number of integration points. Conversely, we can express the optimization problem in terms of the minimization of the number of integration points for a fixed error with e ⋆ K the intended error level.

Remark 1 (Quadrature Refinement Strategy).
In this work, the integration quadrature per sub-cell is refined by increasing the order of the integration scheme. Alternatively, refinement would be possible by subdivision of a sub-cell and then constructing quadrature rules over the refined sub-cells. Since we consider the integration of smooth functions over the cut-cells, improving the quadrature by increasing its order is more efficient than spatial refinement of the sub-cells.

The optimization algorithm
The developed algorithm to obtain the optimal distribution of integration points over all sub-cells is presented in Algorithm 1. This algorithm approximates the minimization problem (60) through the generation of a sequence of refinement schemes, {ı 0 , ı 1 , . . . , ı ⋆ }, where for the initial quadrature rule the lowest order of integration on each sub-cell is considered, i.e., ı 0 = 0, and where the optimized quadrature rule (approximately) satisfying the constraint condition is denoted by i ⋆ . Given the rth iterate in the optimization procedure, ı r , the next integration rule in the sequence, ı r+1 , is determined in three steps: 1. The function corresponding to the maximum integration error, p K ,max , is determined for an integration rule ı r using Eq. (54). The corresponding integration error can be found using (59a). This step corresponds to the lines These steps are repeated until the specified number of integration points is reached, i.e., when #Q ı ℘ ≥ q ⋆ . The algorithm is stopped prematurely when the integration rule sequence (57) is depleted. Specifically, in the considered implementation [56] the maximum quadrature order for the lowest level tessellation is equal to 6 (12 points) for triangles and 7 (31 points) for tetrahedrons. The algorithm termination conditions correspond to line (24) in Algorithm 1.

Sub-cell indicators.
To form the sub-cell indicators we consider a uniform refinement of the integration indices, i.e., ı r+1 = ı r + 1 for each sub-cell. The error reduction per added integration point in the case of this uniform refinement can be expressed as Using the sub-cell integration error according to and the reduction in the number of integration points, 3 the error reduction (62) is bounded as: with the sub-cell indicator function defined as The sub-cell indicator (65) weighs the local integration error by the cost of increasing the integration order. Under the assumption that the sub-cell integration error reduces significantly when increasing the order of the integration scheme 3 #Q ı r+1 by one, the numerator in (65) can be interpreted as the error reduction rather than the error itself. This assumption of local error reduction is, however, in general, not satisfied. As a consequence, in the proposed algorithm sub-cells whose integration error does not decrease significantly (or even increases) by raising the order of integration by one would be underrated in the indicator function. Considering the error rather than the error reduction in the indicator provides a more robust measure of the potential to decrease the error in a sub-cell.
Sub-cell marking strategy. Based on the upper bound (64) of the error reduction per added integration point, sub-cells are marked for refinement. We herein consider two marking strategies, with M denoting the set of sub-cells marked for increasing the integration order: • Sub-cell marking: In this marking strategy we increase the integration order of the single sub-cell with the largest indicator: • Octree level marking: Sub-cell marking involves refinement of individual sub-cell in each refinement step. This generally leads to a large number of optimization steps, which is undesirable from a computational effort point of view. In order to expedite the optimization procedure, it is possible to mark multiple sub-cells simultaneously. A natural way of marking multiple sub-cells is to agglomerate all sub-cells in a particular level. This leads to a marking strategy in which we increase the integration order for the level (in the octree partitioning) with the largest indicator: Remark 2 (Global Adaptive Integration). The algorithm presented above is framed in an element-by-element setting, i.e., each cut element is considered separately in the pre-processing operation to determine the integration schemes. Since the evaluation of the integration error and its indicators involves the repeated computation of solutions to a linear system with the size equal to the number of supported basis functions, i.e., Eq. (56), global application of the algorithm is impractical from a computational effort point of view. It should be noted that the element-by-element setting considered here does not account for the fact that the operatordependent constants in the multiplicative decomposition (45) in principle vary per element. For example, a source term might be negligibly small on an element that is challenging to integrate, and hence improving quadrature on that element would be inefficient from the error approximation point of view. Additional integration points are, however, assigned to such a cut element in the element-by-element strategy employed here, as the operator is not considered in the marking strategy.
Global adaptive integration is feasible, however, by computing the indicators on all elements individually (i.e., solving multiple relatively small local linear systems, rather than a large global system), and then to apply a global marking strategy. Both the sub-cell marking and the level marking strategy can then be applied globally. The potential benefit of such a global marking strategy is that it automatically accounts for the fact that not all cut elements are equally hard to integrate. However, in principle, the operator-dependent constants should then be incorporated in the error indicators. Moreover, one should be aware that parallelization of the pre-processing operation, which is trivial in the per-element setting, is more challenging using the global marking strategy, as communication between the elements is required. The development of a rigorous global integration optimization routine, which would consider the above-mentioned complications, is considered beyond the scope of the current work.

Numerical study of the adaptive integration procedure
To assess the performance of the developed adaptive integration technique, in this section we study its behavior in terms of integration accuracy versus the number of integration points. We here consider the case of a single cut-cell. The effect of the integration accuracy on actual finite-cell simulations will be studied subsequently in Section 5.
We consider a d-dimensional unit cube [0, 1] d in two and three dimensions as a single element of the background mesh. Throughout this section, with the exception of Section 4.3.3, we exclude a circle in the two-dimensional cases and a sphere in the three-dimensional cases. The considered exclusions are centered at the origin of the background element, as illustrated in Fig. 5.

Equal-order degree integration
Before we proceed with the presentation of the results obtained using the adaptive integration procedure, we first consider the case for which the order of the integration scheme is chosen to be same over all integration sub-cells. We consider the case of a circular exclusion with r = 0.6 for two dimensions and spherical exclusion with r = 0.6 for three    In the remainder of this manuscript the integration error is computed with respect to the exact integral over the considered partitioning. Hence, the geometric error corresponding to the octree partitioning is not represented in the studied error results. The rational behind this choice is that immersed finite element methods (specifically the finite cell method) allow to control the geometry error independently of the considered grid, and hence that the geometry error can be made insignificant in comparison to the other error contributions. To illustrate the effect of the geometry error in relation to the integration error, Fig. 6(b) plots the approximation error versus the number of integration points under uniform increase of the integration order, for octree depths ϱ max = 0, 1, 2, 3, 4. In panel 6(a), the reference result is determined with very large ϱ max and quadrature order. The approximation error hence contains both a geometric component and a quadrature component. In panel 6(b), the reference result is determined for each ϱ max with very large quadrature order. Correspondingly, the approximation error comprises only a quadrature component. The isolated integration error converges to machine precision upon increasing the order of the Gauss integration schemes, the geometry error does not. This is observed clearly from Fig. 6(a), where the combined error levels off at the geometry error upon increasing the integration degree. The geometry error, i.e., the value to which the combined error converges upon resolving the integration error, is observed to converge with the theoretical rate of O(2 −2ϱmax ) [7]. Fig. 6(b) conveys that the integration errors are representative of the combined error provided that the octree depth is set sufficiently large. Throughout this manuscript it is assumed that this condition is satisfied, which indeed is the case for the integration schemes that are of practical interest from a computational effort point of view. The reader should be aware, however, that in our analyses we generally improve our integration schemes until the integration errors are fully resolved. In these limiting cases, the geometry error will of course be the dominant factor in the combined error, although it is not observed in our plots. Fig. 7 displays integration errors based on both equal-order Gauss quadrature and equal-order uniform quadrature (i.e., the trapezoidal rule with equidistantly-spaced quadrature points) on all sub-cells. It is plotted for the cases of ϱ max = 3, 4, 5 in two dimensions and ϱ max = 2, 3, 4 in three dimensions. It is observed that for all ϱ max the observed errors for a particular integration order are similar, which is explained by the fact that the integration error is dominated by contributions from the sub-cells that are already present at the lowest ϱ max and that the errors are computed with respect to the considered partitioning. In terms of the number of integration points, the octree depth does, however, have a substantial effect. The observed increase in number of integration points is studied in detail in Fig. 8 for the case of a uniform scheme with 4 2 points per sub-cell in two dimensions and 4 3 points per sub-cell in three dimensions and for a Gauss scheme of order 4 on each sub-cell. This figure conveys that the total number of integration points scales in agreement with the relation (38). In particular the doubling of the number of points in two dimensions and the quadrupling in three dimensions while increasing the octree by a single level is clearly observed.
From Fig. 7 it is observed that for Gauss orders that are significantly below the order required for exact integration there is not a substantial difference with uniform integration. Once the Gauss order reaches that needed for exact integration, the error observed for Gauss integration is substantially smaller than that using uniform integration. When considering the same integration order on all sub-cells, this advantage is, however, only attained at the expense of introducing a large number of integration points, in particular in cases of high octree depths.

Adaptive integration
We now consider the adaptive-integration procedure for the circular exclusion setting considered in Section 4.1. In Fig. 9 various steps in the quadrature-optimization procedure are shown, displaying for each step the number of integration points per sub-cell (left column), error and the function leading to that error (middle column), and the subcell error indicators (56) (right column). The presented results are generated using the per sub-cell marking strategy. The sequence of steps demonstrates that in each step the worst possible function in terms of integration is determined. In the first step, this corresponds to a function that is large in magnitude at the largest sub-cell. This function leads to an indicator that is highest in that sub-cell, and hence the order of integration on that cell is increased. As a result, in the second step a function that is large in magnitude on the sub-cells surrounding the largest sub-cell is found to be the worst possible in terms of integration error. The larger volume of the ϱ = 1 integration cell makes, however, that the largest indicator is still found for that cell. After another increase in integration order on that sub-cell, in the subsequent steps the largest indicators are found on the ϱ = 2 sub-cells, which are therefore gradually increased in order.  The displayed results pertain to ϱ max = 3 in both two and three dimensions, displaying the equal-order results discussed in Section 4.1 for reference. As can be seen, the error per integration point is substantially lower using the adaptive integration procedure. For example, in two dimensions, the error corresponding to the equal second order Gauss scheme is equal to 2.52 × 10 −2 , while that particular integration scheme involves 144 points. For the same number of points, the error corresponding to the optimized quadrature is equal to 1.00 × 10 −3 , i.e., a factor 25 reduction in error. Fig. 11 displays the distribution of the integration points over the sub-cells for the equal-order Gauss scheme and the optimized case, which clearly demonstrates that the significant reduction in error is achieved by assigning more integration points to the larger sub-cells before introducing additional points in the smaller sub-cells.
In general, for a fixed integration error, the number of integration points is substantially lower using the adaptive integration procedure than when using the equal-order schemes. For example, in two dimensions, the error corresponding to the equal-order Gauss scheme with degree 4 is equal to 7.35 × 10 −3 , which involves 303 points. For a similar error, the number of points for the optimized quadrature is equal to 83, i.e., a reduction by a factor of approximately 4. Fig. 12 displays the distribution of integration points over the sub-cells for the equal-order Gauss scheme and the optimized quadrature, which clearly shows the reduction in number of points for a fixed error. This figure also conveys that the observed reduction factor of 4 is significantly influenced by the initial integration order in the adaptive procedure, which in this case has been set to 1 (for this case the minimum number of points is equal to 43). Fig. 10(b) displays the error versus the number of integration points for the three-dimensional spherical exclusion with ϱ max = 3. The integration error corresponding to the equal-order Gauss scheme with degree 2 is equal to approximately 1.14 × 10 −2 and pertains to 7168 integration points. For the same number of points, the adaptive procedure reduces the error to 2.67 × 10 −5 , which constitutes a reduction factor of approximately 450. This reduction factor is significantly higher than that observed in the two-dimensional case, despite the similarity in asymptotic scaling rate between the integration error and the number of integration points in the two-and three-dimensional cases. The larger reduction factor in three dimensions is attributed to the significant reduction in error that is achieved during the first integration point optimization steps, in which the integration order on the ϱ = 1 sub-cells is increased. The difference between the two-dimensional and three-dimensional setting in this regard is the fact that in three dimensions there is a significant number of such sub-cells, whereas in two dimensions there is only one sub-cell of level ϱ = 1. For an integration error 1.139 × 10 −2 , the adaptive procedure requires only 1646 points, which, as in the two-dimensional case, is a reduction by a factor of approximately 4. As for the two-dimensional case, this reduction factor is significantly influenced by the integration order with which the adaptive procedure is initiated.

Influence of algorithmic settings and cut-cell configurations
In the previous section the performance of the adaptive integration routine has been studied using representative settings for the cut-cell geometry and algorithmic parameters. In this section we will study the influence of the most prominent parameters on the performance of the adaptive integration technique.

The marking strategy
From the refinement patterns that emerge from the per sub-cell refinement strategies -such as the ones discussed in the previous section -it is observed that, as a general trend, integration orders are increased on a per-level basis, i.e., the number of integration points is increased from level ϱ to ϱ max + 1. This is explained by the fact that the indicators scale with the volume of the sub-cells. Based on this observation it is anticipated that the level-based marking strategy discussed in can be very efficient, in the sense that it yields a similar refinement pattern as the per-cell marking, but that it needs fewer iterations by virtue of marking a larger number of sub-cells per step.
In Fig. 13 the per-level and per-cell marking strategies are compared for the test case introduced above. In both two and three dimensions it is observed that the per-level marking strategy closely follows the per-cell marking. In particular in the initial steps of the optimization a very close agreement is observed between the marking strategies. In Fig. 14 we inspect the distribution of points corresponding to the two marking strategies with 204 points (in the two-dimensional setting). This figure conveys that, indeed, the distribution of integration points between the two markings is very similar.
Although the observed point distributions in Fig. 14 are similar, the number of iterations required to attain these distributions is significantly different. The per-cell marking strategy requires 54 iterations, whereas the corresponding result using the per-level marking is achieved in 9 iterations. Evidently, this reduction in number of iterations translates into a computational-effort advantage for the per-level marking strategy. Table 1 Analysis of the CPU time to attain the results displayed in Fig. 13 In Table 1 we study the computational effort for both marking strategies corresponding to the optimization procedure to attain the results in Fig. 13(a). The reported CPU times are based on a Python implementation using the open source finite element toolbox Nutils [56], which is executed on a four core CPU with a 2.60 GHz 7th generation Intel Core i5 processor. Although the reported CPU times are highly dependent on a myriad of aspects, it is noted that the computation time of the Gramian matrix (0.4-0.45 s) is representative for the computational effort involved in the construction of the element contribution to the system matrix. By relating the reported CPU times to this element system matrix construction time, a more meaningful notion of the computational effort is obtained, although, of course, also this relative notion of performance is subject to implementation and architecture considerations. Table 1 conveys that the optimization algorithm has a similar initialization time for both marking strategies. This initialization time pertains to the computation of the reference basis function integrals (50) and the Gramian matrix (G). Moreover, it is observed that the CPU time per step is also very similar between the two marking strategies. The majority of the computational effort per step resides in the computation of the basis function integrals,ξ, and in the evaluation of the integration error localized to the sub-cells. In total, the per-level marking strategy is, however, computationally more efficient on account of the significant reduction in number of iterations.
The computation time for the determination of the optimized integration scheme depends primarily on the number of sub-cells and on the dimension of the considered polynomial space. In principle, the computation time of the integral evaluations in the algorithm -i.e., the computation of the basis function integrals, the Gramian, and the localized errors -scales linearly with the number of sub-cells. It should be noted that the number of sub-cells is highly dependent on the bisectioning level and on the number of dimensions; see Eq. (36). For example, for the two-dimensional case reported in Fig. 13(a) and Table 1 the number of sub-cells is equal to 43, whereas the three-dimensional case in Fig. 13(b) has 341 sub-cells. Since the integral evaluations dominate the overall CPU time of the initialization phase and that of an optimization step, these operations become proportionally more expensive. It is important to note, however, that the number of iterations scales with the number of sub-cells in the case of the per sub-cell marking, whereas it scales with the number of octree levels in the case of the per-level marking. This implies that, with a growing number of sub-cells, the per-level marking strategy becomes more favorable from a computational effort point of view.
In terms of computational effort, the dimension of the polynomial space primarily has an effect on the computation time for the basis function integrals and that for the Gramian. In particular the Gramian computation becomes more expensive with an increase in basis size, as the number of terms to be integrated scales quadratically with this size. The linear system solving step involved in the computation of the worst possible function also increases with an increase in system size, but for all considered systems the computational effort involved in the employed direct solver remains negligible compared to the integral evaluations.

Influence of the functional setting
From the functional-setting point of view, the distribution of integration points is influenced by both the approximation order of the polynomial space (47) and by the norm (54) in which the integrands are considered, encoded by the Gramian G.
In Fig. 15 we study the influence of the integration error for various orders of the integrands. Note that the lower bound on the error at approximately 10 −16 corresponds to machine precision errors, and hence these integrands can be interpreted to be exact. As expected, with an increase of the integrand order, a larger number of points is required to attain a certain accuracy. From the relatively low order settings, in particular the case of linear functions, it is observed that the algorithm quickly reaches machine precision, as all sub-cells are marked up until the degree required for exact polynomial integration. Note that in the initial setting, i.e., one point per sub-cell, the k = 1 integrands are not resolved exactly, because the cross term x 1 x 2 (in two dimensions) is not integrated exactly on the triangulated sub-cells. The same applies to the three-dimensional setting. In Fig. 16 the point distributions for the two dimensional case corresponding to approximately 90 points are shown for various orders. It is observed that the distribution of the points over the levels is initially rather uniform for the k = 2 case, but becomes more dispersed as the order increases. The reason for this is that the higher-order integrands encompass functions that are more localized to the bigger sub-cells than the lower-order integrands. For a fixed number of integration points, the distribution of integration points over the elements is observed to converge with increasing integrand orders.
In Fig. 17 the influence of the function space norm (54) is studied. Although the normalization does affect the overall magnitude of the error, the difference in norm is observed to have a minor effect on the distribution of the integration points. This is confirmed in Fig. 18 for the two-dimensional setting, where it is observed that for a total number of 90 points, the distribution using either the L 2 -norm or the H 1 -norm is virtually identical.

Influence of the cut-element geometry
The studies presented above were all conducted for cut-elements with spherical exclusions. To study the influence of the geometry of the cut-elements, we here consider elliptic exclusions with varying ellipticity and orientation. In two dimensions we consider an exclusion of an ellipse with semi-major axis r 1 and semi-minor axis r 2 (with the sphere as the special case r 1 = r 2 ) that is centered at the origin of the background element and with its major axis residing in the x 1 − x 2 plane at an inclination of ϕ with respect to the x 1 axis. In three dimensions we consider an exclusion based on the same ellipse, which is extruded along the x 3 direction. The considered cut-elements are illustrated in Fig. 19.
In Fig. 20 we consider a range of exclusions in a two-dimensional cut-element with differing ellipticity. The corresponding error per integration point using the per-level adaptive integration procedure is displayed in Fig. 21. It is noted that the rate with which the error decreases is similar for all geometries, but that the number of points to attain a particular error is geometry dependent. This behavior is explained by the fact that although the distribution of the orders over the levels is very similar for each of the geometries, the number of sub-cells in each level is geometry dependent. As elaborated in Section 3.1 the number of cells on each level scales with the surface fraction. This scaling is reflected in the results in Fig. 20. For example, the surface-to-volume ratio of Fig. 20(c) is approximately two times that of Fig. 20(a). An increase in number of integration points by a factor of approximately 2 is also observed from the offset of the corresponding curves in Fig. 21. Fig. 22 shows the integration error versus the number of integration points for a cut-cell with r 1 = 0.6 and r 2 = 0.1 for various inclination angles ϕ. This study confirms the scaling relation with the surface to volume ratio as observed for the ellipticity variations considered above.

Manually selected quadrature rules
Although the computational effort involved in the construction of the optimized quadrature rules is in general acceptable when one wants to re-use a quadrature rule multiple times, some computational effort is involved in this construction. In addition, one has to set up a suitable code to determine the optimal distributions for arbitrary cut-cells. Considering this, one may not be interested in obtaining the optimized distributions of the points, but may instead want a simple rule of thumb to select the quadrature on a cut-cell; see, e.g., Refs. [27,40].
The per-level selection of the integration order makes it practical to manually select integration rules that outperform full order integration on all octree levels. We here consider two manual selection strategies based on Gauss integration of the sub-cells: A. Minimal degree lowering: In this strategy we set the order of the Gauss scheme on the level ϱ = 1 sub-cells to k max .
We then reduce the number of Gauss points per direction by one for each level, which implies that for ϱ ≤ ϱ max  the Gauss degree is decreased by two between two octree levels. The integration order of the tessellated sub-cells at level ϱ max + 1 is set to be two orders lower than that at level ϱ max . Once a Gauss degree of zero is reached this value is maintained for the underlying levels. For example, for the case of k = 8 and ϱ max = 3, the Gauss orders over the levels are set to 8 for ϱ = 1, 6 for ϱ = 2, 4 for the ϱ = ϱ max , and 2 for the ϱ = ϱ max + 1 (see Fig. 24(a)). That is, for an octree-depth of ϱ max = 3, the list of levels [1,2,3,4] contains [8,6,4,2] as the corresponding list of the Gauss orders.
B. Uniform degree lowering: This strategy also starts with selecting the ϱ = 1 integration degree as k max , but then decreases the degree between two levels in such a way that the single point (degree is zero) is reached at the levels ϱ max , and ϱ max + 1. For the case considered above, the orders are set to 8 for ϱ = 1, 4 for ϱ = 2, and 0 for ϱ = ϱ max   and ϱ = ϱ max + 1 (see Fig. 24(c)). That is, for an octree-depth of ϱ max = 3, the list of levels [1,2,3,4] contains [8,4,0,0] as the corresponding list of the Gauss orders.
Evidently, alternative quadrature rules can be formulated, but a detailed study of such rules is beyond the scope of the current work.
The adaptive algorithm developed in this work allows us to assess the suitability of the rules of thumb defined above. As a preliminary study on the effectiveness of these rules of thumb, we again consider the case of a two-dimensional circular exclusion with ϱ max = 3 and a polynomial function of order k = 8. In Fig. 23(a) the manually selected schemes are compared to the full order integration schemes and to the per sub-cell optimized integration schemes. This plot conveys that although the manually selected schemes are, as expected, outperformed by the optimized schemes, they generally do provide a dramatic improvement in accuracy per integration point for a fixed number of points. This observed behavior is explained by consideration of the distributions of orders over the sub-cells as shown in Fig. 24, from which it is observed that the per-level decrease of the order of both rules of thumb indeed qualitatively matches the results of the optimization procedure. Evidently, a difference between the two rules of thumb is that the minimal degree lowering strategy (A.) results in a larger number of points and a lower error compared to the uniform degree lowering (B.). Both schemes do, however, reasonably well resemble the optimized distribution patterns. Similar observations for the three-dimensional setting are obtained, as illustrated in Fig. 23(b). The displayed three-dimensional results are based on ϱ max = 3 and k = 5.

Application to immersed isogeometric analysis
In this section we assess the suitability of the optimized integration schemes in the context of immersogeometric analysis. We consider an elasticity problem in two and three dimensions, which we solve using globally defined values of the integration degrees for each octree level. The integration degrees are obtained by application of the optimization procedure to the entire mesh, starting with a single integration point in all sub-cells. The integration errors computed per octree level are summed over all elements, after which the level with the highest global error is refined in terms of integration degree. It is noted that this global optimization strategy ignores operator-dependent variations between the elements, as discussed in Remark 2.
We consider a linear elasticity problem on a computational domain that is immersed into a background mesh of size L d ; see Fig. 25 for d = 2. Displacements are prescribed on the exterior boundary, while the interior boundary is traction free. In the absence of inertia effects and body forces, the boundary value problem reads as with stress tensor σ(u) = λdiv(u)I + 2µ∇ s u, where ∇ s denotes the symmetric gradient operator. Throughout this section, the Lamé parameters are set to λ = 1 2 and µ = 1 2 . The Dirichlet condition on the external boundary is applied strongly, i.e., by constraining the degrees of freedom related to the boundary displacements. The Galerkin problem corresponding to (68) then follows as with the discrete trial space being a subset of H 1 (Ω) and satisfying the Dirichlet boundary conditions. The discrete test space is identical to the trial space, modulo inhomogeneous boundary conditions. All results in this section pertain to full regularity, C p−1 , B-splines of degree p = 2. The polynomial order used to compute the optimized integration schemes is taken equal to k = 4.
As a quantity of interest we consider the constrained modulus where F is the resultant reaction force in normal direction along the displaced boundary, and where A is the area (length in two dimensions) at which this resultant force acts. We normalize the computed constrained modulus by that of a homogeneous square with the same material properties as the considered model, i.e., by M = λ + 2µ. For both the two and three dimensional test cases we consider three levels of bisectioning, i.e., ϱ max = 3. We represent the selection of the per level integration orders as a list of length ϱ max + 2, where the first ϱ max + 1 entries refer to the bisectioning levels 0 to ϱ max . Note that the level ϱ = 0 in fact pertains to cells that are not intersected by the boundary. The final entry in the list of orders pertains to the integration order on the tessellated sub-cells at level ϱ max + 1.

Two-dimensional test case
We consider the two-dimensional test case introduced in Ref. [57], which is illustrated in Fig. 25. The size of the bounding square is set to L = 1, and the radii of the exclusions are taken as R 1 = 0.3 and R 2 = 0.2. The right boundary is displaced byū = 1 2 . Fig. 26 displays the solution to the elasticity problem (69) in terms of the displacement magnitude and the shear stress. The displayed result is based on a 10 × 10 ambient domain mesh with uniformly distributed fourth order Gauss integration scheme. Fig. 27(a) displays the computed constrained modulus for different steps in the integration optimization procedure. It is observed that by selecting a second order Gauss scheme on the untrimmed elements, and a single point in each sub-cell of the trimmed elements, yields a modulus of approximately 0.3484M, which is already within 0.1% of the fully resolved integral result. This computation, however, uses a total of only 1208 points, in comparison to the 6810 points used for the full integral computation. Upon application of the integration optimization procedure the constrained modulus evidently converges to that using full integration. A virtually identical result as the full integration result is obtained with 3654 points. For this result fourth order Gauss quadrature is used on the untrimmed cells and on the first bisection level, and second order Gauss quadrature on all lower levels.

Three-dimensional test case
We consider the elastic deformation of a sintered glass specimen, which is extracted from µCT-scan data as discussed in Ref. [58]. 4 The size of the bounding box is set to L = 1.5, and the right boundary is displaced byū = 1 2 . Fig. 28 displays the solution to the elasticity problem (69) in terms of the displacement magnitude and the σ xy shear stress. The displayed result is based on a 10 × 10 × 10 ambient domain mesh and uniformly distributed fourth order Gauss integration scheme.
As reported in Ref. [7], for coarse meshes the basis functions artificially connect disjoint parts of the geometry. This effect is also observed in the shear stress displayed in Fig. 28(b). For the purpose of illustrating the performance of the integration optimization procedure considered in this work, this effect is, however, not prohibitive.  Fig. 27(b) displays the constrained modulus versus the number of integration points. By using second order quadrature on the untrimmed cells and one point per sub-cell in the trimmed elements, provides a modulus of approximately 0.0925M, which is within 2.5% of the fully resolved integral result. However, this result is obtained at the expense of 1.4 × 10 6 integration points. This is a factor thirteen lower than the 18.4 × 10 6 points used to compute the fully resolved integral result for a reference. A virtually identical result to the reference result is obtained at 6.8 × 10 6 points, which is still approximately a factor of three lower than the full integration scheme.

Conclusion
We have developed an algorithm to construct quadrature rules for cut-elements in which the integration points are distributed over the (cut) element in such a way that the integration error is minimized. Strang's first lemma -which provides an error estimate for the approximate solution including integration errors -provides a theoretical underpinning of the developed integration procedure. In the setting of this approximation theory, it is found that the integration error is bound by the supremum of integration errors over the considered test space. When the integration error is localized to the (cut) elements, this supremum translates to the consideration of a worst possible function to be integrated over a particular element. By discretizing the space of functions over an element by means of polynomials, the worst possible function and its associated integration error become computable.
The ability to evaluate the integration error forms the basis of the proposed integration optimization procedure for octree subdivision. The pivotal idea behind the proposed algorithm is to gradually increase the number of integration points by adding integration points to the sub-cells for which the error reduction per added integration point is highest. As the number of sub-cells increases dramatically with an increase in number of dimensions and with the octree depth -as illustrated by the derived scaling relation for the number of sub-cells -such an optimization procedure has the potential to significantly reduce the computational effort involved in the integration procedures for immersed methods. The presented numerical simulations demonstrate that, for a given number of integration points, the integration error resulting from the optimization algorithm is, in general, significantly lower than that of the integration scheme considering the same integration orders on all sub-cells. Conversely, when fixing the error, the number of integration points required for the optimized quadrature rule is significantly lower than that of the equal order schemes. The considered simulations demonstrate that the developed optimization algorithm efficiently distributes integration points over cut-elements for a wide range of cut-cell configurations.
The developed algorithm provides insight into the way in which integration points should be distributed over the cutelements in order to minimize integration errors. Based on these insights, it is evident that integration rules for which the number of points per sub-cell is decreased with increasing octree depth form a good approximation to the optimized distribution of integration points. The integration point distribution algorithm presented in Ref. [27] therefore can also be expected to yield quadrature rules that are a good approximation to the rules that follow from the theoretical framework considered in this contribution.
To effectively use the quadrature optimization algorithm one must choose various parameters, most importantly the order of the polynomial integrand to be considered and the norm used for the integrand normalization. In particular the order of polynomials is important with respect to the obtained distribution of integration points, in the sense that selecting this order too low will lead to suboptimal integration results. In general it is possible to determine the order of the integrands based on the problem under consideration and the finite element basis functions being considered. The discrete problem setting therefore provides a rational for the selection of the optimization algorithm parameters.
Although the costs involved in the quadrature optimization algorithm are limited on account of the fact that the employed polynomial basis should only be integrated once using the expensive full order integration scheme, the determined optimized integration rules are particularly attractive when they can be re-used for multiple simulations. This is, for example, the case in time-dependent or non-linear problems. The attained reduction in number of integration points is also highly advantageous when data is to be stored in integration points, such as is the case for example with history data in many non-linear constitutive models. Further performance improvements of the presented optimization algorithm can be considered. For example, nested Gauss schemes can improve the computational performance on account of the fact that the integration points are nested in these schemes, which facilitates the re-using of information between steps in the optimization procedure.
In this contribution we have restricted ourselves to optimizing the distribution of integration points over the volumes of the cut-elements. The underlying theoretical framework is, however, more general in the sense that other parameters controlling the integration error can be incorporated as well. In the context of the considered octree integration procedure it would be natural to include the octree depth in the optimization procedure, as this parameter controls the geometry error. Since both the error reduction associated with the increase in octree depth and the associated computational cost are quantifiable, including this depth as an additional parameter in the optimization algorithm is anticipated to be feasible. Inclusion of this parameter is, however, considered beyond the scope of the current work. The same holds for the incorporation of errors associated with the boundary integrals, in particular the Nitsche terms typically considered in finite cell simulations.
The focus in this work has been on the optimization of the integration error contribution in Strang's first lemma, not taking into account the approximation error. It is noted, however, that the need to optimize the integration error depends on the approximation error. If the approximation error is small compared to the overall error, then there is a need to optimize the integration error. However, if the approximation error is the dominating contribution to the overall error, then there is not a strong need to optimize the integration error. Evaluation of the balance between the approximation error and the integration error is an important topic of further study.
The error estimation framework presented in this work is very generic, in the sense that it enables the evaluation of integration errors for a larger class of integration techniques. In order to translate the error estimation concept into an effective optimization strategy, it is required to attain (possibly in an approximate sense) the sensitivity of the error to the numerical parameters of the considered integration technique. In octree-type algorithms, the localization of the error to the sub-cells provides a natural way of attaining these sensitivities. Application of the optimization strategy to alternative sub-cell-based integration techniques is therefore expected to be straightforward. Other integration techniques likely require the development of tools to evaluate the integration error sensitivities.   Fig. A.29(c), and the splitting of the edges is visualized by the edge colors. Note that the (bi-)linear interpolation of the level set function on the sub-cell level will lead to at most one zero level set point per edge, and hence the immersed boundary cannot fold back along a single edge. Folding back along an edge is, however, possible on the background cell level, as unique intersection points can be determined for all octree sub-cells. (d) Along each of the 4 lines connecting the center of the square with its vertices, linearly interpolated zero level set points are determined. For the case considered in Fig. A.29(d), the two additionally obtained zero level set points are indicated by the blue squares on the yellow diagonals. (e) An additional zero level set point is formed by taking the arithmetic mean of the coordinates of the zero level set points along the diagonals. This point, which we refer to as the (approximate) midpoint of the trimmed boundary, is illustrated by the blue circle in the interior of the cell in Fig. A.29(e). The reason for determining the midpoint by the arithmetic mean is that this point must be in the interior of the background cell, which is a property that is of the utmost importance for the robustness of the tessellation procedure. In the case of a planar immersed boundary, the arithmetic mean will yield an exact zero level set point. Relatively poor approximations occur when the immersed boundary has a significant curvature inside the sub-cell, as the arithmetic mean of the diagonal intersections then provides a poor approximation of the zero level set. It is noted, however, that in the context of the finite cell method, at the lowest level of bisectioning, the curvature of the immersed boundary should be small in relation to the sub-cell size, provided that the octree depth is set appropriately.  From this tessellation procedure applied in two dimensions it is evident that it traverses through the dimensions of the problem, in the sense that it first trims the edges (a-c), after which the quadrilateral is trimmed through the extrusion of the edges to the computed mid-point (d-i). This approach directly extends to the three-dimensional case, an illustration of which is presented in Fig. A.30. In Fig. A.30(b) the level set function is evaluated in all vertices of the cube. For each of the six faces of the cube the two-dimensional procedure as described above is applied, as illustrated on the unfolded cubes in Figs. A.30(d)-A.30(e). Finally, the trimmed faces are extruded toward the computed mid-point of the cube, as illustrated in Figs. A.30(g)-A.30(h). The midpoint is again computed by taking the arithmetic mean of the zero level set coordinates computed along all (eight) diagonals connecting the center of the sub-cell to its vertices, with the level set value in the cell center computed by bilinear interpolation of the vertex values. Note that in the three-dimensional case both tetrahedrons and pyramids are created, based on whether or not a face is trimmed. Integration rules are available for both these shapes.