Reliable Residual-Based Error Estimation for the Finite Cell Method

In this work, the reliability of a residual-based error estimator for the Finite Cell method is established. The error estimator is suitable for the application of hp-adaptive finite elements and allows for Neumann boundary conditions on curved boundaries. The reliability proof of the error estimator relies on standard arguments of residual-based a posteriori error control, but includes several modifications with respect to the error contributions associated with the volume residuals as well as the jumps across inner edges and Neumann boundary parts. Important ingredients of the proof are Stein’s extension theorem and a modified trace theorem which estimates the norm of the trace on (curved) boundary parts in terms of the local mesh size and polynomial degree. The efficiency of the error estimator is also considered by discussing an artificial example which yields an efficiency index depending on the mesh-family parameter h. Numerical experiments on more realistic domains, however, suggest global efficiency with the occurrence of a large overestimation on only few cut elements. In the experiments the reliability of the error estimator is demonstrated for h- and p-uniform as well as for hp-geometric and h-adaptive refinements driven by the error estimator. The practical applicability of the error estimator is also studied for a 3D problem with a non-smooth solution.


Introduction
The Finite Cell method (FCM) [12,30] is a popular immersed-boundary method that combines the fictitious domain approach [16,32] with (higher-order) finite elements. Its basic idea consists in embedding the possibly complicated physical domain of the problem into an enclosing domain of a simple shape, which can easily be meshed. The geometry of the physical domain is recovered by multiplying the finite element functions defined on the enclosing domain with an indicator function which is 1 in the domain and a small value in the fictitious domain (the difference of the enclosing and physical domain). Obviously, this proceeding shifts the problem of meshing a complicated domain to the problem of integrating restrictions of functions by means of an appropriate quadrature rule. Since the FCM enables computations on extremely complex geometries, it has been applied to a vast number of involved problems including thermo-elasticity [40], geometrical non-linearities [33], bio-mechanics [31,38], elasto-plasticity [1,37], and foamed materials [18,19].
While a posteriori error control is available for many standard finite element approaches, only few publications deal with a posteriori error control for the FCM or similar cut finite element methods [2,4]. These include recovery-based error estimates [28], goal-oriented error estimates [8,10,11,14,15,29] and implicit error estimates [35]. An error estimate for cut cell finite volume methods has been derived in [13]. To our knowledge, a residual-based error estimate has only been derived for a so-called composite finite element method in [7]. However, these results are limited to domains with small holes, lowest-degree elements and mixed boundary conditions with homogeneous Neumann boundary conditions. A posteriori estimates for some immersed-finite element methods incorporating the domain-approximation error are provided in [5,17].
In this work, we rigorously establish the reliability of a residual-based error estimator for Poisson's problem in the setting of the FCM. Our theoretical findings are suitable for hpadaptive finite elements and deal with mixed boundary conditions, where Dirichlet boundary conditions are assumed to be homogeneous and strongly imposed on boundaries being compatible with edges of the underlying finite element mesh. In contrast, the Neumann boundary conditions are allowed for curved boundaries intersecting the interior of the enclosing domain. We refer to [5] for the application of Nitsche's method which enables weakly imposed Dirichlet conditions on boundaries in the interior of the enclosing domain. The quadrature error is supposed to be negligible; in particular, the domain-approximation error is assumed to be negligibly small. Further investigations would be needed to account for weak Dirichlet conditions as well as the domain-approximation and the quadrature error, which are not taken into consideration in this work. We refer to [11] where the quadrature error is considered in the context of a posteriori error estimates based on the dual-weighted residual approach. The reliability proof of the error estimator relies on standard arguments of residual-based a posteriori error control as, for instance, introduced in [27] for hp-finite elements, but includes several modifications with respect to the error contributions associated with the volume residuals as well as the jumps across inner edges and Neumann boundary parts. These modifications mainly result from appropriate intersections of Neumann boundary parts with the enclosing domain. Important ingredients of the proof are the use of Stein's extension theorem and the application of a trace theorem which estimates the norm of the trace on (curved) boundary parts in terms of the local mesh size and polynomial degree. For the a posteriori error estimate it is necessary to assume that each mesh element near to Neumann boundary parts intersecting the interior of the enclosing domain is completely surrounded by one layer of mesh elements. This condition can easily be satisfied in the context of the FCM. In comparison to residualbased a posteriori error estimators for standard finite elements, the unknown multiplicative constant in the reliability estimate also depends on the constant used in Stein's extension theorem. In addition to the considerations on the reliability, we comment on the efficiency of the error estimator by providing a rather artificial example where the efficiency index (i.e., the ratio of the error estimator and the error) is not independent of the mesh-family parameter h. Therefore, we deem it unlikely that the error estimator can be proven to be efficient in the present form. However, our numerical examples based on more realistic domains suggest global efficiency with the occurrence of a large overestimation only on few cut elements (i.e., mesh elements with a non-empty intersection with the physical and fictitious domain). In the case of increasing polynomial degrees we observe the typical mild dependency of the efficiency index on the polynomial degree as described for hp-finite elements in [27].
The outline of this paper is as follows. In Sect. 2, the setting of the Finite Cell method with some assumptions on the configuration of the hp-FEM meshes and spaces is introduced. Section 3 recalls well-known results on Stein's extension operator and nonsmooth hp-interpolation operators. We summarize some statements on patches of mesh elements and prove the above-mentioned trace theorem. The main part of this section consists in the reliability proof of the error estimator. Moreover, we discuss the artificial example which hints at a theoretical inefficiency of the error estimator. In Sect. 4, numerical examples with nonsmooth solutions on two domains with polygonal holes and a domain with circular holes are provided. The 2D examples demonstrate the reliability for h-and p-uniform as well as for hp-geometric and h-adaptive refinements and confirm the efficiency for these more realistic examples. Although our theoretical results are presented for two-dimensional domains only, a final example demonstrates the practical applicability of the FCM error estimator to 3D problems as well.

