A Posteriori Error Estimates with Boundary Correction for a Cut Finite Element Method

TL;DR

This paper develops residual-based a posteriori error estimates for the CutFEM method applied to elliptic problems with complex boundaries, ensuring robustness and accounting for boundary geometry and data approximation.

Contribution

It introduces a new a posteriori error estimation framework for CutFEM that is reliable, efficient, and robust to boundary cuts, including boundary geometry and data approximation.

Findings

Reliability and efficiency of the error estimator are theoretically proven.

Constants in the estimates are robust with respect to boundary cuts.

The method effectively handles non-polygonal boundaries in elliptic problems.

Abstract

In this work we study a residual based a posteriori error estimation for the CutFEM method applied to an elliptic model problem. We consider the problem with non-polygonal boundary and the analysis takes into account the geometry and data approximation on the boundary. The reliability and efficiency are theoretically proved. Moreover, constants are robust with respect to how the domain boundary cuts the mesh.

Equations

$$-\Delta u$$

$$u$$

$$\Vert u{\Vert }_{H^{s+2}(\Omega )}\lesssim \Vert f{\Vert }_{H^s(\Omega )},-1\le s\le s_0$$

$${\mathcal{T}}_h(\omega )=\{K\in {\mathcal{T}}_{0,h}:\overline{K}\cap \overline{\omega }\ne \varnothing \}$$

$${△}_h(\omega )=\bigcup _{K\in {\mathcal{T}}_h(\omega )}K$$

$${\mathcal{T}}_h:={\mathcal{T}}_h(\Omega \cup {\Omega }_h)$$

$${△}_h:={△}_h(\Omega \cup {\Omega }_h)$$

$${\mathcal{T}}_h^b=\{K\in {\mathcal{T}}_{0,h}:K\cap \partial {\Omega }_h\ne \varnothing \}\subset {\mathcal{T}}_h.$$

$$\Vert {\varrho }_h{\Vert }_{L^{\infty }(\partial {\Omega }_h\cap K)}\le O(h_K)\forall K\in {\mathcal{T}}_h^b.$$

$$\Vert \varrho {\Vert }_{L^{\infty }(\partial {\Omega }_h\cap K)}=O(h_K)\forall K\in {\mathcal{T}}_h^b.$$

$${\mathcal{F}}_h=\{F\in {\mathcal{E}}_I:(K_F\cup K_F^{\prime })\cap \partial {\Omega }_h\ne \varnothing \}$$

$$[[v]]{|}_F:=v_F^{+}-v_F^{-}\text{mbox}and{v}_F^{\pm }(x)=\underset{s\to 0^{+}}{\lim }v(\mathbf{x}\mp s{\mathbf{n}}_F).$$

$$V_h:=\{v\in H^1({△}_h):v{|}_K\in {\mathbb{P}}_1(K)\forall K\in {\mathcal{T}}_h\}$$

$$\begin{array}{rl}{a}_0(v,w)& :=(\nabla v,\nabla w{)}_{{\Omega }_h}-{⟨D_{{\mathbf{n}}_h}v,w⟩}_{\partial {\Omega }_h}-{⟨v,D_{{\mathbf{n}}_h}w⟩}_{\partial {\Omega }_h}+\sum _{K\in {\mathcal{T}}_h^b}\frac{\beta }{h_K}{⟨v,w⟩}_{{\Gamma }_K},\\ j_h(v,w)& :=\gamma \sum _{F\in {\mathcal{F}}_h}{h}_F{⟨[[D_{{\mathbf{n}}_F}v]],[[D_{{\mathbf{n}}_F}w]]⟩}_F,\\ a_h(v,w)& :=a_0(v,w)+j_h(v,w),\\ l_h(w)& :=(f,w{)}_{{\Omega }_h}-{⟨g_h,D_{{\mathbf{n}}_h}w⟩}_{\partial {\Omega }_h}+\sum _{K\in {\mathcal{T}}_h^b}\frac{\beta }{h_K}{⟨g_h,w⟩}_{{\Gamma }_K},\end{array}$$

$$a_h(u_h,v)=l_h(v)\forall v\in V_h$$

$$|\text{mspace}-1.0mu|\text{mspace}-1.0mu|v|\text{mspace}-1.0mu|\text{mspace}-1.0mu{|}^2$$

$$|\text{mspace}-1.0mu|\text{mspace}-1.0mu|v|\text{mspace}-1.0mu|\text{mspace}-1.0mu{|}_h^2$$

$$\Vert v{\Vert }_{\partial K}\le C(h_K^{-1/2}\Vert v{\Vert }_K+h_K^{1/2}\Vert \nabla v{\Vert }_K)\forall v\in H^1(K),$$

$$h_K^{-\frac{1}{2}}\Vert v_h{\Vert }_{\partial K}+\Vert \nabla v_h{\Vert }_K\le C{h}_K^{-1}\Vert v_h{\Vert }_K\forall v_h\in {\mathbb{P}}_1(K)$$

$$\Vert v_h{\Vert }_{{\Gamma }_K}\le C(h_K^{-1/2}\Vert v_h{\Vert }_K+h_K^{1/2}\Vert \nabla v_h{\Vert }_K)\forall v_h\in {\mathbb{P}}_1(K)$$

$$\Vert v{\Vert }_K\lesssim h_K\Vert \nabla v{\Vert }_{{\mathcal{S}}_K},$$

$$\text{mbox}dist(K,\partial \Omega )\le O(h)$$

