Non-negative moment fitting quadrature for cut finite elements and cells undergoing large deformations

Fictitious domain methods, such as the finite cell method, simplify the discretization of a domain significantly. This is because the mesh does not need to conform to the domain of interest. However, because the mesh generation is simplified, broken cells with discontinuous integrands must be integrated using special quadrature schemes. The moment fitting quadrature is a very efficient scheme for integrating broken cells since the number of integration points generated is much lower as compared to the commonly used adaptive octree scheme. However, standard moment fitting rules can lead to integration points with negative weights. Whereas negative weights might not cause any difficulties when solving linear problems, this can change drastically when considering nonlinear problems such as hyperelasticity or elastoplasticity. Then negative weights can lead to a divergence of the Newton-Raphson method applied within the incremental/iterative procedure of the nonlinear computation. In this paper, we extend the moment fitting method with constraints that ensure the generation of positive weights when solving the moment fitting equations. This can be achieved by employing a so-called non-negative least square solver. The performance of the non-negative moment fitting scheme will be illustrated using different numerical examples in hyperelasticity and elastoplasticity.


Introduction
Over the last few years, interest in fictitious domain methods, such as the finite cell method (FCM) [12,14,38,44], Cut-FEM [7,8], or CutIGA [15], has increased rapidly, mainly because the mesh generation process is simpler than that of the standard finite element method (FEM). In fictitious domain methods, the mesh does not need to conform to the geometry, which results in elements/cells that are cut by the boundary of the domain. Consequently, those broken elements/cells lead to discontinuous integrands that need to be evaluated using more sophisticated numerical integration schemes, since the standard Gauss-Legendre quadrature will not perform well anymore. In this paper, we will focus on the B Wadhah Garhuom wadhah.garhuom@tuhh.de Alexander Düster alexander.duester@tuhh.de 1 Numerical Structural Analysis with Application in Ship Technology (M-10), Hamburg University of Technology, Am Schwarzenberg-Campus 4 (C), Hamburg 21073, Germany numerical integration of broken cells in the finite cell method -which is a combination of the fictitious domain approach with high-order finite elements [12,38].
One possibility for integrating broken cells is the use of the adaptive spacetree decomposition [2,3,12] -quadtree or octree in two or three-dimensional space, respectively. Thereby, the integration domain of each broken cell is subdivided into smaller sub-cells until a user-defined tree-depth level is reached. Afterwards, the Gauss-Legendre quadrature scheme is applied on a cell or sub-cell level excluding the points that are located outside of the physical domain. The approach based on spacetree decomposition is very robust especially for nonlinear applications, and it is also fully automatic and easy to implement. However, it results in a large number of integration points, which renders the method very expensive. In this paper, we will use this approach as a reference solution when compared to the newly proposed methods.
It was shown in [39,40] that the approach based on the spacetree decomposition could be further improved by compressing the sub-cells of each broken cell that are located completely in the physical domain, which can lead to a much lower number of integration points. Another improvement was suggested in [31,32], where the octree is combined with a node relocation and high-order polynomials in a so-called smart octree. In doing so, the geometry can be captured with a smaller number of refinement levels, and consequently, less integration points are generated.
A second approach for integrating broken cells is based on equivalent polynomials, as proposed by [1,51,52]. Thereby, the discontinuous integrand is first transformed into a smooth one with the help of equivalent polynomials. Afterwards, the standard Gauss-Legendre scheme can be applied to integrate the broken cells since the integrand is not discontinuous anymore.
Another, yet very promising approach for integrating broken cells is based on the moment fitting quadrature. Here, for each broken cell a separate quadrature rule is constructed by solving the moment fitting equations [35,48]. The moment fitting was applied in the context of the finite cell method where an adaptive scheme was applied to distribute the integration points only in the physical domain [27]. Afterwards, a linear least square solver is utilized to solve the moment fitting equations. In [10,23], the position of the integration points was chosen as the Gauss-Legendre points allowing them to be located also outside of the physical domain. This results in a well-conditioned linear system of equations of the moment fitting even for high-order functions. Furthermore, in [11,22], Lagrange polynomials through Gauss-Legendre points were used as a basis for the moment fitting quadrature. In doing so, taking advantage of the Kronecker delta property, the solving of a linear system of equations can be avoided. Those versions of the moment fitting quadrature work well for linear problems.
On the other hand, for nonlinear problems, standard moment fitting quadratures turn out to deteriorate the stability of the overall nonlinear solution procedure. This is due to the existence of some integration points with negative weights which can lead to the divergence of the Newton-Raphson method during the incremental/iterative procedure of the nonlinear computation. To overcome this problem, [11,22] suggested an adaptive version of the moment fitting scheme. Here, the moment fitting quadrature is combined with the adaptive octree, where the moment fitting is applied on a cell or a sub-cell level depending on the volume fraction of the broken cell. This approach improves the robustness of the moment fitting quite well. However, the number of generated integration points is still too large as compared to the standard moment fitting. In [25] a non-negative moment fitting quadrature was suggested for integrating one-dimensional functions. Here, the moment fitting quadrature was extended to include an additional constraint to ensure the generation of only positive weights when solving the linear least square problem. This approach was recently further extended in [34] to two-dimensional problems in linear elasticity as well as in small strain elastoplasticity showing promising results.
In this paper, we will extend the work done by Huybrechs [25] and Legrain [34], and apply the non-negative moment fitting approach to three-dimensional problems of solid mechanics. Furthermore, we will test the stability of the proposed approach for a large deformation analysis including geometrical and material nonlinearities in both hyperelasticity and finite J 2 elastoplasticity. As an application, foamed types of materials will be considered since they pose a complex geometry. Furthermore, a new reduced version of the adaptive moment fitting will be introduced. Thereby, the order of the quadrature is reduced by one for every refinement level of the octree, as it was done in [3]. In doing so, the number of integration points can be reduced significantly while not affecting the accuracy too much. We want to emphasize the fact that a reduction of the number of integration points for nonlinear problems is very important since the cost per integration point is much higher than for linear problems and the system has to be recomputed many times during the incremental/iterative procedure. Therefore, reducing the number of integration points will be very beneficial for several reasons. The overhead of the moment fitting will be amortized during the many re-computations of the nonlinear problem. Finally, the non-negative moment fitting and the reduced adaptive moment fitting approaches will be compared in terms of robustness and efficiency (number of integration points) to the existing approaches of the adaptive octree, the moment fitting, and the adaptive moment fitting. This paper is organized as follows: Sec. 2 illustrates the basic idea of the finite cell method. In Sec. 3, the different numerical integration schemes used in this paper will be summarized, together with the proposed non-negative moment fitting approach. In Sec. 4, preliminary investigations will be presented for the different moment fitting approaches. In Sec. 5, detailed investigations of the different moment fitting quadrature schemes will be carried out based on a hyperelastic material model -including a large deformation analysis. In Sec. 6, the different moment fitting quadrature schemes will be investigated for a finite J 2 elastoplastic material model, also including a large deformation analysis. Finally, the paper is concluded in Sec. 7.