The Finite Cell Method
In this section we introduce the Finite Cell method (FCM), where we use a similar setting as presented in [9]. We assume that the boundary of the domain Ω ⊂ R 2 is split into a closed, non-empty boundary part Γ D for Dirichlet boundary conditions and a boundary part Γ N for Neumann boundary conditions. We define V := {v ∈ H 1 (Ω); v = 0 on Γ D } with the usual Sobolev space H 1 (Ω) and assume f ∈ L 2 (Ω) and g ∈ L 2 (Γ N ). Poisson's problem consists in finding a solution u ∈ V such that for all v ∈ V . Here, (·, ·) ω denotes the L 2 inner product on ω ⊂ R 2 with the induced norm · 2 ω := (·, ·) ω . For the H 1 -norm we write · 2 1,ω := · 2 ω + ∇· 2 ω . To specify a discrete solution by means of the FCM, we assume an enclosing domainΩ with Ω ⊂Ω. In the practical application of the FCM, the enclosing domain is defined to be of simple shape, e.g. the union of (few) rectangles or triangles, so that a (finite element) mesh ofΩ can easily be constructed. By {T h } h∈H , H ⊂ R + , we denote a family of meshes T h ofΩ with closed triangles or parallelograms so that K ∩ K is empty or a vertex or edge of K and K for all K , K ∈ T h and K ∈T h K coincides with the closure ofΩ. In order to avoid dealing with the weak imposition of Dirichlet boundary conditions, we assume that Γ D is a subset of ∂Ω which is compatible with T h , i.e. Γ D is the union of closed edges of T h for all h ∈ H , so that the Dirichlet data can be imposed in a strong manner. ByΓ N we denote the part of Γ N whose closure is also assumed to be compatible with {T h } h∈H . Since the main application of the FCM consists in the treatment of domains with complicated boundaries, we allow the boundary partΓ N := Γ N \Γ N to have a curved shape (in contrast to the piecewise linear boundary parts Γ D andΓ N ), see Fig. 1. For this purpose, we assume the existence of We assume that the family of meshes {T h } h∈H is κ-shape regular with κ > 0, i.e.
for all K ∈ T h and h ∈ H . Here, F K :K → K is an affine-linear element mapping with F K (K ) = K and h K denotes the diameter of K . The reference elementK is the unit simplex if K is a triangle and it is the unit square if K is a parallelogram. The condition (2) implies that the sizes of neighboring elements are comparable: there exists γ ≥ 1 such that γ −1 h K ≤ h K ≤ γ h K for neighboring elements K , K ∈ T h with K ∩ K = ∅ for all h ∈ H . Moreover, from (2) we easily conclude the quasi-uniformity of {T h } h∈H for a constantγ > 0, i.e. h K ≤γ ρ K for all K ∈ T h and h ∈ H , where ρ K denotes the diameter of the in-circle of K ∈ T h if K is a triangle. If K is a parallelogram, ρ K is the smallest diameter of an in-circle contained in a triangle resulting from the bisection of K along one of the diagonals. Let p h : T h → N be a polynomial degree distribution such that the polynomial degrees of neighboring elements are comparable in the sense that there exists a constantγ ≥ 1 such thatγ −1 p h K ≤ p h K ≤γ p h K . Herewith, we define the hp-finite element space .. p ifK is a rectangle and P p (K ) = span x i y j i+ j= p ifK is a triangle. Given a parameter 0 < ε 1, the application of the FCM consists in finding the discrete solution u h ∈ V h such that for all v h ∈ V h . It is obvious that (3) is a perturbed version of the usual discrete weak formulation of Poisson's problem, where ε serves as a stabilization parameter ensuring that the left hand side of (3) is associated with a positive definite bilinear form. The discrete solution of the FCM uniquely exists for all ε > 0, however in computational practice, the parameter is often chosen in the range of 10 −14 up to 10 −8 . Clearly, ε introduces a modeling error and it should be chosen in such a way that the modeling error and the discretization error are properly balanced.
In the standard application of the FCM, the integrals in (3) are usually approximated by using an appropriate (non-standard) quadrature scheme as the domain Ω may not coincide with any union of elements of T h . Such a quadrature scheme may be based on, for instance, smart octrees [24] or a moment-fitting method [21] as, e.g., described in Sect. 4.3. While some of these quadrature schemes provide exact integration up to machine precision, others introduce a quadrature error which should also be in balance with the discretization error. In the following, we suppose that the integrals are computed exactly. Otherwise, additional error terms resulting from the quadrature need to be taken into account [11].
When the FCM is employed, one is typically interested in the restriction u h | Ω rather than u h . Moreover, only the error u − u h | Ω and norms thereof are relevant. The following Lemma states that the H 1 semi-norm of the discrete solution in the complementΩ \ Ω can become large as it may behave like 1/ε 1/2 . Consequently, the occurrence of this norm should be avoided in error estimates for the FCM.