$$\tilde{e}\in H^1(\Omega )\text{mbox}and\tilde{e}{|}_{\partial \Omega }=(g-u_h)$$

$${\eta }_K^2=h_K^2\Vert f{\Vert }_{K\cap {\Omega }_h}^2+\sum _{F\in {\mathcal{E}}_K\cap {\mathcal{E}}_I}\frac{h_F}{2}\Vert [[D_{{\mathbf{n}}_F}{u}_h]]{\Vert }_F^2+h_K^{-1}\Vert g_h-u_h{\Vert }_{L^2({\Gamma }_K)}^2+\Vert \nabla \tilde{e}{\Vert }_{K\cap \Omega }^2.$$

$$\eta ={\left(\sum _{K\in {\mathcal{T}}_h}{\eta }_K^2\right)}^{1/2}.$$

$$\Vert \nabla e{\Vert }_{\Omega }\le C_r{\left(\sum _{K\in {\mathcal{T}}_h}{\eta }_K^2+\sum _{K\in {\mathcal{T}}_h}{h}_K^2\Vert f{\Vert }_{(\Omega \setminus {\Omega }_h)\cap K}^2\right)}^{1/2}$$

$$\begin{array}{rl}\Vert \nabla e{\Vert }_{\Omega }& \le \Vert \nabla (e-\tilde{e}){\Vert }_{\Omega }+\Vert \nabla \tilde{e}{\Vert }_{\Omega }\\ & =\underset{v\in H_0^1(\Omega )}{\sup }\frac{(\nabla (e-\tilde{e}),\nabla v{)}_{\Omega }}{\Vert \nabla v{\Vert }_{\Omega }}+\Vert \nabla \tilde{e}{\Vert }_{\Omega }\\ & \le \underset{v\in H_0^1(\Omega )}{\sup }\frac{(\nabla e,\nabla v{)}_{\Omega }}{\Vert \nabla v{\Vert }_{\Omega }}+2\Vert \nabla \tilde{e}{\Vert }_{\Omega }.\end{array}$$

$$\underset{v\in H_0^1(\Omega )}{\sup }\frac{(\nabla e,\nabla v{)}_{\Omega }}{\Vert \nabla v{\Vert }_{\Omega }}\lesssim {\left(\sum _{K\in {\mathcal{T}}_h}{\eta }_K^2+\sum _{K\in {\mathcal{T}}_h}{h}_K^2\Vert f{\Vert }_{(\Omega \setminus {\Omega }_h)\cap K}^2\right)}^{1/2},$$

$$\Vert \varrho {\Vert }_{L^{\infty }({\Gamma }_K)}\le o(h_K)\forall K\in {\mathcal{T}}_h^b.$$

$$(\nabla e,\nabla v{)}_{\Omega }={\mathcal{A}}_1+{\mathcal{A}}_2$$

$${\mathcal{A}}_1=(f,v-v_h{)}_{{\Omega }_h}+(f,v{)}_{\Omega \setminus {\Omega }_h}$$

Full text

{AMS}

68Q25, 68R10, 68U05

1 Introduction

Meshing and re-meshing procedures remain an important challenge in computational methods, since they can be very computationally expensive for the finite element methods, particularly for problems with complex geometries and interfaces that move during the computational process. Several methods have been invented to alleviating the meshing procedure of the exact domain. Fictitious domain methods were introduced in the finite element context in the papers by Glowinski et al. Glowinski & Pan, (1992) and unfitted finite element methods in the works by Barrett and Elliott Barrett & Elliott, (1984).

Since these seminal works many approaches have been suggested on how to integrate the geometry data in finite element computations in a way that reduces the meshing effort. For instance the fat boundary method by Bertoluzza et al. Bertoluzza *et al. *, (2005), the fictitious domain methods using Lagrange multipliers inspired by extended finite element methods, pioneered by Haslinger and Renard Haslinger & Renard, (2009), Burman & Hansbo, (2010), or the cut finite element method using Nitsche’s method, introduced by Hansbo and Hansbo and further developed in the fictitious domain framework by various authors Hansbo & Hansbo, (2002), Burman, (2010), Embar *et al. *, (2010), Becker *et al. *, (2011), Burman & Hansbo, (2012).

Cut finite element methods have been extensively studied in the recent years, both in the context of interface problems and for fictitious domain methods Burman *et al. *, (2015b). We restrict the discussion herein to the fictitious domain case. The methodology alleviates the meshing process by employing a background mesh that can be highly structured and letting the domain boundary cut through the elements of the mesh. In cells cut by the boundary the equations are integrated only on the intersection of the element with the physical domain. The boundary conditions are then imposed either using Lagrange multipliers or Nitsche’s method and stability with respect to the cut may be ensured using a ghost penalty stabilization Burman, (2010). The method was originally designed for continuous approximation spaces, but has been adapted for various other methods, e.g, Discontinuous Galerkin method Johansson & Larson, (2013), Gürkan & Massing, (2019), Hybrid High order Methods Burman & Ern, (2019) or Isogeometric analysis Elfverson *et al. *, (2018). It has been applied to a number of different partial differential equations, e.g, incompressible elasticity/Stokes’ equations Burman & Hansbo, (2014), Massing *et al. *, (2014b), Burman *et al. *, (2015a), Guzmán & Olshanskii, (2018), linear elasticity Hansbo *et al. *, (2017), Helmholtz equations Swift, (2018), time dependent parabolic problems on moving domains Hansbo *et al. *, (2015), Oseen’s problem Massing *et al. *, (2018), Winter *et al. *, (2018) and other fluid models Schott & Wall, (2014), Schott *et al. *, (2016). It has also been applied successfully to shape optimization problems Burman *et al. *, (2017, 2018b), Bernland *et al. *, (2018) and other advanced engineering applications Bui *et al. *, (2019).