The finite cell method
The finite cell method is a combination of the fictitious domain approach with high-order finite elements [49]. It has been successfully applied to a variety of problems in solid mechanics -such as thermoelasticity [56], geometrical nonlinearities [20,45], explicit and implicit elastodynamics [16,26], biomechanics [43,53,55], elastoplasticity [3,22,50], modeling of fracture and crack propagation [24,37,41,42], and microstructured materials [13,30]. The basic idea of the FCM is outlined in Fig. 1 for a two-dimensional case. Therein, the physical domain is embedded into a fictitious domain e \ , resulting in an extended domain e that has a simple shape. The extended domain is then discretized using a Cartesian grid excluding the cells that are completely outside of the physical domain.
To account for the geometry of the real problem, the indicator function α is introduced to the weak form with a value of one inside the physical domain and a value of zero elsewhere. Consequently, the nonlinear weak form defined in the initial configuration reads [4,54] Here, P corresponds to the first Piola-Kirchhoff stress tensor, η denotes the test function, andt defines the applied traction on the boundary N 0 . Dirichlet boundary conditionsū are applied on the boundary D 0 . Equation (1) is highly nonlinear since it involves both geometrical and physical nonlinearities. Applying the directional derivative DG α e · u results in the following linearized form which will serve as a basis for the Newton-Raphson method.
Here, u refers to the displacement increment and A is a fourth order tensor defined as follows where F corresponds to the deformation gradient. Finally, the linearized weak form reads which is then discretized using the finite cell method utilizing hierarchical shape functions of high order. As a result, the global system of equations reads which needs to be solved at every Newton iteration i within every load step k. Here, K i T (u i k ) and G i (u i k ) correspond to the global tangent stiffness matrix and the out of balance vector, respectively, resulting from the assembly process The cell stiffness matrix k c,i and the residual vector g c,i are defined as follows The quantities with the upper index "v" are implemented in Voigt notation. The discrete gradient operator G contains derivatives of the shapes functions, while the matrix N contains the shape functions of the corresponding cell c. The solution u is then updated at every iteration i of the current load increment k The integrals in Equation (7) are discontinuous because of the presence of the indicator function α. Therefore, special quadrature schemes need to be applied to integrate the cells that are cut by the boundary of the geometry, as will be explained in detail in Sec. 3. We will start with the adaptive octree quadrature and then focus on the non-negative moment fitting scheme.

Adaptive octree quadrature
In the adaptive octree quadrature scheme (AOT), each broken cell is partitioned into sub-cells utilizing an octree decomposition to capture the boundary of the domain more accurately until a user-defined tree-depth level k AOT is reached. Afterwards, the standard Gauss-Legendre quadrature is applied to each of the sub-cells, see [2,12]. Furthermore, all points that are located outside of the physical domain are neglected. In Fig. 2 (bottom), a quadtree of a tree-depth level of k AOT = 2 is applied to a two-dimensional plate with a circular hole with a quadrature order of p q = 2n 1 d − 1 = 9. Here, n corresponds to the number of integration points and d refer to the space dimension of the problem.

Moment fitting
In the moment fitting quadrature scheme (MF), a new rule is derived for every broken cell by solving the following moment fitting system of equations: Here, f j defines the basis functions that serve to approximate the function we need to integrate, and x i denotes the position of the integration point i with the corresponding weight w i . Furthermore, n and m correspond to the number of integration points and the number of basis functions used, respectively. Equation (9) can also be expressed in matrix notation as follows or in a compact form Here, A is referred to as the coefficient matrix containing the evaluation of the basis functions f j at the integration points x i . The right-hand side of Equation (10) is referred to as the moments b, which can be computed in different ways. In this paper, we follow the approach suggested by [22], where the moments are integrated utilizing the adaptive octree quadrature as explained in Sec. 3.1. Using polynomials as basis functions, the moment fitting can integrate a polynomial exactly up to an order of p q . Since the coefficient matrix A depends on the position of the integration points, the system is nonlinear w.r.t. the points x i but linear w.r.t. the weights w i . The system can be transferred into a linear one by fixing the position of the points a priori. In [19,22], the positions were chosen according to the Gauss-Legendre points. Setting n = m results in a square coefficient matrix A. At this point, one has the freedom to decide what polynomials to use for f j . In this section, as suggested in [22], we choose Lagrange polynomials through the integration points such that f j (x i ) = δ ji . In doing so, taking advantage of the Kronecker delta property, the coefficient matrix is then reduced to the identity matrix. Thus, the solution of the linear system is trivial and the weights are computed by integrating the moments as follows This moment fitting approach has the advantage of being able to compute the weights without having to solve a linear system at all. Nevertheless, one drawback of such an approach is that the weights are not guaranteed to be positive for all of the integration points. Unfortunately, negative weights can become very critical when solving nonlinear problems since the Newton-Raphson scheme tends to diverge during the incremental/iterative solution procedure, as will be shown in Secs. 5 and 6 .