Lemma 1 It holds
Proof Poincaré's inequality and the trace theorem imply Taking the square root and multiplying by ε 1/2 give the assertion.

Preliminaries
In this subsection, we summarize some basic properties of element patches, recall the quasiinterpolation operator for the hp setting as introduced in [26] and present some modifications of two other ingredients of error control, which are tailored to the FCM setting, namely, a modified Galerkin orthogonality property and a trace theorem. Some properties of Stein's extension operator are also recalled. We start by defining element patches for vertices v ∈ V h , edges e ∈ E h , and elements K ∈ T h by From the shape regularity of {T h } h∈H we conclude that each K ∈ T h has at most n(γ ) neighbors. Hence, with ν(γ ) := 1 + n(γ ) it follows by induction that Trivially, we have for some edges e 1 and e 2 of T h . This means that there exists a vertex Q of e 1 or e 2 such that d h

From [27, Theorem 2.2] we recall the quasi-interpolation operator
From this interpolation estimate we easily conclude the following estimates for a K ∈ T h and an edge e and state the stability property We note that the Galerkin orthogonality does not hold in the context of the FCM. However, subtracting (3) from (1) and using the fact One important ingredient in the subsequent section is the use of Stein's extension operator E : V →V for the extension of H 1 functions on Ω to functions onΩ. The operator fulfills (Ev)| Ω = v and for v ∈V [34]. The crucial aspect is that the increase in the norm of the extension depends only on the geometry of Ω: The second inequality in (10) is clear. The first inequality directly results from the extension property (9) and Poincaré's inequality, Ev 1,Ω v 1,Ω ∇v Ω . Finally, we state the following trace theorem for the images T (Γ ) and T (R) of the sets where T :

