Finite Element Error Estimates on Geometrically Perturbed Domains

TL;DR

This paper derives error estimates for finite element methods applied to elliptic PDEs on domains with geometric perturbations, highlighting the impact of domain mismatch on overall accuracy.

Contribution

It provides new $H^1$ and $L_2$ error estimates for the Laplace problem on perturbed domains, emphasizing the significance of domain errors.

Findings

Domain mismatch can dominate finite element errors.

Error estimates are validated through computational examples.

Abstract

We develop error estimates for the finite element approximation of elliptic partial differential equations on perturbed domains, i.e. when the computational domain does not match the real geometry. The result shows that the error related to the domain can be a dominating factor in the finite element discretization error. The main result consists of $H^1-$ and $L_2-$ error estimates for the Laplace problem. Theoretical considerations are validated by a computational example.

Figures

Equations

$$-\Delta u=f\text{ in }\Omega ,u=0\text{ on }\partial \Omega .$$

$$(\nabla u,\nabla \varphi {)}_{\Omega }=(f,\varphi {)}_{\Omega }\forall \varphi \in H_0^1(\Omega ).$$

$$\Vert u{\Vert }_{H^{m+2}(\Omega )}\le c\Vert f{\Vert }_{H^m(\Omega )},$$

$${\rm Y}:=\mathrm{dist}(\partial \Omega ,\partial {\Omega }_r):=\max \{\text{adjustlimits}{\sup }_{x\in \partial \Omega }{\inf }_{y\in \partial {\Omega }_r}|x-y|,\text{adjustlimits}{\sup }_{y\in \partial {\Omega }_r}{\inf }_{x\in \partial \Omega }|x-y|\}.$$

$$S:=({\Omega }_r\setminus \Omega )\cup (\Omega \setminus {\Omega }_r).$$

$$\big{|}\bigcup_{i=1}^{n(S)}R_{i}\cap S\big{|}\leq c_{S,2}|S|,$$

$$\Vert f{\Vert }_{H^m({\Omega }_r)}\le c\Vert f{\Vert }_{H^m(\Omega )}.$$

$$(\nabla u_r,\nabla {\varphi }_r{)}_{{\Omega }_r}=(f,{\varphi }_r{)}_{{\Omega }_r}\forall {\varphi }_r\in H_0^1({\Omega }_r),$$

$$\Vert u_r{\Vert }_{H^{m+2}({\Omega }_r)}\le c\Vert f{\Vert }_{H^m({\Omega }_r)}\le c\Vert f{\Vert }_{H^m(\Omega )}.$$

$$\gamma :=\mathrm{dist}(\partial V,\partial W).$$

$$\begin{split}&\|\psi\|_{\partial W\cap V}\leq c\Big{(}\|\psi\|_{\partial V}+\gamma^{\frac{1}{2}}\|\nabla\psi\|_{V\setminus W}\Big{)},\\ &\|\psi\|_{V\setminus W}\leq c\gamma^{\frac{1}{2}}\Big{(}\|\psi\|_{\partial V}+\gamma^{\frac{1}{2}}\|\nabla\psi\|_{V\setminus W}\Big{)},\end{split}$$

$$\big{|}\psi(x_{\partial W})\big{|}^{2}\leq 2\big{|}\psi(x_{\partial V})\big{|}^{2}+2\big{|}\int_{x_{\partial V}}^{x_{\partial W}}\psi^{\prime}(s)\,\text{d}s\big{|}^{2}.$$

$$|\psi (x_{\partial W}){|}^2\le 2|\psi (x_{\partial V}){|}^2+2c_S\gamma {\int }_{x_{\partial V}}^{x_{\partial W}}|\nabla \psi (s){|}^2\text{d}s.$$

$$\int |\psi (g_W^R(t)){|}^2\text{d}t\le \int |\psi (g_V^R(t)){|}^2\text{d}t+2c_S\gamma \int {\int }_{g_V^R(t)}^{g_W^R(t)}|\nabla \psi (s){|}^2\text{d}s\text{d}t.$$

$$\begin{gathered} \frac{1}{{\max }_t\{1+|\nabla g_W^R(t){|}^2\}}{\int }_{(\partial W\cap V)\cap R}|\psi {|}^2\text{d}s \\ \le \frac{1}{{\min }_t\{1+|\nabla g_V^R(t){|}^2\}}{\int }_{R\cap \partial V}|\psi {|}^2\text{d}s+2c_S\gamma {\int }_{R\cap (V\setminus W)}|\nabla \psi {|}^2\text{d}x. \end{gathered}$$

$$\|\psi\|^{2}_{(\partial W\cap V)\cap R}\leq c(\partial V,\partial W)c_{S}\Big{(}\|\psi\|^{2}_{R\cap\partial V}+\gamma\|\nabla\psi\|^{2}_{R\cap S}\Big{)}.$$

$$\|\psi\|^{2}_{\partial W\cap V}\leq c(\partial V,\partial W)c_{S}\Big{(}\|\psi\|^{2}_{\partial V}+\gamma\|\nabla\psi\|^{2}_{V\setminus W}\Big{)}.$$

$$|\psi (x(t,s)){|}^2\le 2|\psi (g_V^R(t)){|}^2+2c_S\gamma \int {\int }_{x(t)}^{g_V^R(t)}|\nabla \psi {|}^2\text{d}s\text{d}t.$$

$$\Vert \psi (x){\Vert }_{R\cap (V\setminus W)}^2\le \frac{2}{{\min }_t\{1+|\nabla g_V^R(t){|}^2\}}\gamma \Vert \psi {\Vert }_{R\cap \partial V}^2+2c_S{\gamma }^2\Vert \nabla \psi {\Vert }_{R\cap (V\setminus W)}^2.$$

$$\Vert u-u_r{\Vert }_{\Omega }+{{\rm Y}}^{\frac{1}{2}}\Vert \nabla (u-u_r){\Vert }_{\Omega }\le c{\rm Y}\Vert f{\Vert }_{\Omega \cup {\Omega }_r}.$$