Adaptive moment fitting
In [11,22], the so-called adaptive moment fitting (AMF) was introduced to improve the stability of the moment fitting quadrature for nonlinear applications. In the AMF, the moment fitting is combined with the adaptive octree, where for every broken cell the volume fraction v of the physical domain v = volume of the physical domain volume of the extended domain × 100% is computed to identify how badly a cell is cut by the boundary of the domain. Consequently, the moment fitting scheme is applied directly to the cell if the volume fraction v is greater than a user-defined threshold. Otherwise, the cell is partitioned into sub-cells applying the octree decomposition. Subsequently, the moment fitting quadrature is utilized on the broken sub-cells depending on their volume fraction, while the standard Gauss-Legendre quadrature is performed for non-broken sub-cells. Table 1 shows the threshold values of the volume fraction related to each refinement level k MF of the moment fitting, which can reach up to a maximum tree-depth of three. If higher tree-depth levels are needed for the resolution of the geometry, this will be considered for the computation of the moments (13). The tree-depth level used for the computation of the moments is chosen as k AOT -k MF . Note that the adaptive octree refinement for the computation of the moments starts from the existing refinement level k MF , see also the example given below.
In this contribution, the number of integration points used for the moment fitting will be reduced instead of just using the "full" number of points of ( p q +1) d for every refinement level k MF , as was done in [22]. To this end, a reduced version of the adaptive moment fitting (RAMF) is proposed to lower the number of integration points without affecting the solution significantly. Here, a number of ( p q + 1) d integration points is only applied if the cell is not refined, i.e., k MF = 0. If the moment fitting is applied to the first refinement level of k MF = 1, the total number of points is reduced to n = ( p q ) d . However, if the resulting number of points on this level is less than ( p q +1 2 + 2) d , the "full" number of points of ( p q + 1) d is utilized. If the moment fitting is applied to the second or third refinement levels, then a set of ( p q − 1) d or ( p q − 2) d points is used, respectively. Nevertheless, the total number of points on both levels is not allowed to be less than ( p q +1 2 + 1) d to ensure a stable integration scheme, as illustrated in Table 1. The RAMF scheme is only applied for quadrature orders larger than three ( p q > 3).
An example that illustrates the reduction of the integration points is given in Fig. 3 (bottom) for an intended quadrature order of p q = 9 and a tree-depth level of k AOT = 4 in two-dimensions for the sake of simplicity. Since the cell is badly broken, it is first refined into four sub-cells. The moment fitting is performed on the sub-cell on the top left with ( p q ) 2 = 81 integration points instead of ( p q + 1) 2 = 100 points, since the sub-cell is located on a tree-depth level of k MF = 1. As the sub-cell on the bottom left is not broken, the standard Gauss-Legendre quadrature is utilized with ( 2 ) 2 = 25 integration points. The two remaining subcells are still badly broken and are therefore further refined. This process continues in the same manner until a maximum refinement level of k MF = 3 is reached. After that, the moment fitting is utilized -no matter how badly the cell is cut. Higher refinement levels than three are considered in the computation of the right-hand side of the moment fitting equations (13). Let us now examine the sub-cell on the top right tr sc . Here, the moment fitting is applied with only ( p q −1) 2 = 64 integration points instead of ( p q +1) 2 = 100 points since they are located on a tree-depth level of k MF = 2. For computing the moments, the adaptive quadtree is utilized with an integration depth of k AOT − k MF = 4 − 2 = 2, as illustrated in Fig. 3 (top).

Non-negative moment fitting
A non-negative moment fitting approach (NNMF) was recently proposed by Legrain [34] to ensure the construction of a quadrature rule with only positive weights. The approach was applied to one-and two-dimensional problems in linear elasticity and small strain elastoplasticity. Thereby, the construction of a non-negative moment fitting quadrature was realized with the help of a non-negative least square solver (NNLS). In this contribution, the approach of Legrain [34] will be extended to three-dimensional problems including geometrical and physical nonlinearities with hyperelastic and finite J 2 elastoplastic material models.
Starting from the moment fitting equations (9), the type of polynomials to be used needs to be defined. Herein contrast to Sec. 3.2, where Lagrange polynomials were used -integrated Legendre polynomials are utilized for the moment fitting basis f j . The right-hand side of the system (the moments b) can easily be computed by integrating the polynomials over the domain with the help of the adaptive octree scheme. The number of moments is set to m = ( p q + 1) d . To construct the coefficient matrix A, the number of integration points n should be chosen greater than the number of moments m, with all points located in the physical domain [6,34]. This results in an under-determined system of equations that can be solved using a linear least square solver (LLS), as it was done in [27]. Nevertheless, utilizing an LLS solver can lead to integration points with negative weights, which is not desirable in a nonlinear computation, since it can affect the stability of the Newton-Raphson method severely.
To overcome this problem, a non-negative least square solver is applied to the moment fitting system of equations. In doing so, Equation (9) is minimized subject to inequality constraints enforcing the weights to be greater or equal to zero (w i ≥ 0) for every integration point Algorithm 1 Non-negative moment fitting quadrature (NNMF) 9: Apply the AOT with tree-depth k NNMF and 10: (n GL ) d Gauss-Legendre points per broken 11: cell distributed in the physical domain 12: Update n = total number of integration points 13: if ((n < (L × m)) and (i < i max )) then 14: if ((n GL = 50) and (k NNMF < k AOT )) then 15: The way how the integration points are distributed in the physical domain plays an important role. In this paper, we take advantage of the adaptive octree and the Gauss-Legendre points. First, we start by applying the AOT with a tree-depth of k NNMF = 0 with a number of r p q +1 2 d integration points per broken cell where r = 1. The factor r and the depth k NNMF are increased until a minimum of (n ≥ L × m) integration points is obtained, where L is set to a number between three and six, for instance. The integration depth k NNMF is not allowed to surpass the depth k AOT specified to integrate the moments, i.e., k NNMF ≤ k AOT . Once the number of integration points is fixed, it is possible to compute the coefficient matrix A, which has a total size of (m ×n). With the moments b and the coefficient matrix A at hand, the non-negative least square solver is applied to the system Aw = b. Next, the norm of the residual is calculated to check for convergence. If the norm of the residual is larger than a user-defined tolerance and the maximum number of iterations is not exceeded, the factor r is increased by one and the previous steps are repeated, as illustrated in Algorithm 1. It is important to start with a sufficiently large set of integration points to ensure the convergence of the NNLS solver with an acceptable accuracy.
Once the NNLS solver is converged, a significant number of points will have weights that are equal to zero. Those weights and their corresponding integration points are filtered out. The remaining number of integration points, to be used later to integrate the corresponding broken cell, are eventually less or equal to the number of moments (n ≤ m), as shown in Sec. 4.