Lemma 3 It holds
From the fundamental theorem of calculus and the mean-value theorem we conclude the existence of some Thus, Hölder's inequality gives A density argument gives the assertion.
Applying this trace theorem to some specific set R b with b defined as the smallest in-circle diameter of a certain patch-layer divided by the polynomial degree, we obtain the following result.
Proof From the assumption d h 2 (K ) > 0 and Lemma 2 we conclude that

Reliability
We define the residual-based error estimator η as Here, h e is the length of e ∈ E h , p h e := min{ p h K | e ⊂ ∂ K , K ∈ T h } and n is a fixed normal unit vector to e which coincides with the outer normal on Γ N . Furthermore, the set E h (K ) contains the edges of K ∈ T h and [·] e denotes the jump across e.
In this subsection we prove the reliability of η, where we follow the reliability proof of the residual-based error estimator for Poisson's problem in the standard FEM setting [27], but apply several modifications in order to take the non-conformity of the boundary into account.
Proof The main idea of the proof is to use the Stein extension δ * := Eδ of the error δ := u − u h | Ω in order to apply the interpolation operator I h to obtain estimates on whole elements K ∈ T h which then can be estimated by norms on K ∩ Ω. We set w * := δ * − I h δ * and note that w * | Ω ∈ V . Thus, by using the fact that u is a solution of (1) we obtain Note (∇δ, ∇ I h δ * ) Ω = 0 as the Galerkin orthogonality does not hold in the FCM, see (8).
Applying integration by parts on Ω we obtain For the boundary integrals we observe that where all unions are pairwise disjoint except for sets of boundary measure zero. Exploiting w * | Γ D = 0 and using jumps across inner edges e ⊆ Ω we get We split (g, w * ) Γ N into edge contributions onΓ N and onΓ N and combine them with the terms in (15), which gives Inserting (16) in (14) and using Cauchy's inequality, we obtain To estimate the first, second, and third term of (17), we exploit K ∩ Ω ⊆ K , e ∩ Ω ⊆ e and apply (6) to obtain Here, we omit the index h at h K , h e , p K and p e to ease the notation. From (4) we conclude that the number of patches of level m ∈ N containing a specific K ∈ T h is bounded independently of h, more precisely and Hence, we have Applying (18) and Cauchy's inequality give Analogously, we derive II The bound of w * Γ N ∩K in the fourth term of (17) is the most delicate. It holds for each for some m K ∈ N and graphs Γ K i as defined in (11) and some counterclockwise rotations with translation T K i : R 2 → R 2 . From Lemma 4 we conclude where we apply the stability property (7) to estimate ∇w * ℘ h 2 (K ) and exploit the γ -shape regularity and the comparability of the polynomial degrees yielding Hence, we obtain from (20) Since m K is bounded by the total number of segments ofΓ N , we have and, therefore, where we use Cauchy's inequality and which results from (19). To estimate the fifth term in (17) we again use (19) to obtain The stability property (7) then yields so that exploiting the modified Galerkin orthogonality (8) and applying Lemma 1 give Finally, inserting the estimates for I, II, III, IV, and V in (17) and using the extension property (9) we obtain which gives the assertion.

Remark 1
The assumption on the diameter of patch-layers in Theorem 1 means that all K ∈ T h with K ∩Γ K = ∅ have to be completely enclosed by the patch-layer L h 2 (K ). This can be ensured by a sufficiently fine mesh in combination with a sufficiently large enclosing domainΩ in the neighborhood of K ∩Γ N = ∅.