Typically, in cut finite element methods, a domain with curved boundary is approximated using a piecewise affine boundary approximation. This gives a sufficiently good geometry approximation for piecewise affine elements, but if higher order elements are used the geometry approximation must be improved. To alleviate the integration problem resulting from the elements cut by curved boundaries, isoparametric techniques Lehrenfeld, (2017) and so called boundary value correction techniques have been proposed Burman *et al. *, (2018a), Boiveau *et al. *, (2018). For both these cases optimal order a priori error estimates for arbitrary order of the polynomial approximation have been derived. For an analysis of the fitted finite element method on domains with curved boundaries we refer to Bramble & King, (1994).

The purpose of this paper is to design an a posteriori error estimation for the cut finite element method for problems with non-polygonal boundary. We will base our discussion on the cut finite element methods for the Poisson problem introduced by Burman and Hansbo Burman & Hansbo, (2012). This method uses Nitsche’s method Nitsche, (1971) to impose Dirichlet boundary conditions and a ghost penalty term Burman, (2010) to enhance stability in the boundary zone. We restrict the discussion to the fictitious domain problem and piecewise affine approximation. The main motivation for studying the a posteriori error estimation is for the application of adaptive mesh refinement (AMR) procedure. It is well known that AMR is extremely useful for problems with singularities, discontinuities, sharp derivatives, etc. and it has been extensively studied in the last two decades, Verfürth, (1994), Ainsworth & Oden, (2011). However, there is very limited work in the literature that takes the geometry approximation into account. In Dörfler & Rumpf, (1998) the a posteriori error estimation is studied for the conforming finite element method on curved boundary where the boundary vertices of the approximated mesh must be located on the true boundary and it is assumed that data can be requested at any point on the boundary and inside the domain, i.e., there is only approximation of the geometry. In Ainsworth & Rankin, (2017) a fully computable error bound is provided for the conforming linear elements with pure Neumann data. To the author’s knowledge there is no existing work on the a posteriori error estimation for cut finite element methods on problems with curved boundary, to allow for. Another approach for the handling of geometric singularities in the cutFEM framework was recently proposed in Jonsson *et al. *, (2019).

In our error analysis we do not assume that the vertices of the approximation boundary are located on the boundary of the continuous problem, while we do require that the discrete boundary is a sufficiently good approximation of the boundary of the continuous problem. We do not require that data inside the domain can be everywhere requested but only those originally provided in the continuous setting. For the formulation it is necessary to extend the source term data from the domain of the continuous problem to the computational domain, the extension must satisfy certain stability properties. The error estimator comprises two parts. One part is due to the numerical approximation; and the other part is exclusively due to the boundary approximation. We refer to this latter contribution as the boundary correction error. This part must be computed using a locally improved boundary approximation in the boundary zone. The computation of this correction is discussed in the numerical section.

Unlike the methods whose mesh is an exact partition of the computational domain, one major challenge for cut finite element is that the approximated domain cuts the background mesh in an arbitrary way. Because of this the classical efficiency analysis, i.e., by applying the local elementary bubble functions, is not robust. Therefore we divide the efficiency part of the estimates in two parts. The residuals in the bulk, away from the boundary, are estimated in the ordinary fashion. For the residuals in the cut elements and the nonconforming ghost penalty operator, we instead prove a global lower bound of best approximation type, showing that the boundary residuals are bounded by the best approximation of $u$ in the physical domain and oscillations.

This paper is organized as follows. In Section 2 the model problem and the cut finite element method is introduced. The error estimator is introduced in Section 3 and as well its global reliability. The local efficiency is proved in Section 4. Finally, we show the results of several numerical experiment in Section 5.

2 Model Problem and the Cut Finite Element Method

2.1 The Continuous Problem

Let $\Omega$ be a domain in $\mathbb{R}^{d}$ with Lipschitz continuous, piecewise smooth boundary $\partial\Omega$ and exterior unit normal $\boldsymbol{n}$.

We consider the problem: find $u:\Omega\rightarrow\mathbb{R}$ such that

[TABLE]

where $f\in L^{2}(\Omega)$ and $g\in H^{1/2}(\partial\Omega)$ are given data. It follows from the Lax-Milgram lemma that there exists a unique solution $u\in H^{1}(\Omega)$ to this problem. We also have the following elliptic regularity estimate

[TABLE]

for some $s_{0}\geq 1/2$ depending on the domain. Here and below we use the notation $\lesssim$ to denote less or equal up to a generic constant that is independent of the mesh-geometry configuration.

2.2 The Mesh, Discrete Domains, and Finite Element Spaces

Assume that $\partial\Omega$ is composed of a finite number of smooth surfaces $\Gamma_{i}$, such that $\partial\Omega=\cup_{i}\bar{\Gamma}_{i}$. We let $\rho$ be the signed distance function, negative on the inside and positive on the outside, to $\partial\Omega$ and we let $U_{\delta}(\partial\Omega)$, for $\delta>0$, be the tubular neighborhood $\{\boldsymbol{x}\in\mathbb{R}^{d}:|\rho(\boldsymbol{x})|<\delta\}$ of $\partial\Omega$. Let $\Omega_{0}\subset\mathbb{R}^{d}$ be a polygonal domain such that $U_{\delta_{0}}(\Omega)\subset\Omega_{0}$ where $\delta_{0}$ is chosen such that $\rho$ is well defined in $U_{\delta_{0}}(\Omega)$. Let ${\mathcal{T}}_{0,h}$ be a partition of $\Omega_{0}$ into shape regular triangles or tetrahedra. Note that this setting allows meshes with locally dense refinement.