Non-negative least square solver
The non-negative least solver (NNLS) developed by Lawson and Hanson [33] is utilized in this contribution. The method is mainly based on an active set strategy, as can be seen in Algorithm 2. As an input, the algorithm requires the matrix A, the vector b, and a user-defined tolerance ρ. To begin with, the weights w and an auxiliary vector z of size n are initialized. Afterwards, two sets are introduced that contain the indices of the vector w, namely the passive set P and the active set Z. The weights that are indexed in the passive set P are free to have any value. However, they are corrected during the iterations if they become negative. All entries of w that are indexed in the active set Z are set to zero. The dual is introduced. It can be thought of as a Lagrange multiplier that enforces the equality constraints at every integration point. Once the algorithm has converged, we have and Before starting the main loop, the weights are initialized with zeros, and consequently, the active set contains all indices of w, i.e., Z = {1, . . . , n}. On the other hand, the passive set is empty at this point (P = ∅). Next, at the beginning of the main loop, the index j in the active set Z that corresponds to the maximum entry in v is determined and then moved from the active to the passive set P. Afterwards, the inner loop is executed by first extracting the matrix A P from the matrix A, which contains only the columns of A that are indexed in the passive set P, see Algorithm 2 line 11. The restricted linear least square problem is then solved in line 12 for the unknown z P . All values in z that are indexed in the active set are set to zero in line 13. Next, we check whether all values of the solution of the restricted least square problem are greater than zero. If this is the case, the values of the weights are set equal to the auxiliary vector (w = z) before the inner loop is terminated in line 16. Going back to the main loop, the vector v is updated in line 22. Then, a new variable is considered in line 7. The main loop will then continue until either the maximum entry of v becomes less than a tolerance ρ, or the passive set has no more indices left (Z = ∅).
If some of the entries of z P are not greater than zero after solving the restricted system in line 12, the factor β is computed in line 18, resulting in a value in a range between zero and one. Afterwards, the new weights are computed as a linear combination of w and z with the help of β to ensure positive results, see line 19. Finally, each index j which corresponds to a weight equal to zero is then moved from the passive to the active set. If there are any negative weights at this point due to round-off errors, their index is shifted to the active set as well. Once the algorithm is finished, the weights w are returned (containing only positive values).
There are different ways how to solve the restricted least square problem in line 12. In this contribution, the QR decomposition is employed to solve the system, which leads to a good accuracy even for badly conditioned systems. For a more detailed explanation of the NNLS solver, the reader is referred to the work of Lawson and Hanson [33], and also Bro and Jong [5].

Stabilization of the solution in the fictitious domain
Broken cells that obey a small volume fraction deteriorate the condition number of the resulting global equation system. To overcome this problem, the standard Gauss-Legendre quadrature is applied by distributing a set of ( p +1) d integration points in every broken cell, excluding the points located in the physical domain. Afterwards, the weights of these points are multiplied by the small positive value of the indicator function α. Furthermore, a hyperelastic material model is used for those stabilization points. Taking into account these extra integration points (with very small positive weights) helps to improve the condition number of the resulting overall equation system. In this paper, the indicator function in the fictitious domain is set to α = 10 −5 . Fig. 2 (top) shows the additional points that are added to stabilize the solution : Restrict A to the passive set of unknowns 12: Here, the β value satisfies: 0 < β ≤ 1 19:

20:
Move every index j with w j ≤ 0 from the passive set P to the active set Z

Recovery of the gauss-legendre quadrature
In this section, we first test a situation where the hexahedral cell (in d = 3 dimensions) is not cut by the boundary of the domain. In this case, it is common to apply the standard Gauss-Legendre (GL) quadrature to integrate those cells.
Here, however, we will apply the NNMF and the MF quadratures to check whether they can reproduce the same weights as the GL scheme. To this end, we use the following definition to compute the error in the weights The position and the number of the integration points are set to be the same as the Gauss-Legendre quadrature scheme.
Since the cell is not broken, the accuracy of the NNMF and the MF schemes become ( p q = 2n 1 d − 1). To this end, in Fig. 4, the error e w is plotted using different quadrature orders. It can be observed that both the NNMF as well as the

