A posteriori error control for the finite cell method

The paper presents some concepts of the finite cell method and discusses a posteriori error control for this approach. The focus is on the application of the dual weighted residual approach (DWR), which enables the control of the error with respect to a user‐defined quantity of interest. Since both the discretization error and the quadrature error are estimated, the application of the DWR approach provides an adaptive strategy which equilibrates the error contributions resulting from discretization and quadrature. The strategy consists in refining either the finite cell mesh or its associated quadrature mesh. Numerical experiments confirm the performance of the error control and the adaptive scheme for a non‐linear problem in 2D.


Introduction
The Finite Cell Method (FCM) is designed for the simulation of the mechanical behavior of a body by using an unfitted mesh and enables computations on extremely complex shapes. The FCM was originally introduced in [1] as an immersed boundary method and has been successfully applied to a variety of problems in solid mechanics, cf., e.g., [2,3]. Its basic idea is to combine a fictitious domain approach with a finite element method, where the possibly complicated physical domain is replaced by an embedding domain of a geometrically simple shape, for instance, a (paraxial) quadrilateral in 2D or a (paraxial) hexahedron in 3D, which can easily be meshed. The variational formulation of the problem and its finite element discretization are defined on the embedding domain. The geometry of the physical domain is incorporated via an indicator function which necessitates the implementation of an appropriate quadrature scheme. This can be done by using an additional quadrature mesh, which approximates the physical domain and, thus, implies a certain quadrature error. Hence, mainly two computational error sources occur in the finite cell method: the discretization error and the quadrature error.
In this paper, we present some concepts and applications of the finite cell method and discuss a posteriori error control for this method. The focus is on the application of the dual weighted residual approach (DWR) as presented in [5]. It enables the control of the error with respect to a user-defined quantity of interest. The aim is to estimate both the discretization error and the quadrature error. The application of the DWR approach provides an adaptive strategy in which the error contributions resulting from discretization and quadrature are balanced. The strategy consists in refining either the finite cell mesh or its associated quadrature mesh. In several numerical experiments the performance of the error control and the adaptive scheme is demonstrated for a non-linear problem in 2D.

Finite cell method
In this section, we briefly introduce a general abstract setting of the finite cell method (FCM). Typically, one seeks a solution u ∈ V of an operator equation A(u) = 0, where the operator A : V → V * maps each element of a suitable function space V into its dual space V * . A prototypical example for such an operator is A(u)(v) := Ω ∇u ∇v dx − Ω f v dx encoding Poisson's problem in its weak formulation. Here, the operator is defined on the Sobolev space V = H 1 0,Γ D (Ω) with Dirichlet zero data on the part Γ D of the Lipschitz domain Ω, see Fig. 1a. The application of the FCM to such a problem consists of two steps. First, an enclosing domain Ω ⊃ Ω is chosen which allows for a simple discretization, e.g., a rectangle in 2D or a cuboid in 3D, see  second step consists in defining a finite-element space V h ⊂ V on a mesh T of Ω. Due to the choice of Ω, the meshing process is simplified substantially as, e.g., a paraxial mesh may be chosen, see Fig. 1c. The integrals of functions defined on Ω arising in the definition of the operator A cannot be computed as usual since no triangulation of Ω is available. Instead, the possibly complicated geometry of elements cut by the boundary needs to be captured by the numerical integration. For this, several techniques have been developed. A first class of techniques aims at exact integration, i.e., the error attains the level of machine precision. One drawback of these techniques is the possibly high computational effort. Moreover, it is not necessary to perform exact integration if the overall error is dominated by the approximation error. Hence, it is often more convenient to adapt the quadrature such that the quadrature error is at a similar level as the approximation error. Such an adaptation of the quadrature necessitates the use of refinable quadrature schemes which comprise the second class of widelyused FCM quadrature techniques. Here, a quadrature level : T → N 0 is defined on the mesh which translates into individual approximations of each T ∩ Ω. Their union yields an approximate domain Ω ( ) ≈ Ω. These schemes are usually constructed such that an increase in (T ) delivers an improved approximation of the geometry of T ∩ Ω. Prototypical examples of such a refinable quadrature scheme are the quadtree in 2D and the octree in 3D. Here, the complicated geometry of T ∩ Ω is captured by (T ) recursive refinements of T towards the boundary of Ω passing through T , see Fig. 2a. Note that the quadtree delivers only a piecewise constant approximation of ∂Ω, which implies that a large level is required to sufficiently resolve the boundary. It is therefore practical to use a piecewise linear approximation to ∂Ω as shown in Fig. 2b. By using inexact integration, the original operator A is replaced by its perturbation A ( ) where each integral in the definition of A is replaced by the quadrature rule based on Ω ( ) . Thus, in the application of the FCM, a perturbed discrete solution u The most important feature of the FCM is to resolve geometries with very complicated or fine structures at the boundary. Fig. 3a shows an example of such a geometry: a large mountain region with a complex and fine-structured surface (the Lauterbrunnen Valley in Switzerland). Such a terrain's surface is typically represented by a digital elevation model (DEM) as a heightmap. The application of the marching volume polytopes algorithm as introduced in [7] gives a piecewise linear approximation of the surface, see Fig. 3b. The resulting volume mesh can directly be used in the FCM, for instance, to compute stress distributions in order to detect potential sites of rock failure and to deduce erosive processes in mountain peaks, rock slopes and bedrocks. The FCM approach can also be applied to geometries with varying boundaries (in time). For instance, in [6] the varying boundary in an NC-milling process is resolved by using the marching volume polytopes algorithm in order to simulate thermoelastic deformations. We emphasize that using FCM drastically reduce the engineering effort for pre-processing (almost no cost for meshing) and one benefits from the high robustness of the method with respect to imprecise or even flawed CAD models [8]. Modeling flaws occur frequently in B-rep models used in everyday engineering practice. An example of a B-Rep model is given in Fig. 4a. It represents a steel joint and was supplied by ITE@Braunschweig. The B-Rep model is submerged into a non-boundary conforming computational grid and the FCM is applied to compute the spatial distribution of the von Mises stresses depicted in Fig. 4b. A hierarchical mesh refinement of two levels is applied and an octree integration is used as a refinable quadrature scheme with a fixed depth of (T ) = 3. The model contains only 3 360 degrees of freedom and the corresponding system can be solved using a direct solver in only 129 ms on a standard personal computer. However, the octree quadrature of the system matrices takes approx. two minutes. That is, one thousand times more time is spent on integration than on solving. This indicates the potential gain by including an estimation of the integration error (as proposed in the next section) when solving problems of industrial relevance.