Some Remarks on Efficiency and Overestimation
The derivation of the error estimates in Sect. 3.2 are based on the use of a Clément-or Scott-Zhang-type interpolation operator which typically ignores the possibly complicated geometry of cut elements K ∩ Ω. Hence, the estimates contain factors h K instead of factors which are adapted to the size of the cut elements. To show the efficiency of a residual-based error estimator, inverse estimates such as ∇u h K h −1 K p 2 K u h K are typically applied. However, they cannot be utilized on K ∩Ω. Instead, if Ω is for instance polygonal, one may use a subtriangulation T K of K ∩ Ω with K ∩ Ω = T ∈T K T . Then, for any T ∈ T K we would have where C K depends on the shape of K ∩ Ω. For standard finite element approximations, the proof of efficiency exploits the fact that the factor h K gained by the interpolation estimate cancels out the factor h −1 K incurred by the inverse estimate. In the FCM situation the factor h −1 T resulting from the subtriangulation could be significantly larger than the factor h K . Clearly, this may not be a problem if a subtriangulation with comparable element diameters fulfilling c −1 h T ≤ h K ≤ ch T is available for some constant c ≥ 1 independent of h ∈ H . However, elements could be very thinly cut so that each subtriangulation will admit triangles of a very small diameter h T h K and it is not clear whether a subtriangulation with comparable element diameters exists for each possible cut of elements with Ω. In summary, the mismatch between the factors introduced by the interpolation operator and the factor incurred by the inverse estimate may lead to inefficient error estimates or to extremely large overestimations, respectively. The situation may be remedied if more localized interpolation estimates are available.
In order to illustrate that the estimator (13) may lead to large overestimation or may even be inefficient, we study Poisson It is easy to see that {T h } 0<h<1/2 fulfills the shape regularity condition (2). As a finite element space on where ϕ h i , i = 1, 2, 3, 4, are the piecewise bilinear functions associated to the vertices outside of Γ D . Since the basis functions are linearly independent on Ω, we set ε := 0 in the FCM discretization given by (3). In particular, no conditioning issues occur if we solve the resulting linear system exactly. Denoting the discrete solution by u h , we compute the local error on K 2 , where we have to take only the part of K 2 into account which is associated with Ω, i.e. K 2 ∩ Ω = [−1 − h, −1) × (0, 1), see "Appendix A" for the details of this computation. We eventually obtain The error estimator terms related to K h 2 are h 2 where we exploit Δu h = 0 and p = 1. Since h K h we obtain a factor of 1/h 2 for the inner term (23) and 1/h for the Neumann boundary term (24) when we relate these terms with the exact error (22). This means that the error estimator (13) does not yield a lower bound with a constant independent of h in this example and is therefore not efficient. Moreover, the error is largely overestimated by the error estimator for small h. We emphasize that we consider this artificial counter example to demonstrate the possible overestimation and inefficiency of the error estimator. Indeed, the numerical examples of Sect. 4 show that the estimator seems to be efficient with a moderate overestimation in practical applications (at least in the h-version of the FCM). In the p-and hp-version of the FCM we observe the typical p-dependency of the efficiency index.