$$\begin{split}\|\nabla(u-u_{r})\|_{\Omega}^{2}=&\big{(}\nabla(u-u_{r}),\nabla(u-u_{r})\big{)}_{\Omega\cap\Omega_{r}}+\big{(}\nabla(u-u_{r}),\nabla(u-u_{r})\big{)}_{\Omega\setminus\Omega_{r}}\\ =&-\big{(}\Delta(u-u_{r}),u-u_{r}\big{)}_{\Omega}+\langle\partial_{n}(u-u_{r}),u-u_{r}\rangle_{\partial(\Omega\cap\Omega_{r})}\\ &\quad+\langle\partial_{n}(u-u_{r}),u-u_{r}\rangle_{\partial(\Omega\setminus\Omega_{r})}.\end{split}$$

$$\begin{split}\|\nabla(u-u_{r})\|_{\Omega}^{2}=&-\big{(}\Delta(u-u_{r}),u-u_{r}\big{)}_{\Omega}+\langle\partial_{n}(u-u_{r}),u-u_{r}\rangle_{\partial\Omega}\\ &\quad+\langle[\partial_{n}(u-u_{r})],u-u_{r}\rangle_{\partial\Omega_{r}\cap\Omega}.\end{split}$$

$$-(\Delta (u-u_r),u-u_r{)}_{\Omega }=-(\Delta (u-u_r),u-u_r{)}_{\Omega \cap {\Omega }_r}-(\Delta u,u{)}_{\Omega \setminus {\Omega }_r}=(f,u{)}_{\Omega \setminus {\Omega }_r}.$$

$$\begin{gathered} ⟨{\partial }_n(u-u_r),u-u_r{⟩}_{\partial \Omega }+⟨[{\partial }_n(u-u_r)],u-u_r{⟩}_{\partial {\Omega }_r\cap \Omega } \\ =-⟨{\partial }_n(u-u_r),u_r{⟩}_{\partial \Omega \cap {\Omega }_r}-⟨{\partial }_n{u}_r,u{⟩}_{\partial {\Omega }_r\cap \Omega }. \end{gathered}$$

$$\begin{gathered} \Vert \nabla (u-u_r){\Vert }_{\Omega }^2\le \Vert f{\Vert }_{\Omega \setminus {\Omega }_r}\Vert u{\Vert }_{\Omega \setminus {\Omega }_r} \\ +\Vert {\partial }_n(u-u_r){\Vert }_{\partial \Omega \cap {\Omega }_r}\Vert u_r{\Vert }_{\partial \Omega \cap {\Omega }_r}+\Vert {\partial }_{n}u{\Vert }_{\partial {\Omega }_r\cap \Omega }\Vert u{\Vert }_{\partial {\Omega }_r\cap \Omega }. \end{gathered}$$

$$\|\nabla(u-u_{r})\|_{\Omega}^{2}\leq\|f\|_{\Omega\setminus\Omega_{r}}\|u\|_{\Omega\setminus\Omega_{r}}\\ +c\big{(}\|u\|_{H^{2}(\Omega)}+\|u_{r}\|_{H^{2}(\Omega_{r})}\big{)}\big{(}\|u_{r}\|_{\partial\Omega\cap\Omega_{r}}+\|u\|_{\partial\Omega_{r}\cap\Omega}\big{)}.$$

$$\Vert u{\Vert }_{\partial {\Omega }_r\cap \Omega }$$

$$\Vert u_r{\Vert }_{\partial \Omega \cap {\Omega }_r}$$

$$\Vert u{\Vert }_{\partial {\Omega }_r\cap \Omega }\le c{\rm Y}\Vert f{\Vert }_{\Omega },\Vert u_r{\Vert }_{\partial \Omega \cap {\Omega }_r}\le c{\rm Y}\Vert f{\Vert }_{\Omega }.$$

$$\|u\|_{\Omega\setminus\Omega_{r}}\leq c\Upsilon\|\nabla u\|_{\Omega\setminus\Omega_{r}}\leq c\Upsilon^{\frac{3}{2}}\Big{(}\|u\|_{H^{2}(\Omega)}+\Upsilon^{\frac{1}{2}}\|u\|_{H^{2}(\Omega)}\Big{)}\leq c\Upsilon^{\frac{3}{2}}\|f\|_{\Omega}.$$

$$\|\nabla(u-u_{r})\|^{2}_{\Omega}\leq c\Big{(}\Upsilon^{\frac{3}{2}}+\Upsilon\Big{)}\|f\|^{2}_{\Omega},$$

Full text

1 Introduction

The main aim of this work is to develop finite element (FE) error estimates in the case when there is uncertainty with respect to the computational domain. We consider the question of how a domain related error affects the finite element discretization error. We use the conforming finite element method (FEM) which is well established in the scientific computing community and allows for a rigorous analysis of the approximation error ErnGuermond2004 .

Our motivation is as follows. The steps to obtain a mesh for FE computations often come with some uncertainty, for example related to empirical measurements or image processing techniques, e.g. medical image segmentation Oberkampf2006 ; OberkampfBook . Therefore, we often perform computations on a domain which is an approximation of the real geometry, i.e., the computational domain is close to but does not match the real domain. In this work we do not specify the source of the error, but we take the error into account by explicitly using the error laden reconstructed domains.

This theoretical result is of great importance for scientific computations. Vast numbers of engineering branches rely on the results of computational fluid dynamics simulations, where there is often uncertainty connected to the computational domain. A prime example of this is computational based medical diagnostics, where shapes are reconstructed from inverse problems, such as computer tomography. The assessment of error attributed to the limited spatial resolution of magnetic resonance techniques has been discussed in Moore1998 ; Moore1999 . For a survey on computational vascular fluid dynamics, where modeling and reconstruction related issues are discussed, we refer to Quarteroni2000 . Error analysis of computational models is a key factor for assessing the reliability for virtual predictions.

Uncertainties in the computational domain have been studied from the numerical perspective. Rigorous bounds for elliptic problems on random domains have been derived, for approximate problems defined on a sequence of domains that is supposed to converge in the set sense to a limit domain, for both Dirichlet Babuska2003 and Neumann Babuska2002 boundary conditions. Although our techniques are similar, we consider a case where the geometrical error is not small, but where it might dominate the discretization error.