A posteriori error control
While classical error control often aims for estimates in the energy or the L 2 error, we focus on goal-oriented error control. Here, estimates for the error |J(u) − J(u h )| measured with respect to a goal functional J : V → R are sought. The goal functional encodes a user-defined quantity of interest such as point evaluations, averages over subsets, line integrals, or nonlinear quantities such as |∇u(p)| 2 for p ∈ Ω. In the setting of the FCM, error contributions resulting not only from the discretization, but also from the quadrature need to be considered. In this section, an a posteriori error estimate allowing for a separation of these two error sources is presented. To this end, the dual weighted residual method (DWR) as introduced in [4] is extended to the FCM. The main idea is to represent the goal functional J and its approximation J ( ) by means of solutions z ∈ V and z ( ) The following error representation is derived in [5] and contains the evaluation of the approximated operators as well as the exact operators. The evaluation of the latter involves integrals on Ω which cannot be calculated exactly. To overcome this problem, a refineable quadrature scheme : T → R is used. By performing a sufficiently large number k of additional refinements on each element, the error incurred by the quadrature is negligibly small. The operators whose integrals are based on the refinable quadrature schemes are indexed by and + k.
Theorem 3.1 Let u + , z + ∈ V + be approximations of u and z in a higher-order approximation space h is a higher-order remainder.
The terms e in terms of the primal and dual residuals ρ and ρ * given by The notation (y) 0 for a function y onΩ denotes the extension of the restriction y| Ω by zero ontoΩ. The error contributions e N S,· mainly consist of differences of evaluations between the exact operators and their approximations based on improved quadrature rules and can, thus, be assumed to be negligibly small. Abbreviating the primal error e := u − u

Numerical results
In this section, we present the results of two numerical examples from [5]. In the first one, the claim that the terms e N S,· in the error representation are negligibly small compared to the dominating error terms e D and e Q is validated. We consider the semilinear problem −∆u + sin(πx) sin(πy) and z(r, ϕ) := (−r 2 + 2r)(− 16 π 2 ϕ 2 + 8 π ϕ). We use a piecewise linear quadrature scheme and set the quadrature refinement level = 2 as well as the additional level k = 2 and perform uniform mesh refinement of lowest-order finite elements. In Fig. 5a, the decrease of error terms is O(h 2 ) as expected. Moreover, the terms e N S,· are indeed at least one order of magnitude smaller than e D , e Q .
In the second example, automatic mesh adaptation is performed. The usual Solve-Estimate-Mark-Refine loop needs to be modified to account for the adaptation of the quadrature mesh, see Alg. Require: Coarse level (T ) = r ∈ N 0 , T ∈ T , additional quadrature levels k ∈ N, error ratio σ ∈ (0, 1) repeat Setup the quadrature meshes according to , + k Solve the primal and dual problem with the quadrature mesh based on and compute e D , e Q with level + k if e Q > σ · e D then Localize e Q and increase for mesh elements with highest contributions else Localize e D (e.g. by filtering, see [5] and references therein) Mark mesh elements with highest contributions, refine T , update and reduce by 1 on new mesh elements end if until termination criterion is met ∩ Ω is considered. The results comparing uniform refinement with two variants of adaptive refinement are displayed in Fig. 5b, where the parameter σ in Alg. 1 is chosen as σ := 0.01. The variant "adaptive 1" creates a quadrature rule based on the quadtree, while "adaptive 2" performs a piecewise linear approximation of ∂Ω. The quadtree computation stops at approx. N = 10 000 due to the large number of quadrature points, which demonstrates the superiority of the linear approximation. While the adaptive refinement yields an optimal algebraic convergence of O(N −1 ), uniform refinement does not attain the optimal rate. In Fig. 6, adaptive mesh refinements towards the reentrant corner and the point p can be seen. The quadrature levels after 12 steps of the refinement loop are displayed in Fig. 7. In Fig. 8, the effect of the adaptive strategy is visualized. Whenever e Q is too large, only the quadrature mesh is refined while the finite element mesh remains unchanged. This leads to a reduction of the quadrature error at the same N . The overestimation of the error by the estimator, eff, and the overestimation of the error by the sum of the absolute values of the error indicators, ind, is displayed in Fig. 9: These numbers are bounded between 1.5 and 8 for meaningful problem sizes and, therefore, indicate a good representation of the error by the estimator. Finally, Fig. 10 displays the linear increase in the total number of quadrature points for the piecewise linear quadrature and the range of the quadrature level per computation step, which shows a local adaptation of the quadrature. We refer to [5] where further numerical experiments are presented, in particular, in the context of model plasticity with linear isotropic hardining. Moreover, implementation aspects as, for instance, the computation of u + and z + and localization techniques are discussed in more detail.