Given a subset $\omega$ of $\Omega_{0}$, let ${\mathcal{T}}_{h}(\omega)$ be the submesh defined by

[TABLE]

i.e., the submesh consisting of elements that intersect $\overline{\omega}$, and let

[TABLE]

be the union of all elements in ${\mathcal{T}}_{h}(\omega)$.

For each ${\mathcal{T}}_{0,h}$ let $\Omega_{h}$ be a polygonal domain approximating $\Omega$ and we assume that $\partial\Omega_{h}\subset U_{\delta_{0}}(\partial\Omega)$ which implies that $\partial\Omega_{h}$ is within the distance of $\delta_{0}$ to $\partial\Omega$. We assume neither $\Omega_{h}\subset\Omega$ nor $\Omega\subset\Omega_{h}$, instead the maximum distance between the two domains has to be small enough.

Let the active mesh be defined by

[TABLE]

i.e., the submesh consisting of elements that intersect $\Omega_{h}\cup\Omega$, and let

[TABLE]

be the union of all elements in ${\mathcal{T}}_{h}$. Since $\partial{\Omega}_{h}$ does not necessarily fit the mesh we denote by ${\mathcal{T}}_{h}^{b}$ the set of elements that are “cut” by $\partial{\Omega}_{h}$,

[TABLE]

We assume that $\Omega_{h}$ is constructed in such a way that for each $K\in{\mathcal{T}}_{h}^{b}$, the intersection $K\cap\partial\Omega_{h}$ is a subset of a $d-1$ dimensional hyperplane, i.e., a line segment in two dimensions and a subset of a plane in three dimensions. Under the smallness assumption of $h$, for any element $K\in{\mathcal{T}}_{h}^{b}$ there exists an element $K^{\prime}\in{\mathcal{T}}_{h}(\Omega_{h})\setminus{\mathcal{T}}_{h}^{b}$ such that $\mbox{dist}(K,K^{\prime})=O(h)$ where $O(\cdot)$ denotes the big ordo. This is always possible for small enough $h$ since $\Omega$ is Lipschitz.

Also we denote by ${\mathcal{E}}$ the set of all facets on ${\mathcal{T}}_{h}$ and by ${\mathcal{E}}_{I}\subset{\mathcal{E}}$ the set of all interior facets with respect to ${\mathcal{T}}_{h}$. For each $F\in{\mathcal{E}}$ denote by $\boldsymbol{n}_{F}$ an unit vector normal to $F$ and by $h_{F}$ the diameter of $F$. For each $K\in{\mathcal{T}}_{h}$ denote by $h_{K}$ the diameter of $K$ and by ${\mathcal{E}}_{K}$ the set of all facets of $K$.

On the boundary $\partial{\Omega}_{h}$ let $\boldsymbol{n}_{h}$ be the outer normal to $\partial{\Omega}_{h}$. For each $\Omega_{h}$ we assume that, for $\delta_{0}$ small enough, there exist functions $\boldsymbol{\nu}_{h}:\partial\Omega_{h}\rightarrow\mathbb{R}^{d}$, $|\boldsymbol{\nu}_{h}|=1$, and $\varrho_{h}:\partial\Omega_{h}\rightarrow\mathbb{R}$, on $U_{\delta_{0}}(\partial\Omega)$, such that the function $\boldsymbol{p}_{h}(\boldsymbol{x},\varsigma):=\boldsymbol{x}+\varsigma\boldsymbol{\nu}_{h}(\boldsymbol{x})$, is well defined and satisfies $\boldsymbol{p}_{h}(\boldsymbol{x},\varrho_{h}(\boldsymbol{x}))\in\partial\Omega$ for all $\boldsymbol{x}\in\partial\Omega_{h}$. The existence of the vector-valued function $\boldsymbol{\nu}_{h}$ is known to hold on Lipschitz domains, see Grisvard Grisvard, (2011).

We further assume that $\boldsymbol{p}_{h}(\boldsymbol{x},\varsigma)\in U_{\delta_{0}}(\Omega)$ for all $\boldsymbol{x}\in\partial\Omega_{h}$ and all $\varsigma$ between [math] and $\varrho_{h}(\boldsymbol{x})$. For conciseness we will drop the second argument, $\varsigma$, of $\boldsymbol{p}_{h}$ below whenever it takes the value $\varrho_{h}(\boldsymbol{x})$ and let $\boldsymbol{p}_{h}$ denotes the map $\boldsymbol{p}_{h}:\partial\Omega_{h}\rightarrow\partial\Omega$. Moreover, we assume that the following assumption is satisfied

[TABLE]

The above assumption immediately implies that $\partial{\Omega}_{h}$ is within the distance of $O(h)$ of $\partial{\Omega}$. More precisely,

[TABLE]

This assumption is necessary for the constant in the a posteriori error estimates to be independent of the geometry/mesh configuration. It is however not enough to guarantee optimal a priori error estimates, which requires $\|\varrho\|_{L^{\infty}(\partial\Omega_{h}\cap K)}\leq O(h_{K}^{2})$, and $\|\boldsymbol{n}-\boldsymbol{n}_{h}\|_{L^{\infty}(\partial\Omega_{h}\cap K)}\leq O(h_{K})$, see Burman *et al. *, (2018a).