When measurement data is available the accuracy of numerical predictions can be improved by data assimilation techniques. Applications of variational data assimilation in computational hemodynamics have been reviewed in DElia2012 . For recent developments we refer to Funke2019 and Nolte . On the other hand, the treatment of boundary uncertainty can be cast into a probabilistic framework. The domain mapping method is based entirely on stochastic mappings to transform the original deterministic/stochastic problem in a random domain into a stochastic problem in a deterministic domain, see Xiu2006 ; Tartakovsky2006 ; Harbrecht2016 . The perturbation method starts with a prescribed perturbation field at the boundary of a reference configuration and uses a shape Taylor expansion with respect to this perturbation field to represent the solution Harbrecht2008 . In Allaire2015 and Dambrine2015 a similar technique was used to incorporate random perturbations of a given domain in the context of shape optimization. Moreover, the fictitious domain approach and a polynomial chaos expansion have been applied in Canuto2007 . We note, that the probabilistic approach is beyond the scope of this work and the introduction of the boundary uncertainty as random variable increases the complexity of the problem.

The above approaches incorporate additional information on the domain reconstruction, such as measurement data or a probabilistic distribution of the approximation error. In comparison to these approaches our result can be seen as the worst case scenario. We only require that the distance between the two domains is bounded.

The analysis presented in this paper starts with well-known results regarding the finite element approximation on domains with curved boundaries. But in contrast to these estimates we cannot expect the error coming from the approximation of the geometry is small or even converging to zero. Instead we split the error into a geometric approximation error between real domain and perturbed domain and into an error coming from the finite element discretization of the problem on the perturbed domain. A central step is Lemma 4 which estimates the geometry perturbation. Having in mind that this error is not small and cannot be reduced by means of tuning the discretization, the typical application case is to balance both error contributions to efficiently reach the barrier of the geometry error. Theorem 3.2 gives such optimally balanced estimates that include both error contributions.

This paper is organized as follows. After this introduction, in Section 2 we introduce the mathematical setting and some required auxiliary results. Section 3 covers finite element discretization and proves the main results of this work. We illustrate our result with computational examples in Section 4.

2 Mathematical setting and auxiliary result

2.1 Notation

Let $\Omega\subset\mathds{R}^{d}$ be a domain with dimension $d\in\{2,3\}$. By $L^{2}(\Omega)$ we denote the Lebesgue space of square integrable functions equipped with the norm $\|\cdot\|_{\Omega}$. By $H^{1}(\Omega)$ we denote the space of $L^{2}(\Omega)$ functions with first weak derivative in $L^{2}(\Omega)$ and by $H^{m}(\Omega)$ for $m\in\mathds{N}_{0}$ we denote the corresponding generalizations with weak derivatives up to degree $m\in\mathds{N}_{0}$. The norms in $H^{m}(\Omega)$ are denoted by $\|\cdot\|_{H^{m}(\Omega)}$. For convenience we use the notation $H^{0}(\Omega)\coloneqq L^{2}(\Omega)$. By $H^{1}_{0}(\Omega)$ we denote the space of those $H^{1}(\Omega)$ functions that have vanishing trace on the domain’s boundary $\partial\Omega$ and we use the notation $H^{1}_{0}(\Omega;\Gamma)$ if the trace only vanishes on a part of the boundary, $\Gamma\subset\partial\Omega$. Further, by $(\cdot,\cdot)_{\Omega}$ we denote the $L^{2}(\Omega)$-scalar product and $\langle\cdot,\cdot\rangle_{\Gamma}$ the $L^{2}$-scalar product on a $d-1$ dimensional manifold $\Gamma$, e.g. $\Gamma=\partial\Omega$. Moreover, $[\partial_{n}\psi]$ is the jump of the normal derivative of $\psi$, i.e. for $x\in\Gamma$ with normal $\vec{n}$ (that is normal to $\Gamma$) $[\partial_{n}\psi](x):=\lim_{h\searrow 0}\partial_{n}\psi(x+h\vec{n})-\lim_{h\searrow 0}\partial_{n}\psi(x-h\vec{n}).$

2.2 Laplace equation and domain perturbation

On $\Omega\subset\mathds{R}^{d}$ let $f\in L^{2}(\Omega)$ be the given right hand side. We consider the Laplace problem with homogeneous Dirichlet boundary conditions,

[TABLE]

The variational formulation of this problem is given by: find $u\in H^{1}_{0}(\Omega)$, such that

[TABLE]

The boundary $\partial\Omega$ is supposed to have a parametrization in $C^{m+2}$, where $m\in\mathds{N}_{0}$. Given the additional regularity $f\in H^{m}(\Omega)$, $H^{0}(\Omega)\coloneqq L^{2}(\Omega)$, there exists a unique solution satisfying the a-priori estimate

[TABLE]

see e.g. evans .

In the following we assume that the real domain $\Omega$ is not exactly known but only given up to an uncertainty. We hence define a second domain, the reconstructed domain $\Omega_{r}$ with a boundary that allows for $C^{m+2}$ parametrization. The Hausdorff distance between both domains is then denoted by $\Upsilon\in\mathds{R}$,

[TABLE]

This distance $\Upsilon$ is not necessarily small. When it comes to spatial discretization we will be interested in both cases, $h\ll\Upsilon$ as well as $\Upsilon\ll h$, where $h>0$ is the mesh size. The two domains do not match and either domain can protrude from the other, see Figure 1. In order to prove our error estimates we require the following technical assumption on the relation between the two domains $\Omega$ and $\Omega_{r}$.

Assumption 1** (Domains)**

Let $\Omega$ and $\Omega_{r}$ be two domains with $\Omega\cap\Omega_{r}\neq\emptyset$ and with Hausdorff distance $\Upsilon\in\mathds{R}$. Both boundaries allow for a local $C^{m+2}$ parametrization, $m\in\mathds{N}_{0}$. Let

[TABLE]