Numerical Examples
In this section, we illustrate the reliability of the residual-based error estimator of Sect. 3.2. For this purpose we study Poisson's problem with non-smooth solutions in 2D and 3D. We adapt two classical problems so that an application of the Finite Cell method is convenient.
First, we consider a modified version of the L-shaped domain problem where the L-shape domain is equipped with some additional holes. We study the behavior of the error estimator applied to various configurations including h-uniform, p-uniform, hp-geometric, and hadaptive refinements. As the exact solution is known for the L-shaped domain problem, we are able to compute the exact error (up to machine precision) and check the overestimation of the error by the error estimator η from (13). In particular, we study the efficiency index eff := η/ ∇u − ∇u h Ω . To study the local overestimation and the efficiency of the error estimator, we introduce a (maximum) local efficiency index loc := max K ∈T h η K / ∇u − ∇u h K ∩Ω . In light of the considerations of Sect. 3.3, we check whether a large local overestimation occurs and whether it disturbs global efficiency. Note that a large local overestimation may mainly occur on cut cells or elements in the neighborhood of a cut cell. If we apply only p-refinements with increasing polynomial degree, we expect the estimator to be efficient up to a factor O( p), see [27]. In the second example, we consider the L-shaped domain problem from the first example but replace the polygonal holes by circular ones. In particular, we deal with curved boundaries. In the third example, we apply the estimator to drive h-adaptive refinements for the 3D example of Fichera's corner with holes. In all examples, we aim to examine whether optimal algebraic convergence rates can be recovered by the application of h-refinements driven by the error estimator.

L-Shape with Polygonal Holes
The first example is given by Poisson's problem on a domain Ω with several holes. It is obtained by removing 11 polygonal holes ω 1 , . . . , ω 11 of various sizes and shapes from the standard L-shaped domain, which serves as the enclosing domainΩ, see Fig. 4. In accordance with the compatibility requirements on Γ D and Γ N as introduced in Sect. 2, the holes neither touch the outer boundary nor do they touch each other. The weak formulation of the problem reads: Find u ∈ V such that The data g is derived from the nonsmooth solution u(r , ϕ) := r 2/3 (sin(2ϕ − π)/3) given in polar coordinates, which yields a corner singularity in the origin. The FCM discretization consists in seeking u h ∈ V h such that where ε := 10 −12 and V h is a FEM space based on a mesh T h ofΩ. Let the set of cut cells be defined asT h := {K ; K ∈ T h , K ∩ Ω = K }. To obtain a suitable quadrature on a cut K ∩ Ω with K ∈T h , we perform a Delaunay triangulation of K ∩ Ω using the well-known computational geometry library CGAL [39] and apply an appropriate quadrature rule on each triangle of this triangulation.
We test the error estimator in four configurations. In all cases, we start with an initial mesh consisting of three quad elements with edge length 1. In the first configuration, we perform uniform h-refinements of the initial mesh, i.e. each element is repeatedly bisected into four congruent quads. The influence of the polynomial degree on the estimation is studied Fig. 4 Domain Ω resulting from the removal of 11 polygonal holes from the L-shaped domain Ω with initial mesh consisting of three quad elements for p = 1 and p = 2. Figure 5a shows the decay of the exact error and the reliable error estimator (13). It can be seen that the convergence order of O(DOF −1/3 ) = O(h 2/3 ) is attained for p = 1 as well as p = 2 which is expected due to u ∈ H 5/3−ε (Ω), ε > 0 [22]. The global efficiency index is depicted in Fig. 5b and suggests that the estimator captures the error with acceptable indices between 3 and 6. In particular, we observe that the indices lie in a more constricted range for smaller h. In contrast, the local efficiency indices as shown   Fig. 5c exhibit a more fluctuating behavior. Indeed, the large factors of approx. 80 hint at the large local overestimation. However, taking into account the moderate global efficiency index, the overestimation only occurs on elements that bear a small portion of the global error, so that the overestimation has little impact on the overall estimate. In Fig. 6a, a colormap indicates the local overestimation. Obviously, the overestimation is moderate on elements which are not cut by the holes. The largest overestimation of 76.8 occurs on an element K in the lower-left. Zooming into K in Fig. 6b and c, we see that only a very small portion of the element is contained in Ω. Indeed, the exact error on K ∩ Ω is approx. 2.7 × 10 −7 and the estimation is 2.1 × 10 −5 . However, the average estimator value per element is computed to be approx. 1.6 × 10 −4 meaning that the estimation on K ∩ Ω is well below average and hence has no detrimental effect on the overall estimation.
In the second configuration, we keep the mesh fixed at three elements and perform prefinements only. Here, we make use of a spline-based finite element space V h ⊂ C 1 described in [23]. The differentiable discrete functions of V h have the advantage that the jump terms [∂ n u h ] e∩Ω across inner edges disappear. In the context of the FCM, the use of these functions greatly simplifies the evaluation of the estimator since it avoids the generation of inner edge cuts e ∩ Ω which are actually not needed during the solution of the discrete problem. Figure 7a shows the errors and estimates resulting from p = 3, . . . , 20 indicating that the convergence order of O(DOF −2/3 ) = O( p −4/3 ) is attained [36]. To assess the efficiency of the estimation, we study the global and local efficiency indices which are shown in Fig. 7b. We expect global and local efficiency that may depend on a factor O( p), see [27]. Indeed, both eff and loc exhibit a p-dependent growth. Therefore, we also plot eff / p, loc / p. Since the scaled indices appear to be constant, we conclude that the multiplicative factor of O( p) does occur in this example.
In the third configuration, we perform hp-geometric refinements in which the degree as well as the mesh size are decreased towards the corner singularity, see Fig. 8a. To this end, we utilize the C 1 -conforming finite element space from [23] as well, which allows for multilevel hanging nodes and arbitrary degree distributions. It is well-known that the geometric refinements lead to an exponential convergence of the error of the form O(exp(−bDOF 1/3 )), b > 0 [36]. Indeed, the straight lines in the semilogarithmic plot in Fig. 8b confirm this type of convergence. The global and local efficiency indices are visualized in Fig. 8c. We observe that the indices increase in the beginning, but after a certain number of refinements, the indices seem to become constant. The reason for this may be the fact that only the elements touching the reentrant corner are refined. Hence, the geometry of the cuts K ∩ Ω only varies during the first few refinements where the new elements intersect the holes. In the fourth configuration, we test h-adaptive refinements driven by the estimator based on C 0 elements of degree p = 1 and p = 2, respectively. Figure 10a shows the decay of the exact error and the reliable estimation (13). We observe that optimal algebraic convergence rate of O(DOF − p/2 ) is recovered. The convergence rates computed using two consecutive errors are listed in the Tables 1 and 2. The resulting meshes are visualized in Fig. 9. We see a high resolution of the corner singularity by the refinements. The meshes exhibit only very few unnecessary refinements near the boundaries of the holes. The global efficiency index is depicted in Fig. 10b. As the index ranges between 2 and 6, the overestimation seems to be acceptable. The large local efficiency indices of up to 60 in Fig. 10c hint at a large local overestimation. However, the overestimation seems to occur on elements with a small portion of the overall error since the global efficiency index is small and the optimal convergence rate is attained.