2.3 The Cut Finite Element Method

In this subsection we recall the CutFEM introduced in Burman & Hansbo, (2012). We begin with some necessary notation. Let $\mathcal{F}_{h}$ be the set of faces intersecting the approximate boundary ${\Omega}_{h}$:

[TABLE]

where $K_{F}$ and $K^{\prime}_{F}$ are those two elements sharing $F$ as a common facet, and for any discontinuous function $v$ define its jump on the facet $F$ by

[TABLE]

Next define the finite element space

[TABLE]

and the forms

[TABLE]

where $D_{\boldsymbol{n}_{h}}v:=\boldsymbol{n}_{h}\cdot\nabla v$, $\gamma$ and $\beta$ are positive constants, $g_{h}$ is an approximation of $g$ defined on $\partial\Omega_{h}$, typically $g_{h}(\boldsymbol{x})=g\circ\boldsymbol{p}_{h}$, and $\Gamma_{K}=K\cap\partial{\Omega}_{h}$ and $f|_{\Omega_{h}\setminus\Omega}$ is defined by some suitable extension (in the low order case considered here it can be taken to be zero).

Remark 2.1**.**

*The stabilizing term $j_{h}(v,w)$, which is the so called ghost penalty term, is introduced to extend the coercivity of $a_{h}(\cdot,\cdot)$ to all of $\triangle_{h}$, see Burman, (2010), Massing *et al. *, (2014a). Thanks to this property one may prove that the condition number of the linear system is uniformly bounded independent of how $\Omega_{h}$ is oriented compared to the mesh. *

Remark 2.2**.**

*In order to guarantee the coercivity of Eq. 2.11 $\beta$ has to be chosen large enough. *

The finite element method is then to find $u_{h}\in V_{h}$ such that

[TABLE]

where $a_{h}$ and $l_{h}$ are defined in Eq. 2.11.

For $v\in H^{1}({\Omega}_{0})$ define the continuous and discrete energy norm respectively by

[TABLE]

and

[TABLE]

2.4 Some Important Inequalities

We have the following well known inverse and trace inequalities (Di Pietro & Ern,, 2012, Section 1.4.3),

[TABLE]

[TABLE]

and Hansbo & Hansbo, (2002)

[TABLE]

where the constant $C$ is independent of the relative location of ${\Omega}_{h}$.

Lemma 1**.**

Let $v\in H_{0}^{1}({\Omega})$. Then for any $K$ such that $K\in{\mathcal{T}}_{h}^{b}$ or $K\cap({\Omega}\setminus{\Omega}_{h})\neq\emptyset$ there exists a local convex neighborhood $\mathcal{S}_{K}$ of $K$ such that $v$ vanishes on a nonzero subset of $\partial S_{K}$ and

[TABLE]

*where we defined $v$ outside ${\Omega}$ using the trivial extension $v|_{{\Omega}^{c}}=0$. *

Proof 2.3**.**

If $K\subset{\Omega}^{c}$, the estimate Eq. 2.18 is obvious since $v$ is taken uniformly zero outside ${\Omega}$. If $K\cap\partial{\Omega}\neq\emptyset$, then we let $x\in K\cap\partial{\Omega}$ and let $S_{K}=B_{\delta}(x)$ be the open ball centered at $x$ with radius $\delta\sim h$. The estimate (2.18) is a direct consequence of the Poincaré inequality on $S_{K}$. Otherwise $K\subset{\Omega}$. From Eq. 2.8 it follows that

[TABLE]

*and thus there is $x\in\partial\Omega$ and $\delta\sim h$ such that $K\subset B_{\delta}(x)$ and now Eq. 2.18 follows with $S_{K}=B_{\delta}(x)$ from the Poincaré inequality on $S_{K}$. This completes the proof of the lemma. *

3 Global Reliability

Let $e:=u-u_{h}$ and $\tilde{e}$ be any function such that

[TABLE]

where $u_{h}$ is the solution to Eq. 2.12.

For each element $K\in{\mathcal{T}}_{h}$, define the local element error indicator $\eta_{K}$ by

[TABLE]

and let the global error estimator be defined by

[TABLE]

Theorem 1**.**

Let $\tilde{e}\in H^{1}({\Omega})$ such that $\tilde{e}|_{\partial{\Omega}}=(g-u_{h})|_{\partial{\Omega}}$. We have the following reliability bound

[TABLE]

*where the constant $C_{r}$ does not depend on the location of domain-mesh intersection nor the mesh size. *

Proof 3.1**.**

Note that $e-\tilde{e}\in H_{0}^{1}({\Omega})$. By the triangle inequality we have that

[TABLE]

It then suffices to show that

[TABLE]

*which is a direct result from LABEL:{lem:error-rep} and 4 below. This completes the proof of the theorem. *

Remark 2**.**

We omit the second term in Eq. 3.3 in the algorithm computation in order to avoid integration on the curved domain. What’s more, one can easily make it a higher order term by satisfying

[TABLE]

Lemma 3**.**

For any $v\in H_{0}^{1}({\Omega})$ and $v_{h}\in V_{h}$ the following equality holds:

[TABLE]

where

[TABLE]

and

[TABLE]

Proof 3.2**.**

Extending $v\in H^{1}_{0}(\Omega)$ to $\mathbb{R}^{d}$ by setting $v=0$ outside of $\Omega$ and using integration by parts gives

[TABLE]

For the second term in Eq. 3.8, by Eq. 2.12 we have

[TABLE]