We assume that there exists a cover of $S$ by a finite number of open rectangles (or rectangular cuboids) $\{R_{1},\dots,R_{n(S)}\}$. Each rectangle $R_{i}$ is given as translation and rotation of $(0,h_{i})\times(0,t_{i})$ for $d=2$, or $(0,h_{i})\times(0,t_{i}^{1})\times(0,t_{i}^{2})$ for $d=3$, where the height $h_{i}$ is bounded by $h_{i}\leq c_{S,1}\,\Upsilon$ with a constant $c_{S,1}\geq 0$. Following conditions hold:

1. 1.a)

On each rectangle $R$, the boundary lines $\partial\Omega\cap R$ and $\partial\Omega_{r}\cap R$ allow for unique parametrizations $g^{R}_{\Omega}(t)$ and $g^{R}_{\Omega_{r}}(t)$ over the base $t$, or $(t^{1},t^{2})$ for $d=3$, respectively. 2. 1.b)

The area of the cover is bounded by the area of the remainder $S$, i.e.

[TABLE]

where $c_{S,2}>0$ is a constant.

For the following we set $c_{S}:=\max\{c_{S,1},c_{S,2}\}$.

Figure 1 shows such a cover for different domain remainders. From Assumption A1 we deduce that each line through the height of the rectangle (marked in red in the figure) cuts each of the two boundaries exactly one time. The second assumption limits the overlap of the rectangles. These are shown as the shaded in the left sketch in Figure 1. Both assumptions on the domain are required for the proof of Lemma 3 that is based on Fubini’s integral theorem. A more flexible framework that allows for a wider variety of domains, e.g. with boundaries that feature hooks, could be based on the construction of a map between two boundary segments on $\partial\Omega_{r}$ and $\partial\Omega$. Such approaches play an important role in isogeometric analysis. We refer to Xu2011 and Xu2013 for examples on the construction of such maps.

To formulate the Laplace equation on the reconstructed domain $\Omega_{r}$ we must face the technical difficulty that the right hand side $f\in H^{m}(\Omega)$ is not necessarily defined on $\Omega_{r}$. We therefore weaken the assumptions on the right hand side.

Assumption 2** (Right hand side)**

Let $f\in H^{m}_{\text{loc}}(\mathds{R}^{d})$, i.e. $f\in H^{m}(G)$ for each compact subset $G\subset\mathds{R}^{d}$. In addition we assume that the right hand side on $\Omega_{r}$ can be bounded by the right hand side on $\Omega$, i.e.

[TABLE]

An alternative would be to use Sobolev extension theorems to extend functions $f\in H^{m}(\Omega)$ from $\Omega$ to $\Omega_{r}$, see Calderon1961 .

On $\Omega_{r}$ we define the solution $u_{r}\in H^{1}_{0}(\Omega_{r})$ to the perturbed Laplace problem

[TABLE]

The unique solution to (5) satisfies the bound

[TABLE]

Remark 1 (Extension of the solutions)

A difficulty for deriving error estimates is that $u$ is defined on $\Omega$ and $u_{r}$ on $\Omega_{r}\neq\Omega$. Since the domains do not match, $u$ may not be defined on all of $\Omega_{r}$ and vice versa. To give the expression $u-u_{r}$ a meaning on all domains we extend both solutions by zero outside their defining domains, i.e. $u:=0$ on $\Omega_{r}\setminus\Omega$ and $u_{r}:=0$ on $\Omega\setminus\Omega_{r}$. Globally, both functions still have the regularity $u,u_{r}\in H^{1}(\Omega\cup\Omega_{r})$. We will use the same notation for discrete functions $u_{h}\in V_{h}$ defined on a mesh $\Omega_{h}$ and extend them by zero to $\mathds{R}^{d}$.  $\triangle$

The following preliminary results are necessary in the proof of the main estimates. They can be considered as variants of the trace inequality and of Poincaré’s estimate, respectively.

Lemma 3

Let $\gamma\in\mathds{R},\gamma>0$, $V\subset\mathds{R}^{d}$ and $W\subset\mathds{R}^{d}$ for $d\in\{2,3\}$ be two domains with boundaries $\partial V$ and $\partial W$ that satisfy Assumption 1 with distance

[TABLE]

For $\psi\in C^{1}(V)\cap C(\bar{V})$ it holds

[TABLE]

where the constants $c>0$ depend on $c_{S}$ from Assumption 1 and the curvature of the domain boundaries.

Proof

Let $R$ be one rectangle of the cover and let $x_{\partial W}=g^{R}_{W}(t)\in\partial(W\cap V)\cap R$, see Figure 2. By $x_{\partial V}=g^{R}_{V}(t)\in\partial(W\cap V)\cap R$ we denote the corresponding unique point on $\partial W\cap R$. The connecting line segment $\overline{x_{\partial V}x_{\partial W}}$ completely runs through $V\setminus W$, as, if the line would leave this remainder, it would cut each line more than once which opposes Assumption 1.b). Integrating the function $\psi$ along this line gives

[TABLE]

Applying Hölder’s inequality to the second term on the right hand side, with the length of the line bounded by $c_{S}\gamma$, we obtain

[TABLE]

Using the parametrizations $x_{\partial W}=g^{R}_{W}(t)$ and $x_{\partial V}=g^{R}_{V}(t)$ we integrate over $t$ which gives

[TABLE]

The volume integral on the right hand side is exactly the integral over $R\cap(V\setminus W)$. The boundary integrals can be interpreted as path integrals and therefore be estimated by

[TABLE]

As the boundaries allow for a $C^{2}$ parametrization, we estimate

[TABLE]

Summation over all rectangles and estimation of all overlaps by means of Assumption 1 gives

[TABLE]

For $\psi\in H^{1}_{0}(V)$, the term on $\partial V$ vanishes.

To show the second estimate on $V\setminus W$ we again pick one rectangle $R$ and consider a point $x\in V\setminus W$ on the line connecting $x_{\partial V}=g^{R}_{V}(t)$ and $x_{\partial W}=g^{R}_{W}(t)$ such that we introduce the notation $x=x(t,s)$. By the same arguments as above it holds

[TABLE]

We integrate over $s$ and $t$ to obtain

[TABLE]

Summing over all rectangles gives the desired result.

The above lemma is later used in such a way that $V$ and $W$ can be substituted as both $\Omega$ and $\Omega_{r}$, specifically to the case of use.