L-Shape with Circular Holes
In this example, we assess the error estimator in the presence of curved boundaries and study the influence of the parameter ε on the estimator. To this end, the problem setup from  11 . The domain Ω, the holes, and an initial mesh are displayed in Fig. 11. Recall that, in the FCM, the problem of constructing a mesh for Ω is shifted to the numerical integration. In particular, the domain-approximation error is negligibly small if integrals on Ω are computed with machine precision. For the computation of the integrals on a cut K ∩ Ω for K ∈T h up to machine precision, we employ the moment-fitting approach from [20], in which a customized quadrature rule for each cut is determined. Since the computation of the integrals on cuts is one of the crucial issues in the FCM, we will discuss the moment-fitting approach in more detail. The quadrature rule obtained by moment fitting , the function Φ i can easily be obtained by computing the antiderivative of one of the functions ϕ i,k . The computation of the surface integral in (25) is much simpler and computationally cheaper than the computation of the volume integral since only parameterizations of surfaces are involved. In this example, the integrals on parts of circles have the form t 2 t 1 ϕ i (cos(t), sin(t)) dt and, thus, the integrand is a univariate polynomial in cos(t) and sin(t). To obtain a sufficiently accurate quadrature rule, we distribute p + 60 Gauss points on (t 1 , t 2 ), where p is the degree of the polynomial ϕ i . Note that the increased number of Gauss points is only used for the computation of the surface integral in (25) and, hence, the number of Gauss points used for the assembly is ( p + 1) 2 per element as expected.
In this numerical example, we use C 0 finite elements of degree p := 10 and perform mesh adaptivity using the error estimator. In order to additionally assess the influence of ε, we discuss four configurations with the parameter ε ∈ 10 −1 , 10 −2 , 10 −4 , 10 −6 . The decay of the error and the estimator is displayed in Fig. 12.
The upper bound ∇u − ∇u h Ω η + ε 1/2 from (13) suggests that the error may be as large as O(ε 1/2 ). Thus, the error may become stationary after reaching the magnitude of O(ε 1/2 ), i.e., no further reduction of the error can be achieved even if the DOF are increased. However, in this example, we notice that the error becomes stationary when ∇u − ∇u h Ω ≈ O(ε). This is considerably better than the worst-case bound O(ε 1/2 ) that results from the use of the bound ε ∇u h Ω \Ω ε 1/2 ( f Ω + g Γ N ) in (21). Here, f = 0 gives which explains the improved behavior in this example.
The experiment confirms that ε should be small in order to obtain a sufficiently accurate solution. In fact, the errors are much larger for ε ∈ 10 −1 , 10 −2 than for ε ∈ 10 −4 , 10 −6 , see Fig. 12. In the latter case, the errors are even quite similar (until the error for ε = 10 −4 becomes stationary). We see that the estimator has approximately the same value for each ε until its decay rate reduces. The dependencies of the error and the estimator on ε seem to differ. At a certain number of DOF the error becomes stationary whereas the decay rate of the estimator is only reduced. However, the numbers of DOF at which the behavior of the error and the estimator changes are quite similar for small ε. This may also be seen in Fig. 13 where the efficiency indices are depicted.