For the last two terms in (3.8) applying the integration by parts gives

[TABLE]

In the last equality we used the following identity

[TABLE]

*thanks to the fact that $v=0$ on $\partial{\Omega}$ and in ${\Omega}_{h}\setminus{\Omega}$, where $\boldsymbol{n}$ without subscript denotes the outer normal to the boundary of domain being integrated. The desired identity Eq. 3.5 is then a direct consequence of Eq. 3.8–Eq. 3.10. This completes the proof of the lemma. *

Remark 3.3**.**

*Note that $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ in 3 does not contain the inconsistency error caused by the boundary approximation. In other words, the boundary correction error has been independently isolated by the term $\|\nabla\tilde{e}\|$. Since any $\tilde{e}$ satisfying the boundary condition will yield an upper bound for the error, an inappropriate construction of $\tilde{e}$ will potentially over estimate the error. *

Lemma 4**.**

Let $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ be given in 3. Then we have the following estimates:

[TABLE]

and

[TABLE]

Proof 3.4**.**

The first term of $\mathcal{A}_{1}$ in (3.5) can be bounded directly using the Cauchy-Schwarz inequality and the approximation property of the Scott-Zhang interpolation Scott & Zhang, (1990),

[TABLE]

where $\triangle_{K}$ is the union of elements in ${\mathcal{T}}_{h}$ that shares at least one vertex with $K$.

By Eq. 2.18 the second term of $\mathcal{A}_{1}$ can now be bounded as follows:

[TABLE]

(3.11) is then a direct result of (3.13) and (3.14).

We now proceed to bound terms in $\mathcal{A}_{2}$. The first term of $\mathcal{A}_{2}$ can be bounded directly by applying the Cauchy-Schwarz inequality and the approximation property of the Scott-Zhang interpolation,

[TABLE]

To bound the second term in $\mathcal{A}_{2}$ we apply the trace inequality and Eq. 2.18 that

[TABLE]

The jump penalty term in $\mathcal{A}_{2}$ can be bounded by applying the Cauchy-Schwarz, triangle, trace, and the inverse inequalities and the stability of the interpolator:

[TABLE]

where $K_{F}^{+}(K_{F}^{-})$ is the element with $\boldsymbol{n}_{F}(-\boldsymbol{n}_{F})$ being its outer unit normal on $F$ and $\triangle_{F}:=\triangle_{K_{F}^{+}}\cup\triangle_{K_{F}^{-}}$.

The remaining two terms in $\mathcal{A}_{2}$ can be bounded by applying the Cauchy-Schwarz and trace inequalities, Eq. 2.18 and the stability of the interpolator:

[TABLE]

*Finally, combining (3.15)-(3.18) and the Young’s inequality yields (3.12). This completes the proof of the lemma. *

4 Efficiency

In this section we prove the efficiency of the error indicator introduced in Eq. 3.1. For the classical finite element method the efficiency is usually proved in a local fashion by applying the local facet and element bubble functions. For elements that are not intersected with the boundary and whose facets do not belong to the ghost penalty set we are able to prove the efficiency using the classical bubble technique.

Lemma 1**.**

Let $K$ be a given element in ${\mathcal{T}}_{h}$ such that $K\not\in{\mathcal{T}}_{h}^{b}$ and ${\mathcal{E}}_{K}\cap\mathcal{F}_{h}=\emptyset$. Then the following local efficiency result holds:

[TABLE]

*where $\tilde{\Delta}_{K}$ is a local neighborhood of $K$ and the efficiency constant $C_{e}$ does not depend on the mesh size nor the domain-mesh intersection. *

Proof 4.1**.**

*The proof is classical and we refer Verfürth, (1994), Cai *et al. *, (2017). *

Remark 4.2**.**

*For the regular elements, we take the boundary correction error $\|\nabla\tilde{e}\|_{K\cap{\Omega}}$ to be [math] since we can always design $\tilde{e}$ in such a way that it vanishes inside regular elements in order to avoid over-estimation (see Section 5.1). *

For an element $K$ that is irregular, i.e., $K\in{\mathcal{T}}_{h}^{b}$ or ${\mathcal{E}}_{K}\cap\mathcal{F}_{h}\neq\emptyset$, applying the same bubble technique for cut elements unfortunately will result in the dependence on the location of domain-mesh intersection. As an alternative we prove the efficiency for the term of ghost penalty as a whole. The efficiency bounds for the element residual on irregular elements can be then also proved with the aid of the efficiency result for the ghost penalty and the numerical scheme.

For each $K$ define $E_{K}=\partial K\cap((\Omega\cup\Omega_{h})\setminus(\Omega\cap\Omega_{h}))$ and $D_{K}=K\cap((\Omega\cup\Omega_{h})\setminus(\Omega\cap\Omega_{h}))$, see Fig. 1 for an illustration.

We first prove a discrete Poincare inequality that will be useful in the following efficiency proof.

Lemma 2**.**

For any $v\in V_{h}$ the following estimate is true:

[TABLE]

*where $|\mspace{-1.0mu}|\mspace{-1.0mu}|v_{h}|\mspace{-1.0mu}|\mspace{-1.0mu}|_{j_{h}}^{2}=j_{h}(v_{h},v_{h})$. *

Proof 4.3**.**

First, for each $K$ such that $E_{K},D_{K}\neq\emptyset$, note that using the trace inequality Eq. 2.16 we see that

[TABLE]

It follows that

[TABLE]