We continue by estimating the difference between the solutions of the Laplace equations on $\Omega$ and on $\Omega_{r}$.

Lemma 4

Let $\Omega,\Omega_{r}\in\mathds{R}^{d}$ with $\partial\Omega,\partial\Omega_{r}\in C^{m+2}$ satisfying $\operatorname{dist}(\partial\Omega,\partial\Omega_{r})<\Upsilon$ as well as Assumption 1. Furthermore, let $f\in L^{2}_{\text{loc}}(\mathds{R}^{d})$ satisfy Assumption 2 and let $f_{r}\coloneqq f|_{\Omega_{r}}$. For the solutions $u\in H^{1}_{0}(\Omega)\cap H^{2}(\Omega)$ and $u_{r}\in H^{1}_{0}(\Omega_{r})\cap H^{2}(\Omega_{r})$ to (2) and (5) respectively, it holds that

[TABLE]

Proof

(i) We continuously extend $u$ and $u_{r}$ by zero to $\mathds{R}^{d}$, c.f. Remark 1, such that $u-u_{r}\in H^{1}(\Omega\cup\Omega_{r})$ is well defined. We separate the domains of integration and integrate by parts

[TABLE]

Combining the boundary terms on the right-hand side of (11), into an integral over $\partial\Omega$ and a jump term over $\partial\Omega_{r}\cap\Omega$, we obtain

[TABLE]

In $\Omega\cap\Omega_{r}$ it holds $f=f_{r}$ and hence (weakly) $-\Delta(u-u_{r})=0$, such that

[TABLE]

On $\partial\Omega$ it holds $u=0$ and on $\partial\Omega_{r}\cap\Omega$ it holds $u_{r}=0$. Further, since $u\in H^{2}(\Omega)$ it holds that $[\partial_{n}u]=0$ on $\partial\Omega_{r}\cap\Omega$. Finally, $u_{r}=0$ on $\Omega\setminus\Omega_{r}$, such that the boundary terms reduce to

[TABLE]

Combining (12)-(14) and using the Cauchy-Schwarz inequality, we estimate

[TABLE]

Since $u,u_{r}\in H^{2}(\Omega\cap\Omega_{r})$, the trace inequality gives

[TABLE]

Applying Lemma 3 twice: to $\psi=u$ and to $\psi=\nabla u$ (same for $u_{r})$, and extending the norms from $\Omega\setminus\Omega_{r}$ to $\Omega$ and from $\Omega_{r}\setminus\Omega$ to $\Omega_{r}$ give the bounds

[TABLE]

With the trace inequality and the a priori estimates $\|u\|_{H^{2}(\Omega)}\leq c\|f\|_{\Omega}$ and $\|u_{r}\|_{H^{2}(\Omega_{r})}\leq c\|f_{r}\|_{\Omega_{r}}\leq c\|f\|_{\Omega}$ we obtain the bounds

[TABLE]

Using the fact that $u=0$ on $\partial\Omega$ we apply (7) twice and use the trace inequality to get the estimate

[TABLE]

We can then estimate $\|f\|_{\Omega\setminus\Omega_{r}}\leq\|f\|_{\Omega}$ by extending to the complete domain. Combining (16) with (18) and (19) we obtain the estimate

[TABLE]

which concludes the $H^{1}$-norm bound.

(ii) For the $L^{2}$-estimate we introduce the adjoint problem

[TABLE]

which allows for a unique solution satisfying the a-priori bound $\|z\|_{H^{2}(\Omega)}\leq c_{s}$ with the stability constant $c_{s}<\infty$. Testing with $u-u_{r}$ and integrating by parts twice gives

[TABLE]

It holds $z=0$ and $u=0$ on $\partial\Omega$, $[\partial_{n}u]=0$ on $\partial\Omega_{r}\cap\Omega$ and $-\Delta(u-u_{r})=0$ in $\Omega\cap\Omega_{r}$ such that we get

[TABLE]

The boundary terms $\|z\|_{\partial\Omega_{r}\cap\Omega}$ and $\|u_{r}\|_{\partial\Omega}$ are estimated with Lemma 3, the normal derivatives by the trace inequality and the terms on $\Omega\setminus\Omega_{r}$ by (7)

[TABLE]

The $L^{2}$-norm estimate follows by using the bounds $\|u\|_{H^{2}(\Omega)}\leq c\|f\|_{\Omega}$, $\|u_{r}\|_{H^{2}(\Omega_{r})}\leq c\|f\|_{\Omega}$ and $\|z\|_{H^{2}(\Omega)}\leq c$.

Remark 2

The estimate $\|f\|_{\Omega\setminus\Omega_{r}}\leq c\|f\|_{\Omega}$ is not optimal. Further powers of $\Upsilon$ are easily generated at the cost of a higher right hand side regularity. Also, the estimate $\|\partial_{n}(u-u_{r})\|\leq c(\|u\|_{H^{2}(\Omega)}+\|u_{r}\|_{H^{2}(\Omega_{r})})$ by Cauchy Schwarz and the trace inequality could be enhanced to produce powers of $\Upsilon$. The limiting term in (12) however is the boundary integral $|\langle\partial_{n}u_{r},u\rangle_{\partial\Omega_{r}\cap\Omega}|=\mathcal{O}(\Upsilon^{\frac{1}{2}})$ which is optimal in the $H^{1}$-estimate. In Remark 4 and Corollary 1 we present an estimate that focuses on the intersection $\Omega\cap\Omega_{r}$ only and that allows us to improve the order to $\mathcal{O}(\Upsilon)$ in the $H^{1}$-case by avoiding exactly this boundary integral.  $\triangle$

3 Discretization

The starting point of a finite element discretization is the mesh of the domain $\Omega$. In our setting we do not mesh $\Omega$ directly, because the domain $\Omega$ is not exactly known. Instead, we consider a mesh of the reconstructed domain $\Omega_{r}$.