Cell cut by a sphere with radius r = 1.55 mm
In this section, the proposed NNMF and MF schemes will be investigated and compared for the integration of polynomials of order p i . This is done by considering one hexahedral finite cell with a size of (1 × 1 × 1) mm 3 , cut by a sphere described Here, the center coordinates x c , y c , and z c are set to zero corresponding to one corner of the cell, see Fig. 5. In order to compute the moments (right-hand side of Equation (9)), the adaptive octree is applied with an integration depth of k AOT = 7. Furthermore, the adaptive octree is utilized as a reference solution. To this end, the radius of the sphere is set to r = 1.55 mm, which represents a cell that has a rather large volume fraction with v = 99.42%, see Fig. 5 (bottom). Since the volume fraction of the cell is larger than 85%, utilizing the AMF here will result in using the MF scheme directly without any refinement. The quadrature points generated using the NNMF are distributed in the physical domain, as can be seen in Fig. 5 (top) for a quadrature order of p q = 9. Next, the integration points are generated for quadrature orders p q = 1 to p q = 10 utilizing both the NNMF and the MF schemes. In Fig. 6 (bottom), the condition number κ κ = σ max σ min (22) of the coefficient matrix A is plotted against the quadrature order. Here, σ max and σ min refer to the maximum and minimum singular values, respectively. It can be seen that the condition number increases with higher quadrature orders. Since the cell is not cut too severely, the condition number is not excessively large. For the MF, the coefficient matrix A is equal to the identity matrix because the Lagrange polynomials are defined based on the integration points. Next, we investigate the conditioning of the quadrature rule κ q defined as the sum of the absolute values of the weights The condition number κ q can also be normalized by dividing it by the volume V of the domain Here, a value ofκ q = 1 corresponds to a well-conditioned quadrature rule, while a value ofκ q > 1 indicates the presence of negative weights. This is of great importance since the existence of negative weights severely affects the convergence of the Newton-Raphson scheme during the incremental/iterative procedure, as will be shown in Secs. 5 and 6 . It can be seen in Fig. 6 (top) that the condition number of the NNMF is equal to one for all quadrature orders, while it is greater than one in the MF scheme for quadrature orders bigger than three. Subsequently, the relative error e r e r = I AOT − I q I AOT , is investigated for the integration of polynomials of order p i , where I AOT corresponds to the value of the integral based on the adaptive octree, while I q refers to the value of the integral based on the moment fitting schemes. It can be seen in Fig. 7 (bottom) that both methods produce a very low error of about 1e-14 for polynomial orders up to four. The error is then increased to about 1e-12 for the NNMF applying higher polynomial orders because of the increase in the condition number of the coefficient matrix A. In Fig. 7 (top), the number of integration points are plotted for the different quadrature orders. Here, both the NNMF and the MF schemes result in the same number of points.

Cell cut by a sphere with radius r = 0.3 mm
In this section, the radius of the sphere is set to r = 0.3 mmwhich corresponds to a cell that is badly cut, with a volume fraction of only v = 1.41%, see Fig. 8 (bottom). The quadrature points generated using the NNMF are distributed in the physical domain, as can be seen in Fig. 8 (top) for a quadrature order of p q = 9. To study the different approaches, the integration points are set up for quadrature orders p q = 1 to p q = 10 utilizing the different moment fitting quadratures.
In Fig. 9 (bottom), the condition number κ of the coefficient matrix A (see Equation (22)) is plotted for different quadrature orders. For the NNMF, the condition number evidently increases with increasing quadrature order. Since the cell is badly cut, the condition number is much higher as compared to the case where the cell has a large volume fraction. For the MF, AMF, and RAMF schemes, the coefficient matrix A is equal to the identity matrix since the Lagrange polynomials fulfill the Kronecker delta property at the integration points.
Next, we take a look at the conditioning of the quadrature ruleκ q (see Equation (24)). To this end, Fig. 9 (top) reveals that the condition number of the NNMF is equal to one for all quadrature orders, while it is much larger than one in the MF scheme. The AMF and the RAMF schemes produce a much better conditioning as compared to the MF with about κ q = 1.03. That is why they are more robust for nonlinear problems, see Secs. 5 and 6.
Subsequently, we investigate the relative error e r for the task of integrating polynomials at different orders. Here, the MF, AMF, and the RAMF produce similar accuracy of about 1e-13. Although the condition number of A is much larger for such a badly cut cell, the NNMF scheme produces an excellent accuracy of about 1e-12, see Fig. 10 (bottom).
Finally, the total number of integration points is plotted in Fig. 10 (top). From this, it is evident that both the NNMF and the MF schemes produce a much lower number of integration points (with about a factor of ten) as compared to the AMF. The RAMF scheme shows a significant reduction in the quadrature points as the quadrature order increases. The non-negative least square solver in the NNMF generates a lower number of integration points as compared to MF when the cell is cut by a small volume fraction of the domain. It is clear that if the volume of the physical domain in a cell is too small, a lower number of points should be sufficient for integration as compared to the scenario where the volume fraction is larger.

Application to problems in hyperelasticity
In this section, a detailed investigation of the different moment fitting quadrature schemes is carried out, taking large deformations into account. To this end, a hyperelastic material model is considered with the following strain energy density function [9]  Here, C refers to the right Cauchy-Green deformation tensor and J = √ det(C). The material properties are set to λ = 28.846 N/mm 2 and μ = 19.231 N/mm 2 .
Because of the complexity of the material model, the integrands in Equation (7) are not based on polynomials only. Also, in elastoplasticity, the integrands are discontinuous as the plasticity evolves through the elements.