Now observe that there exists a set denoted by $S_{K}\subset\triangle_{h}$ such that $E_{K},D_{K}\subset S_{K}$, $\mbox{diam}(S_{K}\cap{\Omega}_{h})$ is of order $h_{K}$ and for which using Poincaré’s inequality

[TABLE]

By the equivalence of norms theorem we also have

[TABLE]

*Combining Eq. 4.3 and Eq. 4.4 gives Eq. 4.2. This completes the proof of the lemma. *

We assume that for each $K\in{\mathcal{T}}_{h}^{b}$ there exists at least one vertex, say $z$, of $K$ such that $\mbox{diam}({\omega}_{z}\cap{\Omega}_{h})=O(h_{K})$ where ${\omega}_{z}$ is the union of all all elements sharing $z$ as a vertex. Define

[TABLE]

where $f_{z}$ is chosen such that $f_{z}=\mathop{\mathrm{argmin}}_{c\in\mathcal{R}}h_{K}\|f-c\|_{{\omega}_{z}\cap{\Omega}_{h}}$. We will now proceed and prove a best approximation result for the element residuals of the cut elements. First we give a Lemma showing that the element residual of a cut element can be bounded by the ghost penalty term, the boundary residual and oscillation in data.

Lemma 4.4**.**

For each $K\in{\mathcal{T}}_{h}^{b}$ its element residual has the following bound:

[TABLE]

Proof 4.5**.**

By the triangle inequality we firstly have

[TABLE]

where $f_{z}$ at this point can be chosen as any arbitrary constant. Let $\lambda_{z}$ be the barycentric hat function associated with $z$ then we also have

[TABLE]

Let $w_{z}=f_{z}\lambda_{z}$. Applying Eq. 2.12 and integration by parts yields

[TABLE]

The last inequality utilizes the fact that $f_{z}$ is a constant and the trace and inverse inequalities. Combining Eq. 4.7, Eq. 4.8 and the Cauchy Schwarz inequality gives

[TABLE]

*which, combining with Eq. 4.6, yields Eq. 4.5. *

We now prove the main bound for the residuals in the cut elements.

Theorem 3**.**

Let $u$ and $u_{h}$ be the solution to Eq. 2.1 and Eq. 2.12, respectively. Then the following best approximation result holds:

[TABLE]

*where the constant $C_{e}$ does not depend on the mesh size nor on the domain-mesh intersection. *

Proof 4.6**.**

By the triangle inequality we have

[TABLE]

and

[TABLE]

Denote by $\phi=u_{h}-v_{h}$. Then

[TABLE]

*Applying the coercivity result (see Burman *et al. , (2015b)), Eq. 2.11 and Eq. 2.12 gives

[TABLE]

where

[TABLE]

To estimate $\mathcal{A}_{1}$, applying the Cauchy-Schwarz inequality and 2 gives

[TABLE]

$\mathcal{A}_{2}$* can be estimated directly using the Cauchy-Schwarz inequality,*

[TABLE]

To estimate $\mathcal{A}_{3}$, add and subtract suitable terms to obtain,

[TABLE]

Observe that using integration by parts we have

[TABLE]

where we used Lemma 2 for the last inequality. Then applying the Cauchy-Schwarz inequality and 2 gives

[TABLE]

Here we also used the inequality

[TABLE]

Collecting the bounds (4.13)–(4.16) we see that

[TABLE]

*Finally, combining Lemma 4.4, Eq. 4.11, Eq. 4.12 and Eq. 4.17 yields Eq. 4.10. This completes the proof of the theorem. *

Remark 4.7**.**

*Regarding the efficiency for the indicator of the boundary correction error, we have the following efficiency result: *

[TABLE]

*This also indicates that $\tilde{e}$ needs to be designed properly to avoid over estimation. In Section 5 we propose an approach to design and compute $\tilde{e}$. *

5 Numerical Results

In this section we present several numerical examples to validate the performance of the a posteriori error estimator in the adaptive mesh refinement procedure. The adaptive mesh refinement procedure is set as follows:

[TABLE]

For the penalty parameters in the finite element method, we set $\beta=10$ and $\gamma=0.1$. For the refinement strategy, we use the Dörfler marking strategy and the refine rate is set to be ten percent. Regarding the domain approximation, for a given $\phi$ being the level set function that satisfies $\phi=0$ on the boundary and negative (positive) inside (outside) the domain ${\Omega}$, let $\phi_{h}$ be the nodal interpolation of $\phi$ with respect to ${\mathcal{T}}_{h}$. And we define

[TABLE]

5.1 Computation of $\|\nabla\tilde{e}\|_{{\Omega}}$

In this subsection we introduce one way to compute the boundary correction error $\|\nabla\tilde{e}\|_{{\Omega}}$. We firstly construct a boundary correction mesh, denoted by ${\mathcal{T}}_{h}^{bc}$, which is finer than the current mesh ${\mathcal{T}}_{h}$ in order to more accurately track the boundary $\partial{\Omega}$. More precisely, the new mesh is built inside the union of all intersection elements, i.e.,

[TABLE]

For each element $K\in{\mathcal{T}}$, without loss of generality, in two dimensions, we assume that $\phi(z_{0})\leq\phi(z_{1})\leq\phi(z_{2})$ where $z_{i},i=0,1,2$ are vertices of $K$. For each $K\in{\mathcal{T}}_{h}^{b}$ we assume the following must be true: there exists at least one facet in $K$ such that the nodes connected to that facet have both positive and negative values, i.e.,

[TABLE]