We partition $\Omega_{r}$ into a parametric triangulation $\Omega_{h}$, consisting of open elements $T\subset\mathds{R}^{d}$. Each element $T\in\Omega_{h}$ stems from a unique reference element $\hat{T}$ which is a simple geometric structure such as a triangle, quadrilateral or tetrahedron. The numerical examples in Section 4 are based on quadrilateral meshes. The map $T_{T}:\hat{T}\to T$ is a polynomial of degree $r\in\mathds{N}$. We will consider iso-parametric finite element spaces, that are based on polynomials of the same degree $r\geq 1$. In the following we assume structural and shape regularity of the mesh such that standard interpolation estimates

[TABLE]

will hold for all elements $T\in\Omega_{h}$, c.f. Bernardi1989 . The discretization parameter $h$ represents the size of the largest element in the mesh. See (Richter2017, , Section 4.2.2) for a detailed description.

On the reference element $\hat{T}$ let $\hat{P}$ be a polynomial space of degree $r$, e.g.

[TABLE]

on quadrilateral and hexahedral meshes. Then, the finite element space $V_{h}^{r}$ on the mesh $\Omega_{h}$ is defined as

[TABLE]

This parametric finite element space does not exactly match the domain $\Omega_{r}$. Given an iso-parametric mapping of degree $r$ it holds $\operatorname{dist}(\partial\Omega_{r},\partial\Omega_{h})=\mathcal{O}(h^{r+1})$ and finite element approximation error and geometry approximation error are balanced. Iso-parametric finite elements for the approximation on domains with curved boundaries are well established ErgatoudisIronsZienkiewicz1968 , optimal interpolation and finite element error estimates have been presented in (Ciarlet2002, , Section 4.4). The case of higher order elements with optimal order energy norm estimates is covered in Lenoir1986 . From (Richter2017, , Theorem 4.37) we cite the following approximation result for the iso-parametric approximation of the Laplace equation that also covers the $L^{2}$-error and which is formulated in a similar notation.

Theorem 3.1

Let $m\in\mathds{N}_{0}$ and let $\Omega_{r}$ be a domain with a boundary that allows for a parametrization of degree $m+2$. Let $f_{r}\in H^{m}(\Omega_{r})$ and $u_{h}\in V_{h}^{r}\cap H^{1}_{0}(\Omega_{h})$ be the iso-parametric finite element discretization of degree $1\leq r\leq m+1$

[TABLE]

It holds

[TABLE]

We formulated the error estimate on the domain $\Omega_{r}$ although the finite element functions are given on $\Omega_{h}$ only. To give Theorem 3.1 meaning, we consider all functions extended by zero as described in Remark 1. Combining these preliminary results directly yields the a priori error estimates.

Theorem 3.2

Let $m\in\mathds{N}_{0}$, $\Omega$ and $\Omega_{r}$ be domains with $C^{m+2}$ boundary, distance $\Upsilon$ and that satisfy Assumption 1. Let $\Omega_{h}$ be the iso-parametric mesh of $\Omega_{r}$ with degree $1\leq r\leq m+1$ and let $f\in H^{r-1}_{loc}(\mathds{R}^{d})$ satisfy Assumption 2. For the finite element error between the fully discrete solution $u_{h}\in V_{h}^{r}$

[TABLE]

and the true solution $u\in H^{1}_{0}(\Omega)\cap H^{m+2}(\Omega)$ it holds

[TABLE]

as well as

[TABLE]

Proof

(i) We start with the $H^{1}$ error. Inserting $\pm u_{r}$ and extending the finite element error $u_{r}-u_{h}$ from $\Omega$ to $\Omega_{r}$, where a small remainder appears, we have

[TABLE]

The first and the second term on the right hand side are estimated by Lemma 4 and Theorem 3.1 and, since $u_{r}=0$ on $\Omega\setminus\Omega_{r}$, we obtain

[TABLE]

We continue with the remainder $\nabla u_{h}$ on $\Omega\setminus\Omega_{r}$, which is non-zero on $\Omega_{h}$ only

[TABLE]

This remaining stripe has the width

[TABLE]

and we apply Lemma 3 to get

[TABLE]

where the second derivative $\nabla^{2}u_{h}$ is understood element wise. This term is extended to $\Omega_{h}$ and with the inverse estimate and the a priori estimate for the discrete solution we obtain with $\gamma_{h,\Upsilon}^{2}={\cal O}(h^{2r+2})$ that

[TABLE]

To the first term on the right hand side of (23) we add $\pm u_{r}$ and $\pm I_{h}u_{r}$, the nodal interpolation of $u_{r}$ into the finite element space

[TABLE]

Here, the first and last terms are estimated with the trace inequalities and, in the case of the discrete term with the inverse inequality111We refer to (DiPietro2012, , Chapter 1.4.3) or Barrenechea2017 ; Cangiani2018 for recent developments on the local trace inequality and the inverse estimate on meshes with curved boundaries., followed by adding $\pm u_{r}$ we get

[TABLE]

We used both $\gamma_{h,\Upsilon}={\cal O}(h^{r+1})$ and $\gamma_{h,\Upsilon}={\cal O}(\Upsilon)$. Then, collecting all terms in (22)-(26) and using the interpolation estimates as well as Theorem 3.1 we finally get

[TABLE]

which shows the a priori estimate since $3r-1\leq 2r$ for all $r\geq 1$.

(ii) For the $L^{2}$-error we proceed in the same way, but the remainder appearing in (21) does not carry any derivative, such that, instead of (23) the optimal order variant of Lemma 3 with integration to the boundary $\partial\Omega_{h}$, where $u_{h}=0$, can be applied, i.e.

[TABLE]

The $L^{2}$-estimate directly follows with Lemma 4, Theorem 3.1 and by the a priori estimate $\|\nabla u_{h}\|^{2}_{\Omega\setminus\Omega_{r}}\leq c\|\nabla u_{h}\|^{2}_{\Omega_{h}}\leq c\|f\|^{2}_{\Omega_{r}}$.

Remark 3 (Polygonal domains)

