Reliable Residual-Based Error Estimation for the Finite Cell Method

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 hp-adaptive 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 residual-based 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.

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 \(\varOmega \subset \mathbb {R}^2\) is split into a closed, non-empty boundary part \(\varGamma _D\) for Dirichlet boundary conditions and a boundary part \(\varGamma _N\) for Neumann boundary conditions. We define \(V :=\{v\in H^1(\varOmega ) ;\ v=0 \text { on } \varGamma _D\}\) with the usual Sobolev space \(H^1(\varOmega )\) and assume \(f \in L_2(\varOmega )\) and \(g \in L_2(\varGamma _N)\). Poisson’s problem consists in finding a solution \(u\in V\) such that

for all \(v\in V\). Here, \((\cdot ,\cdot )_\omega \) denotes the \(L_2\) inner product on \(\omega \subset \mathbb {R}^2\) with the induced norm \(\left\Vert \cdot \right\Vert _\omega ^2 :=(\cdot , \cdot )_\omega \). For the \(H^1\)-norm we write \(\left\Vert \cdot \right\Vert _{1,\omega }^2 :=\left\Vert \cdot \right\Vert _\omega ^2 + \left\Vert \nabla \cdot \right\Vert _\omega ^2\).

To specify a discrete solution by means of the FCM, we assume an enclosing domain \({\hat{\varOmega }}\) with \(\varOmega \subset {\hat{\varOmega }}\). 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 \({\hat{\varOmega }}\) can easily be constructed. By \(\{\mathscr {T}_h\}_{h\in {\mathscr {H}}}\), \({\mathscr {H}}\subset \mathbb {R}_+\), we denote a family of meshes \(\mathscr {T}_h\) of \({\hat{\varOmega }}\) with closed triangles or parallelograms so that \(K\cap K'\) is empty or a vertex or edge of K and \(K'\) for all \(K,K'\in \mathscr {T}_h\) and \(\bigcup _{K\in \mathscr {T}_h} K\) coincides with the closure of \({\hat{\varOmega }}\). In order to avoid dealing with the weak imposition of Dirichlet boundary conditions, we assume that \(\varGamma _D\) is a subset of \(\partial {\hat{\varOmega }}\) which is compatible with \(\mathscr {T}_h\), i.e. \(\varGamma _D\) is the union of closed edges of \(\mathscr {T}_h\) for all \(h\in {\mathscr {H}}\), so that the Dirichlet data can be imposed in a strong manner. By \({\bar{\varGamma }}_N\) we denote the part of \(\varGamma _N\) whose closure is also assumed to be compatible with \(\{\mathscr {T}_h\}_{h\in {\mathscr {H}}}\). Since the main application of the FCM consists in the treatment of domains with complicated boundaries, we allow the boundary part \({\hat{\varGamma }}_N :=\varGamma _N \setminus {\bar{\varGamma }}_N\) to have a curved shape (in contrast to the piecewise linear boundary parts \(\varGamma _D\) and \({\hat{\varGamma }}_N\)), see Fig. 1. For this purpose, we assume the existence of functions \(\varphi _i\in C^1([a_1^i,a_2^i])\) with \(a_1^i<a_2^i\), \(i=1,\ldots ,m\) and counterclockwise rotations with translation \(T_i: \mathbb {R}^2 \rightarrow \mathbb {R}^2\) such that the closure of \({\hat{\varGamma }}_N\) coincides with \(\bigcup _{i=1}^m T_i(\varGamma _i)\) where \(\varGamma _i\) is the graph of \(\varphi _i\), i.e. \(\varGamma _i :={\left\{ (x,\varphi (x)) ;\ a_1^i \le x \le a_2^i \right\} }\).

Configuration of \(\varOmega \), \({\hat{\varOmega }}\), and a Finite Cell mesh

Full size image

We assume that the family of meshes \(\{\mathscr {T}_h\}_{h\in {\mathscr {H}}}\) is \(\kappa \)-shape regular with \(\kappa >0\), i.e.

for all \(K\in \mathscr {T}_h\) and \(h\in {\mathscr {H}}\). Here, \(F_K:{\hat{K}}\rightarrow K\) is an affine-linear element mapping with \(F_K({\hat{K}})=K\) and \(h_K\) denotes the diameter of K. The reference element \({\hat{K}}\) 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 \(\gamma \ge 1\) such that \(\gamma ^{-1}h_K\le h_{K'}\le \gamma h_K\) for neighboring elements \(K, K' \in \mathscr {T}_h\) with \(K \cap K' \ne \varnothing \) for all \(h\in {\mathscr {H}}\). Moreover, from (2) we easily conclude the quasi-uniformity of \(\{\mathscr {T}_h\}_{h\in {\mathscr {H}}}\) for a constant \({\hat{\gamma }}>0\), i.e. \(h_K\le {\hat{\gamma }} \rho _K\) for all \(K\in \mathscr {T}_h\) and \(h\in {\mathscr {H}}\), where \(\rho _K\) denotes the diameter of the in-circle of \(K\in T_h\) if K is a triangle. If K is a parallelogram, \(\rho _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:\mathscr {T}_h\rightarrow \mathbb {N}\) be a polynomial degree distribution such that the polynomial degrees of neighboring elements are comparable in the sense that there exists a constant \({\tilde{\gamma }}\ge 1\) such that \({\tilde{\gamma }}^{-1} p_K^h\le p_{K'}^h\le {\tilde{\gamma }} p_K^h\). Herewith, we define the hp-finite element space

with \({\hat{V}} :={\left\{ v\in H^1({\hat{\varOmega }}) ;\ v=0 \text { on } \varGamma _D \right\} }\), where \(P_p({\hat{K}})={\text {span}} {\left\{ x^iy^j \right\} }_{i,j=0,\ldots p}\) if \({\hat{K}}\) is a rectangle and \(P_p({\hat{K}})={\text {span}} {\left\{ x^iy^j \right\} }_{i+j=p}\) if \({\hat{K}}\) is a triangle.

Given a parameter \(0 < \varepsilon \ll 1\), the application of the FCM consists in finding the discrete solution \(u_h \in V_h\) such that

for all \(v_h \in V_h\). It is obvious that (3) is a perturbed version of the usual discrete weak formulation of Poisson’s problem, where \(\varepsilon \) 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 \(\varepsilon > 0\), however in computational practice, the parameter is often chosen in the range of \(10^{-14}\) up to \(10^{-8}\). Clearly, \(\varepsilon \) 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 \(\varOmega \) may not coincide with any union of elements of \(\mathscr {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|_{\varOmega }\) rather than \(u_h\). Moreover, only the error \(u - u_h|_{\varOmega }\) and norms thereof are relevant. The following Lemma states that the \(H^1\) semi-norm of the discrete solution in the complement \({\hat{\varOmega }} \setminus \varOmega \) can become large as it may behave like \(1/\varepsilon ^{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

which yields \(\left\Vert \nabla u_h\right\Vert _{\varOmega } \le \left\Vert f\right\Vert _{\varOmega } + \left\Vert g\right\Vert _{\varGamma _N}\) and, thus, \(\varepsilon \left\Vert \nabla u_h\right\Vert _{\widehat{\varOmega } \setminus \varOmega }^2 \le {(\left\Vert f\right\Vert _{\varOmega } + \left\Vert g\right\Vert _{\varGamma _N})}^2\). Taking the square root and multiplying by \(\varepsilon ^{1/2}\) give the assertion. \(\square \)

3.1 Preliminaries

In this subsection, we summarize some basic properties of element patches, recall the quasi-interpolation 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 \in {\mathscr {V}}_h\), edges \(e \in {\mathscr {E}}_h\), and elements \(K \in \mathscr {T}_h\) by

and recursively define

for a node \(N \in {\mathscr {V}}_h\cup {\mathscr {E}}_h \cup \mathscr {T}_h\) and \(j \ge 2\). We set \({\overline{\wp }}_j^h(N) :=\bigcup _{K\in \wp _j^h(N)} K\), \(j\ge 1\). From the shape regularity of \(\{\mathscr {T}_h\}_{h\in \mathscr {H}}\) we conclude that each \(K\in \mathscr {T}_h\) has at most \(n(\gamma )\) neighbors. Hence, with \(\nu (\gamma ) :=1 + n(\gamma )\) it follows by induction that

Trivially, we have

Let \(\rho _m^h(K) :=\min _{K'\in {\mathscr {L}}_m^h(K)} \rho _{K'}\) be the smallest in-circle diameter of the patch-layer \({\mathscr {L}}_m^h(K) :=\wp _m^h(K)\backslash \wp _{m-1}^h(K)\) for \(K\in \mathscr {T}_h\). In the next Lemma, we observe that \(\rho _m^h(K)\) is smaller than the diameter

of \({\mathscr {L}}_m^h(K)\), see Fig. 2a.

Lemma 2

Let \(K\in \mathscr {T}_h\) and \(d_m^h(K)>0\), then \(\rho _m^h(K)\le d_m^h(K)\).

Proof

Since \(\partial {\overline{\wp }}_m^h(K)\) and \(\partial {\overline{\wp }}_{m-1}^h(K)\) consist of edges of \(\mathscr {T}_h\), we have

for some edges \(e_1\) and \(e_2\) of \(\mathscr {T}_h\). This means that there exists a vertex Q of \(e_1\) or \(e_2\) such that \(d_m^h(K)=\min _{y\in e} \Vert Q-y\Vert \) with \(e\in \{e_1,e_2\}\). Since \(d_m^h(K)>0\) we conclude that Q is a vertex and e is an edge of some \({\tilde{K}}\in {\mathscr {L}}_m^h(K)\), where Q is opposite to e. Thus, \(d_m^h(K)\ge \rho _{{\tilde{K}}}\ge \rho _m^h(K)\). \(\square \)

Geometric constructions for a triangle K (marked in bold) in a mesh \(\mathscr {T}_h\) of 10 elements

Full size image

From [27, Theorem 2.2] we recall the quasi-interpolation operator \(I_h: {\hat{V}} \rightarrow V_h\), which has the properties \({(I_h z)}|_{\varGamma _D} = z|_{\varGamma _D} = 0\) for \(z\in {\hat{V}}\) and

for all \(v \in {\mathscr {V}}_h\) and all edges \(e \in {\mathscr {E}}_h\) incident to v, where \(h_v^h :=\min _{K\in \wp ^h_1(v)} h_K\) and \(p_v^h :=\min _{K\in \wp _1^h(v)} p_K^h\). From this interpolation estimate we easily conclude the following estimates for a \(K\in \mathscr {T}_h\) and an edge \(e\in {\mathscr {E}}_h\)

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 that \(v_h|_{\varOmega } \in V\) for \(v_h\in V_h\) yield the modified Galerkin orthogonality

for all \(v_h \in V_h\).

One important ingredient in the subsequent section is the use of Stein’s extension operator \(\mathrm {E}: V \rightarrow {\hat{V}}\) for the extension of \(H^1\) functions on \(\varOmega \) to functions on \({\hat{\varOmega }}\). The operator fulfills \({(\mathrm {E}v)}|_{\varOmega } = v\) and

for \(v \in {\hat{V}}\) [34]. The crucial aspect is that the increase in the norm of the extension depends only on the geometry of \(\varOmega \): For \(v\in V\) we have

The second inequality in (10) is clear. The first inequality directly results from the extension property (9) and Poincaré’s inequality, \( \left\Vert \mathrm {E}v\right\Vert _{1,{\hat{\varOmega }}} \lesssim \left\Vert v\right\Vert _{1,\varOmega } \lesssim \left\Vert \nabla v\right\Vert _{\varOmega } \).

Finally, we state the following trace theorem for the images \(T(\varGamma )\) and T(R) of the sets

where \(T:\mathbb {R}^2\rightarrow \mathbb {R}^2\) is a counterclockwise rotation with translation, \(b>0\), and \(\varphi \in C^1([a_1,a_2])\) with \(a_1<a_2\).

Lemma 3

It holds

for all \(v \in H^1(T(R_b))\).

Proof

Let \(v \in C^1(T(R_b))\) and \(a_1 \le x \le a_2\). From the fundamental theorem of calculus and the mean-value theorem we conclude the existence of some \(\eta \in (\varphi (x),\varphi (x)+b)\) such that

Thus, Hölder’s inequality gives

where \(0\le \alpha <2\pi \) is the rotation angle of T. With \(M :=\max _{x\in [a_1,a_2]} \left\Vert \varphi '(x)\right\Vert \) we have

A density argument gives the assertion. \(\square \)

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.

Lemma 4

Let \(K\in \mathscr {T}_h\), \(d_2^h(K)>0\) and \(T(\varGamma )\subset K\). It holds

for all \(v \in H^1({\overline{\wp }}_2^h(K))\).

Proof

From the assumption \(d_2^h(K) > 0\) and Lemma 2 we conclude that \(R_b\subset {\overline{\wp }}_2^h(K)\) with \(b :=\frac{\rho _2^h(K)}{p_K}\), see Fig. 2b. From Lemma 3 and

with \({\tilde{K}}\in {\mathscr {L}}_2^h(K)\) and \(\rho _{{\tilde{K}}}=\rho _2^h(K)\), we have

\(\square \)

3.2 Reliability

We define the residual-based error estimator \(\eta \) as

with

Here, \(h_e\) is the length of \(e\in {\mathscr {E}}_h\), \(p_e^h :=\min \{p_K^h\mid e\subset \partial K,\; K\in \mathscr {T}_h\}\) and n is a fixed normal unit vector to e which coincides with the outer normal on \(\varGamma _N\). Furthermore, the set \({\mathscr {E}}_h(K)\) contains the edges of \(K\in \mathscr {T}_h\) and \(\left[ \cdot \right] _e\) denotes the jump across e.

In this subsection we prove the reliability of \(\eta \), 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.

Theorem 1

Assume that \(d_2^h(K)>0\) for all \(K\in \mathscr {T}_h\) with \(K\cap {\hat{\varGamma }}_N\ne \emptyset \). Then,

Proof

The main idea of the proof is to use the Stein extension \(\delta ^* :=\mathrm {E}\delta \) of the error \(\delta :=u - u_h|_{\varOmega }\) in order to apply the interpolation operator \(I_h\) to obtain estimates on whole elements \(K \in \mathscr {T}_h\) which then can be estimated by norms on \(K \cap \varOmega \). We set \(w^* :=\delta ^* - I_h \delta ^*\) and note that \(w^*|_{\varOmega } \in V\). Thus, by using the fact that u is a solution of (1) we obtain

Note \({\left( \nabla \delta , \nabla I_h \delta ^*\right) }_{\varOmega } \ne 0\) as the Galerkin orthogonality does not hold in the FCM, see (8). Applying integration by parts on \(\varOmega \) we obtain

For the boundary integrals we observe that

where all unions are pairwise disjoint except for sets of boundary measure zero. Exploiting \(w^*|_{\varGamma _D}=0\) and using jumps across inner edges \(e \subseteq \varOmega \) we get

We split \((g,w^*)_{\varGamma _N}\) into edge contributions on \({\bar{\varGamma }}_N\) and on \({\hat{\varGamma }}_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 \cap \varOmega \subseteq K\), \(e \cap \varOmega \subseteq 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\in \mathbb {N}\) containing a specific \(K\in \mathscr {T}_h\) is bounded independently of h, more precisely

and

Hence, we have

Applying (18) and Cauchy’s inequality give

Analogously, we derive

The bound of \(\left\Vert w^*\right\Vert _{{\hat{\varGamma }}_N \cap K}\) in the fourth term of (17) is the most delicate. It holds for each \(K\in T_h\) that \(K\cap {\hat{\varGamma }}_N= \bigcup _{i = 1}^{m_K} T_i^K(\varGamma _i^K)\) for some \(m_K \in \mathbb {N}\) and graphs \(\varGamma _i^K\) as defined in (11) and some counterclockwise rotations with translation \(T_i^K:\mathbb {R}^2\rightarrow \mathbb {R}^2\). From Lemma 4 we conclude

where we apply the stability property (7) to estimate \(\left\Vert \nabla w^*\right\Vert _{{\overline{\wp }}^h_2(K)}\) and exploit the \(\gamma \)-shape regularity and the comparability of the polynomial degrees yielding

for all \(K' \in \wp ^h_2{(K)}\). Using Lemma 5 gives

Hence, we obtain from (20)

Since \(m_K\) is bounded by the total number of segments of \({\hat{\varGamma }}_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. \(\square \)

Remark 1

The assumption on the diameter of patch-layers in Theorem 1 means that all \(K\in \mathscr {T}_h\) with \(K\cap {\hat{\varGamma }}_K\ne \emptyset \) have to be completely enclosed by the patch-layer \({\mathscr {L}}_{2}^h(K)\). This can be ensured by a sufficiently fine mesh in combination with a sufficiently large enclosing domain \({\hat{\varOmega }}\) in the neighborhood of \(K\cap {\hat{\varGamma }}_N\ne \emptyset \).

3.3 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 \cap \varOmega \). 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

are typically applied. However, they cannot be utilized on \(K \cap \varOmega \). Instead, if \(\varOmega \) is for instance polygonal, one may use a subtriangulation \(\mathscr {T}_K\) of \(K \cap \varOmega \) with \(K \cap \varOmega = \bigcup _{T \in \mathscr {T}_K} T\). Then, for any \(T \in \mathscr {T}_K\) we would have

where \(C_K\) depends on the shape of \(K \cap \varOmega \). 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_K^{-1}\) incurred by the inverse estimate. In the FCM situation the factor \(h_T^{-1}\) 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 \le h_K \le c h_T\) is available for some constant \(c\ge 1\) independent of \(h\in {\mathscr {H}}\). However, elements could be very thinly cut so that each subtriangulation will admit triangles of a very small diameter \(h_T \ll h_K\) and it is not clear whether a subtriangulation with comparable element diameters exists for each possible cut of elements with \(\varOmega \). 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’s problem on the rectangular domain \(\varOmega :=(-2, -1) \times (0,1)\) with the solution \(u(x,y) :=(x+2)^2 \in V\). The Dirichlet boundary is set to \(\varGamma _D :={\left\{ -2 \right\} } \times [0,1]\), see Fig. 3. The Neumann boundary is given by \(\varGamma _N :=\varGamma _1 \cup \varGamma _2 \cup \varGamma _3\) with \(\varGamma _1 :=(0,1] \times {\left\{ 0 \right\} }\), \(\varGamma _2 :={\left\{ 1 \right\} } \times [0,1]\) and \(\varGamma _3 :=(0,1] \times {\left\{ 1 \right\} }\).

Setup of an artificial example. The domain \(\varOmega = (-2,-1) \times (0,1)\) is shaded

Full size image

With \(f :=-\varDelta u = -2\) and \(\nabla u(x,y) = (2(x+2),0)\) the weak formulation reads

for all \(v \in V\). The enclosing domain \({\hat{\varOmega }} :=(-2,0) \times (0,1)\) is bisected into a family of meshes \(\{\mathscr {T}_h\}_{0<h<1/2}\) with \(\mathscr {T}_h :=\{K_{1}^h, K_{2}^h\}\), where \(K_{1}^h :=[-2, -1-h] \times [0,1]\) and \(K_2^h :=[-1-h, 0] \times [0,1]\). It is easy to see that \(\{\mathscr {T}_h\}_{0<h<1/2}\) fulfills the shape regularity condition (2). As a finite element space on \(\mathscr {T}_h\), we choose \(V_h :={\text {span}}{\left\{ \varphi _{1}^h, \varphi _{2}^h, \varphi _{3}^h, \varphi _{4}^h \right\} }\) where \(\varphi _i^h\), \(i=1,2,3,4\), are the piecewise bilinear functions associated to the vertices outside of \(\varGamma _D\). Since the basis functions are linearly independent on \(\varOmega \), we set \(\varepsilon :=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 \(\varOmega \), i.e. \(K_2 \cap \varOmega = [-1-h, -1) \times (0,1)\), see “Appendix A” for the details of this computation. We eventually obtain

The error estimator terms related to \(K_2^h\) are

where we exploit \(\varDelta u_h = 0\) and \(p=1\). Since \(h_{K_2^h} = \sqrt{(1+h)^2 + 1} \in (\sqrt{2}, \frac{\sqrt{13}}{2}]\) 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.

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 h-adaptive 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 \(\eta \) from (13). In particular, we study the efficiency index \({\text {eff}} :=\eta / \left\Vert \nabla u - \nabla u_h\right\Vert _{\varOmega }\). To study the local overestimation and the efficiency of the error estimator, we introduce a (maximum) local efficiency index \({\mathrm {loc}} :=\max _{K \in \mathscr {T}_h} \eta _K / \left\Vert \nabla u - \nabla u_h\right\Vert _{K \cap \varOmega }\). 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.

Domain \(\varOmega \) resulting from the removal of 11 polygonal holes from the L-shaped domain \(\widehat{\varOmega }\) with initial mesh consisting of three quad elements

Full size image

L-shaped domain problem: h-uniform refinements

Full size image

L-shaped domain problem: overestimation of the error by the estimator

Full size image

4.1 L-Shape with Polygonal Holes

The first example is given by Poisson’s problem on a domain \(\varOmega \) with several holes. It is obtained by removing 11 polygonal holes \(\omega _1, \dots , \omega _{11}\) of various sizes and shapes from the standard L-shaped domain, which serves as the enclosing domain \({\hat{\varOmega }}\), see Fig. 4. In accordance with the compatibility requirements on \(\varGamma _D\) and \(\varGamma _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\in V\) such that

for all \(v \in V\) with Dirichlet boundary \(\varGamma _D :=({\left\{ 0 \right\} } \times [0,1]) \cup ([0,1] \times {\left\{ 0 \right\} })\) and Neumann boundaries \({\bar{\varGamma }}_N :=\partial {\hat{\varOmega }} \setminus \varGamma _D\) and \({\hat{\varGamma }}_N :=\bigcup _{i=1}^{11} \partial \omega _i\). The data g is derived from the nonsmooth solution \(u(r,\varphi ) :=r^{2/3}{(\sin (2\varphi - \pi )/3)}\) given in polar coordinates, which yields a corner singularity in the origin. The FCM discretization consists in seeking \(u_h \in V_h\) such that

for all \(v_h \in V_h\) where \(\varepsilon :=10^{-12}\) and \(V_h\) is a FEM space based on a mesh \(\mathscr {T}_h\) of \({\hat{\varOmega }}\). Let the set of cut cells be defined as \(\tilde{\mathscr {T}_h} :={\left\{ K ;\ K \in \mathscr {T}_h, \, K \cap \varOmega \ne K \right\} }\). To obtain a suitable quadrature on a cut \(K \cap \varOmega \) with \(K \in {\tilde{\mathscr {T}}}_{h}\), we perform a Delaunay triangulation of \(K \cap \varOmega \) 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 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(\mathrm {DOF}^{-1/3}) = O(h^{2/3})\) is attained for \(p = 1\) as well as \(p = 2\) which is expected due to \(u\in H^{5/3-\varepsilon }(\varOmega )\), \(\varepsilon >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 in 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 \(\varOmega \). Indeed, the exact error on \(K \cap \varOmega \) is approx. \(2.7 \times 10^{-7}\) and the estimation is \(2.1 \times 10^{-5}\). However, the average estimator value per element is computed to be approx. \(1.6 \times 10^{-4}\) meaning that the estimation on \(K \cap \varOmega \) 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 p-refinements only. Here, we make use of a spline-based finite element space \(V_h \subset C^1\) described in [23]. The differentiable discrete functions of \(V_h\) have the advantage that the jump terms \(\left\Vert \left[ \partial _n u_h\right] \right\Vert _{e \cap \varOmega }\) 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 \cap \varOmega \) which are actually not needed during the solution of the discrete problem. Figure 7a shows the errors and estimates resulting from \(p = 3, \dots , 20\) indicating that the convergence order of \(O(\mathrm {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 \({\text {eff}}\) and \({\text {loc}}\) exhibit a p-dependent growth. Therefore, we also plot \({\text {eff}}/p\), \({\text {loc}}/p\). Since the scaled indices appear to be constant, we conclude that the multiplicative factor of O(p) does occur in this example.

L-shaped domain problem: p-uniform refinements

Full size image

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 multi-level 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 (-b \mathrm {DOF}^{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 \cap \varOmega \) only varies during the first few refinements where the new elements intersect the holes.

L-shaped domain problem: geometric refinements

Full size image

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(\mathrm {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-shaped domain problem: h-adaptive refinements, meshes for \(p = 1, 2\)

Full size image

L-shaped domain problem: h-adaptive refinements

Full size image

4.2 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 \(\varepsilon \) on the estimator. To this end, the problem setup from Sect. 4.1 is used, where the polygonal holes are replaced by circular holes \(\omega _1, \dots , \omega _{11}\). The domain \(\varOmega \), the holes, and an initial mesh are displayed in Fig. 11.

Recall that, in the FCM, the problem of constructing a mesh for \(\varOmega \) is shifted to the numerical integration. In particular, the domain-approximation error is negligibly small if integrals on \(\varOmega \) are computed with machine precision. For the computation of the integrals on a cut \(K \cap \varOmega \) for \(K \in {\tilde{\mathscr {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 is constructed such that a set of polynomial basis functions \(\varphi _1, \dots , \varphi _m\) is integrated exactly. This is done by solving the linear system

Domain \(\varOmega \) resulting from the removal of 11 circular holes from the L-shaped domain \(\widehat{\varOmega }\) with initial mesh consisting of three quad elements

Full size image

for the unknown coefficients \(\alpha _1,\ldots ,\alpha _m\) with some fixed (transformed) Gauss points \(x_1, \dots \), \(x_m\) \( \in \) K. The unknown integral on the right-hand side can be computed by applying integration by parts, which gives

for a function \(\varPhi _i : \mathbb {R}^d \rightarrow \mathbb {R}^d\) with \({\text {div}} \varPhi _i = \varphi _i\). For a tensor-product polynomial function \(\varphi _i\) with \(\varphi _i(x) = \prod _{k=1}^d \varphi _{i,k}(x_k)\), the function \(\varPhi _i\) can easily be obtained by computing the antiderivative of one of the functions \(\varphi _{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 \(\int _{t_1}^{t_2} \varphi _{i}(\cos (t), \sin (t)) \,\mathrm {d}t\) 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 \(\varphi _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 \(\varepsilon \), we discuss four configurations with the parameter \(\varepsilon \in {\left\{ 10^{-1}, 10^{-2}, 10^{-4}, 10^{-6} \right\} }\). The decay of the error and the estimator is displayed in Fig. 12.

The upper bound \(\left\Vert \nabla u - \nabla u_h\right\Vert _{\varOmega } \lesssim \eta + \varepsilon ^{1/2}\) from (13) suggests that the error may be as large as \(O(\varepsilon ^{1/2})\). Thus, the error may become stationary after reaching the magnitude of \(O(\varepsilon ^{1/2})\), i.e., no further reduction of the error can be achieved even if the \({\text {DOF}}\) are increased. However, in this example, we notice that the error becomes stationary when \(\left\Vert \nabla u - \nabla u_h\right\Vert _{\varOmega } \approx O(\varepsilon )\). This is considerably better than the worst-case bound \(O(\varepsilon ^{1/2})\) that results from the use of the bound \(\varepsilon \left\Vert \nabla u_h\right\Vert _{{\hat{\varOmega }} \setminus \varOmega } \lesssim \varepsilon ^{1/2}{(\left\Vert f\right\Vert _{\varOmega } + \left\Vert g\right\Vert _{\varGamma _N})}\) in (21). Here, \(f = 0\) gives

which explains the improved behavior in this example.

The experiment confirms that \(\varepsilon \) should be small in order to obtain a sufficiently accurate solution. In fact, the errors are much larger for \(\varepsilon \in {\left\{ 10^{-1}, 10^{-2} \right\} }\) than for \(\varepsilon \in {\left\{ 10^{-4}, 10^{-6} \right\} }\), see Fig. 12. In the latter case, the errors are even quite similar (until the error for \(\varepsilon = 10^{-4}\) becomes stationary). We see that the estimator has approximately the same value for each \(\varepsilon \) until its decay rate reduces. The dependencies of the error and the estimator on \(\varepsilon \) 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 \(\varepsilon \). This may also be seen in Fig. 13 where the efficiency indices are depicted.

L-shaped domain problem with circular holes: decay of error and estimate for h-adaptive refinements (\(p = 10\))

Full size image

L-shaped domain problem with circular holes: efficiency indices (\(p = 10\))

Full size image

4.3 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\) and take it as the embedding domain, i.e. \({\hat{\varOmega }} :=(-1,1)^3 \setminus [0,1]^3\). The domain of interest \(\varOmega \) is

where each \(\omega _i\) is a polyhedral hole with a triangular surface in one of the seven subcubes \({\hat{\varOmega }}\) of equal size, see Fig. 14.

Initial mesh \(\mathscr {T}_1\) of seven cubes. \(\varOmega \) is obtained by removing seven polyhedrons from Fichera’s corner

Full size image

We solve Poisson’s problem with homogeneous Dirichlet conditions and inhomogeneous Neumann conditions

where \(\varGamma _D :={({\left\{ 0 \right\} } \times [0,1]^2)} \cup {([0,1] \times {\left\{ 0 \right\} } \times [0,1])} \cup {([0,1]^2 \times {\left\{ 0 \right\} })}\) are the faces of Fichera’s corner at the origin and \(\varGamma _N :=\partial \varOmega \setminus \varGamma _D\) is the union of the remaining faces and all boundaries of the holes \(\omega _i\). The weak formulation of the problem consists in finding \(u\in V\) such that

for all \(v \in V\), where the data \(\widetilde{u}\) is derived from the superposition of three solutions of the L-shape domain problem as suggested in [6, 25], i.e.

with the polar coordinates \((r(w_1, w_2), \varphi (w_1,w_2))\in \mathbb {R}_+\times [0,2\pi )\) of \((w_1, w_2) \in \mathbb {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 \in V_h\) such that

for all \(v_h \in V_h\), 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 \(\varepsilon :=10^{-12}\). The initial coarse mesh consists of the seven cubes as introduced above.

In order to compute the integrals on a cut \(K \cap \varOmega \) for some \(K \in {\tilde{\mathscr {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 \cap \varOmega \). 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(\mathrm {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.

Estimate of the energy error \(\left\Vert \nabla u - \nabla u_h\right\Vert _{\varOmega }\)

Full size image

Adaptive mesh and a parallel cut (elements of degree \(p = 1\))

Full size image