Plate with a cylindrical hole
In this section, we consider a plate with a cylindrical hole. The goal is to compare the performance and accuracy of the different moment fitting quadrature schemes compared to the adaptive octree as well as to a reference solution based on the p-FEM. The geometry of the plate, depicted in Fig. 11, has a side length of a = 100 mm, a thickness of t = 10 mm, and a radius of r = 60 mm. Symmetry boundary conditions are applied together with a prescribed displacement ofū z = 200 mm at the top surface, leading to a stretch of the plate. For Fig. 11 Plate with a cylindrical hole. Geometry and boundary conditions this, a displacement increment ofū z = 0.5 mm is applied, which requires a total of 400 load steps.
In order to create a reference solution, the geometry is discretized with 3200 curved elements to generate an overkill solution based on the p-FEM, see Fig. 12 (bottom). For the FCM, a mesh of 78 cells (10×1×10 subdivisions) is considered, see Fig. 12 (top). The polynomial degree of the shape functions is set to p = 4. For the adaptive octree in the FCM, an integration depth of k AOT = 4 is used. The Cauchy stress σ zz will be computed at point A with coordinates (x = 35.0, y = 5.0, z = 4.0). Table 2 lists the number of quadrature points generated for the different integration schemes used in this paper. Here, we will study the reduction η and the ratio γ of the total number of integration points compared to the AOT as follows where n AOT corresponds to the total number of integration points generated by the AOT, while n q refers to the total number of integration points produced by the other quadrature schemes. The AOT quadrature produces a total number of 402 915 integration points (IPs) utilizing ( p + 1) 3 IPs in every subcell. In the AMF, we will apply and investigate different numbers of quadrature points to see how the accuracy and stability is affected. Thereby, we start first off with the "full" set of (2 p + 1) 3 integration points resulting in a total number of 378 938 IPs. Here, the number of integration points is reduced by only 5.95% as compared to the AOT. Additionally, we investigate the AMF utilizing a number of ( p + 3) 3 IPs as well as ( p + 2) 3 IPs on each broken cell or sub-cell. In doing so, a reduction of η = 53.6% and η = 69.2%, respectively, can be achieved. Those two versions of the AMF use the same set of points (either ( p + 3) 3 or ( p + 2) 3 ) independently of the refinement level of the moment fitting k MF .
Furthermore, we introduce a reduced version of the adaptive moment fitting (RAMF), where a number of (2 p + 1) 3 IPs is applied only if the cell is not refined. Then, for every new tree-depth level, the number of points is reduced by one in each direction so that the sub-cells on the highest level have the least number of points, as explained in Sec. 3.3. In doing so, we obtain a total number of 130 774 IPs, which is a reduction of around η = 67.5%. The MF scheme pro-  Fig. 13 (bottom), the strain energy is plotted versus the displacement at every load step. It can be seen that all quadrature schemes produce similar values. However, the simulation based on the MF scheme failed to converge at an early load step -already at a displacement ofū z = 8 mmbecause of the presence of negative weights that lead to a high condition number of the quadrature. Additionally, the stress σ zz at point A is plotted in Fig. 13 (top) for the different quadrature schemes. Here, a good agreement between the different methods can be observed as well.
Next, we take a look at the relative error in the strain energy as well as the relative error in the stress σ zz compared to the AOT, see Fig. 14. Here, differences in the relative error can be observed between the individual quadrature schemes. The AMF with a number of (2 p + 1) 3 IPs produces the lowest relative error as compared to the AOT. However, with a ratio of only about γ = 1.06, the number of integration points generated is too large. It is interesting to observe the great improvement in the RAMF -which achieves almost the same accuracy as the AMF while leading to a significant reduction in the number of integration points, with a ratio of γ = 3.08. Here, the MF is the least stable scheme, with the largest relative error compared to the AOT. Nevertheless, the error is still within an acceptable range of about 1e-10 for the strain energy and about 1e-8 for the stress. The reduction in the integration points for the MF is remarkable -with a ratio of γ = 22.89. Finally, the NNMF scheme heals the stability problem of the MF and produces an even lower number of integration points, with a much smaller error compared to the MF. In this case, the NNMF shows a good accuracy as well as a huge reduction in the number of integration points, with a ratio of γ = 24.89. Subsequently, we investigate the relative error in the strain energy as well as in the stress σ zz as compared to an overkill solution from a conforming mesh. Here, all methods show good agreement when compared to a reference solution, see Fig. 15. This is because when compared to the overkill solution, other types of errors play a role as well -such as, for instance, discretization error, introducing a soft material in the fictitious domain, and applying an integration depth of only four (integration error). Since those errors are more dominant, the error introduced from the different integration schemes presented earlier is considered negligible compared to an overkill solution. In the end, it is important to have a stable integration scheme that performs well for nonlinear applications and can generate a low number of integration points for broken cells. Fig. 16 shows a contour plot of the Cauchy stress for the plate with a cylindrical hole at the last load step discretized with the FCM.

Single pore of a foam
In this section, we consider a more complex structure of a single pore of a foam. The geometry is obtained from a CTscan and then converted into a triangulated surface [21], as shown in Fig. 17 (bottom). The bottom surface of the foam is fixed in all directions, while the top surface is fixed in x and z-directions only. Furthermore, a displacement ofū y = 1.5 mm is applied on the top surface to compress the foam. To achieve the final deformation, a displacement increment ofū y = 0.05 mm is used, resulting in a total of 30 load steps. The geometry is discretized using 2721 finite cells (25 × 25 × 25 subdivisions) with a polynomial degree of the shape functions of p = 2, see Fig. 17 (top). Table 3 shows the number of quadrature points created by the various integration strategies. Using a tree-depth level of k AOT = 4 and a number of ( p + 1) 3 IPs in each subcell, the AOT quadrature generates a total of 21 138 178 IPs. The AMF produces a total of 13 072 543 IPs utilizing (2 p + 1) 3 integration points in each broken cell or sub-cell. In comparison to the AOT, the number of integration points is reduced by only η = 38.2%. Further, we study the reduction in the number of points in the AMF by the same amount in each refinement level. In doing so, a number of ( p + 2) 3 and ( p + 1) 3 IPs are considered, resulting in a reduction of η = 67.9% and η = 85.9%, respectively.
Furthermore, by applying the RAMF, we obtain a total of 6 957 842 IPs, which is a reduction of η = 67.1%. The MF scheme creates a significantly lower number of quadrature points, with only 319 153 IPs (or around 98.5% fewer integration points). Finally, the NNMF quadrature generates the lowest number of points: only 283 180 IPs, which is a total reduction of η = 98.7% as compared to the AOT. As indicated in Sec. 3.5, a total of 41 869 IPs are assigned in the fictitious domain for stabilization purposes in all of the quadrature schemes discussed.
The strain energy is plotted against the displacement at each load step in Fig. 18 (bottom). As can be observed, all quadrature techniques generate comparable energy values. However, owing to the existence of negative weights, the MF scheme failed at an early load step with a displacement ofū y = 0.25 mm. Furthermore, the AMF simulation with ( p + 1) 3 IPs could only be carried out up to a displacement ofū y = 1.15 mm. Fig. 18 (top) shows the relative error in strain energy compared to the AOT. Differences in accuracy between the various quadrature techniques can be seen here. The AMF with (2 p + 1) 3 IPs achieves high accuracy, but the number of integration points generated is still too large, with a ratio of only γ = 1.62. The RAMF, on the other hand, shows a great accuracy compared to the AMF, with a significant reduction Table 3 Total number of integration points for the different quadrature schemes with tree-depth k AOT = 4 and shape function order p = 2 of the single pore of a foam  in integration points and a ratio of γ = 3.04. The AMF with ( p + 1) 3 IPs results in the largest error as compared to the AOT, failing to converge at an earlier load step. Finally, with a ratio of γ = 74.65, the NNMF demonstrates an excellent reduction in the number of integration points while keeping a very good accuracy as well. Subsequently, Table 4 lists the simulation and quadrature setup computation times for the various integration schemes. The simulations are performed on a machine with two CPUs, each with ten cores (40 threads) and a clock speed of 2.4 GHz. Although the time required to set up the quadrature points using the AOT was the shortest (about 2 minutes), the simula- Fig. 19 The von Mises stress of the single pore of a foam at the last load step using p = 2 tion computation time was the highest one (almost 2 hours). The AMF using (2 p+1) 3 required around 6 minutes to generate the points, followed by 1 hour and 21 minutes to complete the simulation, which is roughly 40 minutes less than for the AOT. We see a good saving in computation time when utilizing the RAMF scheme, with only 46 minutes needed to finish the simulation. The NNMF scheme outperformed the other quadrature methods in computational time, requiring just 3 minutes to create the integration points and around 7 minutes to complete the entire simulation. Finally, Fig. 19 shows the contour plots of the von Mises stress for a single pore of a foam at the last load step. Here, higher values of the von Mises stress can be noticed at the foam struts.