In two dimensions, the extension of the error estimates to the case of convex polygonal domains, where $u\in H^{2}(\Omega)$ and $u_{r}\in H^{2}(\Omega_{r})$, is relatively straightforward. In this case, $\Omega_{h}$ fits $\Omega_{r}$ such that the finite element error $u_{r}-u_{h}$ can be estimated with the standard a priori result $\|u_{r}-u_{h}\|+h\|\nabla(u_{r}-u_{h})\|\leq c\|f\|$. The extension of Lemma 3, which locally requires smoothness of the parametrizations $g^{R}_{W}(\cdot)$ and $g^{R}_{V}(\cdot)$, see steps (8)-(10), can be accomplished by refining the cover of the domain which is described in Assumption 1, see also Figure 1: All rectangles are split in such a way that the corners of $\partial\Omega$ and $\Omega_{r}$ are cut by the edges of rectangles. This allows to derive the optimal error estimates $\|u-u_{h}\|_{H^{1}(\Omega)}={\cal O}(\Upsilon^{\frac{1}{2}}+h)$ and $\|u-u_{h}\|_{H^{1}(\Omega)}={\cal O}(\Upsilon+h^{2})$. In three dimensions, such a simple refinement of the cover is not possible and the extension to polygonal domains is more involved.  $\triangle$

Remark 4 (Optimality of the estimates)

Two ingredients govern the error estimates:

A geometrical error of order $\mathcal{O}(\Upsilon^{\frac{1}{2}})$ and $\mathcal{O}(\Upsilon)$, that describes the discrepancy between $\Omega$ and $\Omega_{r}$, in the $H^{1}$ and $L^{2}$ norms respectively. This term is optimal which is easily understood by considering a simple example illustrated in Figure 3, namely $-\Delta u=4$ on the unit disc $\Omega=B_{1}(0)$ and $-\Delta u_{r}=4$ on the shifted domain $\Omega_{r}=B_{1}(\Upsilon)$. The errors in $H^{1}$ norm and $L^{2}$ norms expressed on the complete domain $\Omega$ are estimated by

[TABLE]

A closer analysis shows that the main error – in the $H^{1}$-case – occurs on the small shaded stripe $\Omega\setminus\Omega_{r}$ such that

[TABLE]

while the $L^{2}$-error in $\Omega\cap\Omega_{r}$ is optimal

[TABLE] 2. 2.

The usual Galerkin error $\|u_{r}-u_{h}\|_{\Omega_{r}}+h\|\nabla(u-u_{r})\|_{\Omega_{r}}=\mathcal{O}(h^{r+1})$ of iso-parametric finite element approximations contributes to the overall error. For $\Omega=\Omega_{r}$, i.e. $\Upsilon=0$, this would be the complete error. This estimate is optimal, as it shows the same order as usual finite element bounds on meshes that resolve the geometry.

$\triangle$

In Section 4 we discuss the difficulty of measuring errors on an unknown domain $\Omega$. The optimality of the error estimates is difficult to verify which is mainly due to the technical problems in evaluating norms on the domain remainders $\Omega\setminus\Omega_{r}$, where no finite element mesh is given. These remainders contribute the lowest order parts $\Upsilon^{\frac{1}{2}}$ in the overall error. The following corollary is closer to the setting of the numerical examples and it yields the approximation of order $\Upsilon$ in the $H^{1}$-norm error. In addition to the previous setting we require a regular map $T_{r}:\Omega\to\Omega_{r}$ between the two domains. By pulling back $\Omega_{r}$ to $\Omega$ via this map a Jacobian arises that controls the geometrical error and that hence has to be controllable by $\Upsilon$.

Corollary 1

In addition to the assumptions of Theorem 3.2 let there be a $C^{1}$-diffeomorphism

[TABLE]

satisfying

[TABLE]

Further, let the following regularity of problem data hold in addition to Assumption 2

[TABLE]

and let the solution satisfy

[TABLE]

Then, it holds

[TABLE]

Proof

We start by splitting the error into domain approximation and finite element approximation errors

[TABLE]

An optimal order estimate of the finite element error

[TABLE]

is given in Theorem 3.1. To estimate the first term of the right hand side of (31) we introduce the function

[TABLE]

which satisfies $\hat{u}_{r}\in H^{1}_{0}(\Omega)$ and solves the problem

[TABLE]

where $\hat{f}_{r}(x):=f(T_{r}(x))$ and where $F_{r}:=\nabla T_{r}$ and $J_{r}:=\operatorname{det}(F_{r})$. See (Richter2017, , Section 2.1.2) for details of this transformation of the variational formulation. To estimate the domain approximation error in (31) we introduce $\pm\hat{u}_{r}$ to obtain

[TABLE]

We introduce the notation $e_{r}:=u-\hat{u}_{r}$, extend the first term from ${\Omega\cap\Omega_{r}}$ to $\Omega$ and insert $\pm J_{r}F_{r}^{-1}F_{r}^{-T}\nabla\hat{u}_{r}$ which gives

[TABLE]

where we also used Poincaré’s estimate. For bounding $f-\hat{f}_{r}$ we consider a point $x\in\Omega\cap\Omega_{r}$, use the higher regularity of the right hand side (29) to estimate by a Taylor expansion

[TABLE]

where $\xi\in\Omega$ is some point on the line from $x$ to $T_{r}(x)$. We take the square and integrate over $\Omega$ to get the estimate

[TABLE]

where $\Omega_{\Upsilon}$ is a enlargement of $\Omega$ by at most ${\cal O}(\Upsilon)$, since intermediate values $\xi$ used in (35) are not necessarily part of $\Omega\cup\Omega_{r}$. This argument is also applicable to the second term on the right hand side of (33) such that it holds

[TABLE]

Combining this with (31), (32), (33) and (34) finishes the proof.

Unfortunately this corollary can not be applied universally as the existence of a suitable map $T_{r}:\Omega\to\Omega_{r}$ depends on the given application. Here a construction, corresponding to the ALE map, can be realised by means of a domain deformation $\hat{d}:\Omega\to\mathds{R}^{2}$

[TABLE]

Such a construction is common in fluid-structure interactions, see (Richter2017, , Section 2.5.2). Given that $|\hat{d}|,|\nabla\hat{d}|=\mathcal{O}(\Upsilon)$ it holds

[TABLE]