We then further partition the elements in ${\mathcal{T}}_{h}^{b}$ based on the intersection of the domain and the mesh. The pseudo code for the algorithm is provided in Algorithm 1. An example boundary correction mesh for Example 1 is shown in Fig. 3.

With the mesh ${\mathcal{T}}^{bc}$ available, we define $\tilde{e}$ to be a piecewise linear conforming function with respect to ${\mathcal{T}}_{h}^{bc}$ such that

[TABLE]

Remark 5.1**.**

It is critical that $\tilde{e}$ is properly designed so that the boundary correction does not over-estimate the error due to the underresolved geometry. In our construction, it is easy to see that $\tilde{e}=\tilde{g}-u_{h}$ on $\partial{\Omega}$ where $\tilde{g}$ is the nodal interpolation of $g$ with respect to ${\mathcal{T}}_{h}^{bc}$ on $\partial{\Omega}$. Then we have that

[TABLE]

The second inequality follows from the equivalence of norms.

5.2 Numerical Examples

Example 5.2**.**

The level set of the problem has a flower shape (see Fig. 3) that has the following representation:

[TABLE]

with

[TABLE]

for $i=1,\cdots,8$, $x_{i}=r(\cos(\pi/8)+\sin(\pi/8))\cos(i*\pi/4)$ and $y_{i}=r(\cos(\pi/8)+\sin(\pi/8))\cos(i*\pi/4)$. The domain ${\Omega}$ is defined as $x\in\mathbb{R}^{2}$ such that $\phi(x)\leq 0$. The data are given such that $g=0$ on $\partial{\Omega}$ and

[TABLE]

In the numerical scheme, we take $g_{h}=0$. With the stopping criteria that the total number of degree of freedoms be not greater than $7000$, the final meshes obtained without and with adding the boundary correction term are given $\|\nabla\tilde{e}\|$ in Fig. 4a and Fig. 4b, respectively. We observe that slightly more degree of freedoms are added around the concave corners in Fig. 4b. The corresponding convergence performances of each adaptive procedure are given in Fig. 5. We observe that both estimators converge optimally when the mesh become fine enough. Adding the boundary correction term does slightly increases the estimator at the initial stage. However, the weight is diminishing as the mesh gets finer.

Example 5.3**.**

The level set of this problem has the following representation:

[TABLE]

with $\phi_{0}$ to $\phi_{8}$ defined the same as in Example 5.2. The datum are chosen such that

[TABLE]

In the numerical scheme we approximate the boundary by ${\Omega}_{h}$ defined in Eq. 5.1. For the Dirichlet date we take $g_{h}$ to be the conforming piecewise linear interpolation of $g=y^{2}$. With the stopping criteria that the total number of degree of freedoms be not greater than $7000$, the final meshes obtained without and with adding the boundary correction term $\|\nabla\tilde{e}\|_{K}$ are given in Fig. 6a and Fig. 6b, respectively. The corresponding convergence performance of each adaptive procedure are given in Fig. 7. We observe that in both cases the estimator converges optimally. In Fig. 6b note that there is no more obvious dense refinement comparing to Fig. 6a around the boundary including the corners since in this case the approximation error $g-g_{h}$ is of uniform order everywhere.

5.3 Example 3

In this example we consider the reentrant problem whose solution has the following polar representation:

[TABLE]

with $\alpha=\pi/\omega$ and $\omega$ being the angle of the reentrant corner. In this example, we test two values for $\omega$, i.e., $31\pi/8$ and $63\pi/16$. The stopping criteria is set such that the maximal number of degrees of freedom does not exceed $5000$. In the numerical scheme we approximate the Dirichlet datum $g$ by $g_{h}$ using conforming piecewise linear interpolation. The final meshes generated without adding the boundary correction error are given in Fig. 8a–Fig. 8b and with the boundary correction error are given in Fig. 9a–Fig. 9b. The corresponding convergence rate of estimators are presented in Fig. 10a–Fig. 10b. We again observe the optimal convergence performance for the estimator and true error in all cases. This indicates that the estimator, with or without the boundary correction error, works equivalently effective for problem even with singularity on the boundary.

We now test the same procedure but with $g_{h}$ to be the piecewise constant interpolation of $g$ for $w=31\pi/16$. From the final meshes generated in Fig. 11a and Fig. 11b we observe that more degree of freedoms are added around the boundary in both cases due to the poorer approximation of the boundary data. The estimators, however, still converges optimally in both cases when the meshes are fine enough.

5.4 Example 4

In this example, we consider the problem that has the following representation:

[TABLE]

with $r=\sqrt{x^{2}+y^{2}}$, $\theta=\tan^{-1}(y/x)$, $\alpha=\pi/\omega$ and $\omega=31\pi/16$. This problem has multiple singularities with one reentrant corner at $(0.0)$ and one peak at $(0.5,0.5)$. In the numerical scheme, $g_{h}$ and $f_{h}$ are taken as the linear nodal interpolation of $g$ and $f$, respectively. The stopping criteria is set to not exceed the maximal refinement step of $50$ and the maximal number of degrees of freedom of $7500$. The final meshes are given in Fig. 13a and Fig. 13b and their corresponding convergence performance are presented in Fig. 14.

This example shows that the algorithm works also effectively when multiple singularities occurs, no matter the singularity happens on the boundary or inside the domain.

Acknowledgement

EB and CH were supported by EPSRC, UK, Grant No.

EP/P01576X/1. ML was supported by the Swedish Foundation for Strategic Research Grant No. AM13-0029, the Swedish Research Council Grants No. 2017-03911, and the Swedish strategic research programme eSSENCE.