Application to problems in elastoplasticity
In this section, a detailed investigation of the different moment fitting quadrature schemes will be carried out, taking large deformations into account with a finite J 2 elastoplastic material model. First, a short summary of the governing equations for the material model will be given. To this end, the elastic part of the deformation is described based on a hyperelastic neo-Hooke model with the following strain energy density function Adaptive moment fitting (2 p + 1) 3 1h , 21min 5min , 48s Reduced adaptive moment fitting 46min , 51s 5min , 28s Non-negative moment fitting (2 p + 1) 3 6min , 57s 3min , 01s where I b e = tr b e and I I I b e = det b e correspond to the first and third invariant of the elastic left Cauchy-Green tensor b e , respectively. The von Mises yield criterion is defined as Here, τ refers to the Kirchhoff stress tensor with its corresponding deviatoric part s. Furthermore, a nonlinear isotropic hardening is considered which is defined as follows whereᾱ defines the equivalent plastic strain. All material parameters are listed in Table 5. For more details about the material model, the reader is referred to [28,29,46,47].

Plate with a cylindrical hole
In this section, we once again consider the plate with a cylindrical hole as described in (Sec. The displacement is then increased by 0.5 mm at every load step until the complete deformation is reached, which requires 47 load steps. Table 6 shows the total number of quadrature points created by the various integration schemes. Using a tree-depth of k AOT = 4 and a number of ( p + 1) 3 IPs for each sub-cell, the AOT quadrature generates a total of 710 265 integration points. Starting with (2 p + 1) 3 integration points in each refinement level, the AMF produces a total number of 485 979 IPs. In comparison to the AOT, the number of integration points is reduced by η = 31.6%. Furthermore, we apply the AMF with a number of ( p +3) 3 and ( p +2) 3 IPs in each broken cell or sub-cell. This results in a decrease of η = 65.6% and η = 76.8%, respectively. Moreover, we obtain a total of 178 243 IPs using the RAMF, which is around 74.9% less points than the AOT. The MF scheme yields a significantly lower number of points, with just 35 476 IPs, or 95% fewer integration points. Finally, the NNMF quadrature generates the lowest number of points, with only 32 285 IPs, with a total reduction of η = 95.5%. A total of 1060 IPs are distributed in the fictitious domain for stabilization purposes in all of the quadrature schemes mentioned above, as explained in Sec. 3.5.
The load displacement curves are plotted in Fig. 20 (bottom) at each load step. It can be observed that all quadrature schemes are in good agreement with the reference solution obtained from a conforming curved mesh (see Fig. 12 (bottom)). However, at a displacement of merelyū z = −0.1 mm, the simulation based on the MF technique failed to converge very quickly. Furthermore, the computation based on the AOT scheme failed to reach the last load step, resulting in a deformation ofū z = −16.5 mm. The final load step ofū z = −20 mm was achieved by all of the remaining moment fitting quadratures. Additionally, the relative error in the reaction force is plotted in Fig. 20 (top). In this case, there is a relatively good agreement between the various quadrature rules as compared to the reference solution. Finally, the von Mises stress σ v M and the equivalent plastic strainᾱ are plotted in Fig. 21 at the last load step. It is clear from this that the plastic zone begins to evolve in a diagonal direction across the plate.

Single pore of a foam
In our last application, we once again consider the pore of a foam described in Sec. 5.2, but now with a finite J 2 elastoplastic material model. A prescribed displacement ofū y = 1.0 mm is applied to compress the foam with a polynomial degree of the shape functions of p = 3. Small load steps are applied, beginning withū y = {0.0001, 0.001, 0.003, 0.005} mm. After that, the displacement is raised by 0.01 mm at a time until the full deformation is obtained, which takes a total of 104 load steps. Table 7 lists the total number of quadrature points generated for each integration scheme. Using a tree-depth of k AOT = 4 and a number of ( p + 1) 3 IPs on each sub-cell, the AOT quadrature gives a total number of 50 142 769 IPs. We apply and analyze different numbers of quadrature points using the AMF, just as we did in the previous sections. The AMF generates a total of 35 793 794 IPs with (2 p + 1) 3 integration points in every broken cell or sub-cell, yielding a decrease of only 28.6% as compared to the AOT. In addition, we apply a set of ( p + 3) 3 and ( p + 2) 3 IPs in every broken cell or sub-cell of the AMF, resulting in η = 54.7% and η = 73.4% reductions, respectively, as compared to the AOT. Furthermore, if the RAMF is applied, a total of 13 656 014 IPs is obtained, which is a 72.8% reduction as compared to the AOT. The MF and NNMF schemes produce the lowest number of integration points with only 873 597 and 675 420 IPs, respectively, which is a reduction of about η = 98%. Also, a total of 99 181 IPs are distributed in the fictitious domain for stabilization purposes, as explained in Sec. 3.5. In Fig. 22 (bottom), the load displacement curves are plotted utilizing the different integration methods. It can be observed that all integration schemes show a good agreement in the reaction forces. However, the simulation based on the MF scheme failed to converge at an early load increment of u y = 0.003 mm. Also, the simulations based on the AMF with ( p + 3) 3 and ( p + 2) 3 IPs only achieved a deformation state ofū y = 0.61 mm andū y = 0.68 mm, respectively. Here, the NNMF and the RAMF are the only two schemes that were able to reach the last load step ofū y = 1.0 mm. The AOT and the AMF with (2 p +1) 3 reached a deformation state ofū y = 0.92 mm. Additionally, Fig. 22 (top) shows the relative error in the reaction force at different load steps as compared to the adaptive octree. It can be noticed that all integration schemes produce approximately a relative error of about 1e-5. Next, the computation time of the simulation and the quadrature setup are listed in Table 8 for a displacement of u y = 0.92 mm. The simulations are run on a computer with two CPUs, each with ten cores (40 threads) and a 2.4 GHz clock speed. Although it took the least amount of time to set up the quadrature points using the AOT (about 2 minutes), the computation time for the entire simulation was the most costly one (approximately 14.7 hours). The AMF took around 6 minutes to generate the points and 10.7 hours to conduct the entire simulation. The RAMF scheme reduced the computation time to just 4.57 hours. The NNMF scheme needed around 5 minutes to generate the integration points, and then used about 0.63 hours to finish the whole simulation. The NNLS algorithm can be terminated after a certain number of iterations to allow for a fast quadrature setup while not affecting the quality of the solution largely. In Fig. 22, we observe a high accuracy of the NNMF scheme although the number of iterations is cut earlier to reduce the setup time of the quadrature.
It is remarkable to see how much computation time can be saved by utilizing the moment fitting quadrature for nonlinear applications. Generally, the time spent to generate the points and weights in the moment fitting is amortized by reusing the quadrature rule numerous times during the incremental/iterative procedure of the nonlinear computation. Finally, the von Mises stress σ v M and the equivalent plastic strainᾱ are plotted in Fig. 23 at the load step ofū y = 1.0 mm. The plastic zone starts to evolve from the struts of the pore of the foam, attaining greater values of the equivalent plastic strain and the von Mises stress.

Conclusions
In this paper, we proposed a non-negative moment fitting quadrature scheme to integrate broken cells using the finite cell method. Thereby, all integration points are constructed to have only positive weights. This is possible by utilizing a non-negative least square solver, where the moment fitting equations are solved including inequality constraints to make sure that the weights are greater or equal to zero. To ensure high accuracy of the weights, the number of integration points should be chosen larger than the number of moments, which results in an under-determined linear system of equations. The Lawson and Hanson non-negative least square solver [33] produces a large number of zero-weight integration points, which are subsequently removed from the set of points to be utilized for integrating the broken cells. Consequently, we showed that the resulting number of integration points is eventually less or equal to the number of moments. In this work, the adaptive octree was used for the moment fitting to accomplish two goals: first, to capture the domain of interest and distribute the integration points exclusively in the physical domain for the non-negative moment fittingand second, to integrate the right-hand side of the moment fitting equations. Furthermore, we proposed a reduced version of the adaptive moment fitting, where the number of quadrature points is reduced for the higher refinement levels.
To this end, we investigated different numerical examples including large deformation analysis in both hyperelasticity and elastoplasticity. We showed that the reduced adaptive moment fitting scheme reduced the number of integration points significantly while keeping the same accuracy and stability as observed for the adaptive moment fitting. In addition, we showed that the non-negative moment fitting scheme is very robust for nonlinear applications, generating the lowest amount of integration points with high accuracy as compared to the other integration methods.

Future outlook
The computation time for the setup of the non-negative moment fitting scheme can be divided into two main parts. The first part is related to the computation of the right-hand side of the moment fitting equations. In this paper, the adaptive octree was utilized to integrate the moments, which can become very expensive when using high tree-depth levels to capture the geometry. Since the moments are based on the integral of polynomials, they can be integrated a priori and stored in a tabulated form, see [1]. After that, computing the moments will simply involve reading the results from a table, leading to a significant reduction in computation time. The second part is related to solving the moment fitting system of equations. Within the non-negative least square algorithm, a restricted linear system of equations is solved in several iterations, where at every iteration one integration point is added to the system or removed from the system. In this paper, a QR decomposition is employed to solve the restricted system on a single thread. To further speed up the non-negative least square solver, it is possible to apply some of the high-performance and well-parallelized linear algebra libraries -such as, for example, OpenBLAS and Intel MKL. Furthermore, in this paper, we applied the non-negative least square algorithm from Lawson and Hanson [33], which is based on an active set strategy. It would also be possible to consider more advanced solvers -such as the fast nonnegative least square algorithm (FNNLS) [5] or the recently developed TNT-NN algorithm [36] -for which the authors claimed that it is several times faster as compared to the nonnegative least square solver from [33].
To improve the accuracy of the elastoplastic simulations, a hp-adaptivity can be applied to capture the plastic front dynamically [58,59]. Alternatively, the polynomials that define the moment fitting basis can be enriched to account for the discontinuity introduced by the elastoplastic front [10].