While the assumption $|\hat{d}|=\mathcal{O}(\Upsilon)$ is easy to satisfy since $\operatorname{dist}(\partial\Omega,\partial\Omega_{r})\leq\Upsilon$, the condition $|\nabla\hat{d}|=\mathcal{O}(\Upsilon)$ will strongly depend on the shape and regularity of the boundary.

We conclude by discussing a simple application of this corollary. Figure 4 illustrates the setting. Let $\Omega$ be the unit sphere, $\Omega_{r}$ be an ellipse

[TABLE]

It holds $\operatorname{dist}(\partial\Omega,\partial\Omega_{r})\leq\Upsilon$ and we define the map $T_{r}:\Omega\to\Omega_{r}$ by

[TABLE]

This map satisfies the assumptions of the corollary

[TABLE]

4 Numerical illustration

In this section we illustrate the theoretical results from the previous section. We compute the Laplace problem on a family of domains representing different values of $\Upsilon$. Moreover, we numerically extend the analytical predictions and show that a similar behavior holds for the Stokes system.

We consider $\Omega$ to be a unit ball in two and three dimensions and define a family of perturbed domains $\Omega_{\Upsilon}$, with the amplitude of the perturbation being dependent on the coefficient $\Upsilon$, cf. Figure 5.

In two dimensions, the boundary of the domain $\Omega_{\Upsilon}$ is given in polar coordinates $(\rho,\varphi)$ by

[TABLE]

and in three dimensions in spherical coordinates $(\rho,\theta,\varphi)$ by

[TABLE]

For computations we take

[TABLE]

In order to illustrate the convergence result from Theorem 3.2, we compute the model problem on a series of uniformly refined meshes. The dependence between the mesh size $h$ and the refinement level $L$ reads $h=2^{-L}$. We denote the mesh approximating $\Omega_{\Upsilon}$, with a mesh size $h$, by $\Omega_{h,\Upsilon}$.

The numerical implementation is realized in the software library Gascoigne 3D Gascoigne3d , using iso-parametric finite elements of degree $1$ and $2$. A detailed description of the underlying numerical methods is given in Richter2017 .

4.1 Laplace equation in two and three dimensions

We consider the following problem

[TABLE]

where $\Omega$ is the unit ball in two dimensions and the unit sphere in three dimensions.

To compute errors we choose a rotationally symmetric analytical solution to (37) as

[TABLE]

with $r=\sqrt{x^{2}+y^{2}}$ in two and $r=\sqrt{x^{2}+y^{2}+z^{2}}$ in three dimensions, respectively, which results in the right hand sides

[TABLE]

For the ease of evaluations the errors, the $H^{1}$- and $L^{2}$-norms will be computed on the truncated domains

[TABLE]

see also Remark 4. We hence do not compute the errors $\|\nabla(u-u_{h})\|$ and $\|u-u_{h}\|$ on the remainders $\Omega\setminus\Omega_{r}$. Therefore we expect optimal order convergence in the spirit of Corollary 1. The restriction of the domain to an area within $\Omega_{h}$ is also by technical reasons, as the evaluation of integrals outside of the meshed area is not easily possible.

In Figures 6 and 7 we see the resulting $L^{2}$- and $H^{1}$-errors. We observe that for finer meshes, $\Upsilon$ becomes the dominating factor of the error. In particular the use of quadratic finite elements shows a strong imbalance between FE error and geometric error, which quickly dominates as seen in the left part of Fig. 7. The result is consistent with Corollary 1. As soon as the FE error is smaller than the geometry perturbation $\Upsilon$, we do not observe any further improvement of the error. In Fig. 8 we show the convergence in both norms in terms of the geometry parameter $\Upsilon$. Linear convergence is clearly observed. The apparent decay of convergence rate in case of the $L^{2}$-error in three dimensions is due to the still dominating FE error in this case.

4.2 Stokes system in two dimensions

To go beyond the Laplace problem, we investigate the behavior of the solution to the Stokes system with respect to the domain variation in two spatial dimensions. The problem is to find the velocity u and the pressure $p$ such that

[TABLE]

with homogeneous Dirichlet condition $\textbf{u}=0$ on the boundary $\partial\Omega$ and a right hand side vector f. System (38) is solved with equal-order iso-parametric finite elements using pressure stabilization by local projections, see BeckerBraack2001 .

We prescribe an analytical solution for comparison with the finite element approximation

[TABLE]

where the corresponding forcing term reads

[TABLE]

In Figure 9 we see the resulting $L^{2}$- and $H^{1}$-errors. Again we observe that $\Upsilon$ becomes the dominant factor for finer meshes. This result is not covered by the theoretical findings, however it shows that geometric uncertainty should be taken into account for the simulations of flow models.

5 Conclusions

We have demonstrated that small boundary variations have crucial impact on the result of the finite element simulations. The developed error estimates are linear with respect to the maximal distance $\Upsilon$ between the real and the approximated domains, cf. Theorem 3.2. We have illustrated the sharp nature of this bound in the computations performed in Section 4.

Particularly, in the case of first and second order approximation we observe how the relation between the mesh size $h$ and aforementioned $\Upsilon$ impact the resulting $L^{2}$- and $H^{1}$-errors. The same behavior has been demonstrated numerically for the Stokes system.

In practice we do not have control on the accuracy of the domain reconstruction. This has shown that it is worth to take into account the geometric uncertainty when deciding on the mesh-size in order to avoid unnecessary computational effort.

In this work we have focused on the Laplace problem (2). Additionally, the Stokes system has been treated numerically and it exhibits similar features. In future work we will extend this consideration to flow models, in particular the Navier-Stokes equations MiMiRi2019 . Among the additional challenges in extending the present work to the Navier-Stokes system are the consideration of the typical saddle-point structure of incompressible flow models introducing a pressure variable RannacherRichter2018 and the difficulty of nonlinearities introduced by the convective term, and thus the non-uniqueness of solutions HeywoodRannacher1990 .

Acknowledgements.

Both authors acknowledge the financial support by the Federal Ministry of Education and Research of Germany, grant number 05M16NMA as well as the GRK 2297 MathCoRe, funded by the Deutsche Forschungsgemeinschaft, grant number 314838170. Finally we thank the anonymous reviewers. Their time and effort helped us to significantly improve the manuscript.