Fichera's Corner with Polyhedral Holes
The error estimator (13) is stated for the 2D case. The extension of the individual error contributions to 3D is straight forward. However, it is not clear whether the arguments for proving reliability can easily be transfered from the two-dimensional case. Nevertheless, we study the 3D extension of the error estimator (13) in a numerical experiment. For this purpose, we consider Fichera's corner which results from removing the unit cube [0, 1] 3 from the cube (−1, 1) 3   with the polar coordinates (r (w 1 , w 2 ), ϕ(w 1 , w 2 )) ∈ R + × [0, 2π) of (w 1 , w 2 ) ∈ R 2 . The exact solution of the problem is unknown. It features edge singularities at the edges emanating from the origin and a corner singularity at the origin. The discrete formulation in the FCM version reads: Find u h ∈ V h such that where V h is a finite element space based on tensor-product hexahedral elements of fixed degree p = 1 or p = 2. Furthermore, we set ε := 10 −12 . The initial coarse mesh consists of the seven cubes as introduced above.
In order to compute the integrals on a cut K ∩Ω for some K ∈T h up to machine precision, we employ the moment-fitting approach described in Sect. 4.2. Here, the computation of the surface integrals is much cheaper than the computation of a volume integral since only surface meshes are involved which essentially amounts to the intersection of 3D triangle meshes (e.g. implemented in CGAL [3]), while a volume mesher is needed to obtain a tetrahedralization of K ∩ Ω. Moreover, the number of tetrahedrons in a volume mesh is usually much larger than the number of triangles on the surface.
As the exact solution is unknown, we only study the convergence of the estimated error, which is shown in Fig. 15 for h-adaptive refinements with p = 1 and p = 2. It can be seen that the optimal algebraic convergence rate of O(DOF − p/2 ) is obtained in both cases. The adaptive mesh and a parallel cut through the mesh are shown in Fig. 16. Refinements occur in the neighborhood of the corner and edge singularities. Moreover, there are no refinements inside the holes and only few refinements appear on the boundaries of the holes. The optimal convergence rate indicates that a large overestimation of the error on cut cells which also cause a large global overestimation does not occur.
Funding Open access funding provided by Paris Lodron University of Salzburg.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.