Finite Element Error Estimates on Geometrically Perturbed Domains

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.

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  [15].

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 [26, 27]. 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 [23, 24]. For a survey on computational vascular fluid dynamics, where modeling and reconstruction related issues are discussed, we refer to [29]. 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 [3] and Neumann [2] 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 [13]. For recent developments we refer to [17, 25]. 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 [18, 32, 33]. 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 [19]. In [1, 12] 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 [10]. 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 2 gives such optimally balanced estimates that include both error contributions.

This paper is organized as follows. After this introduction, in Sect. 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 Sect. 4.

2 Mathematical Setting and Auxiliary Result

2.1 Notation

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

2.2 Laplace Equation and Domain Perturbation

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

$$\begin{aligned} -\Delta u= f \text { in }\varOmega ,\quad u=0\text { on }\partial \varOmega . \end{aligned}$$

(1)

The variational formulation of this problem is given by: find $u\in H^1_0(\varOmega )$, such that

$$\begin{aligned} (\nabla u,\nabla \phi )_\varOmega = (f,\phi )_\varOmega \quad \forall \phi \in H^1_0(\varOmega ). \end{aligned}$$

(2)

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

$$\begin{aligned} \Vert u\Vert _{H^{m+2}(\varOmega )} \le c \Vert f\Vert _{H^m(\varOmega )}, \end{aligned}$$

(3)

see e.g.  [16].

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

$$\begin{aligned} \varUpsilon :={\text {dist}}(\partial \varOmega ,\partial \varOmega _r) :=\max \{\sup _{x\in \partial \varOmega }\inf _{y\in \partial \varOmega _r}|x-y|,\sup _{y\in \partial \varOmega _r}\inf _{x\in \partial \varOmega }|x-y|\}. \end{aligned}$$

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

Assumption 1

(Domains) Let $\varOmega$ and $\varOmega _r$ be two domains with $\varOmega \cap \varOmega _r\ne \emptyset$ and with Hausdorff distance $\varUpsilon \in \mathbb {R}$. Both boundaries allow for a local $C^{m+2}$ parametrization, $m\in \mathbb {N}_0$. Let

$$\begin{aligned} S := (\varOmega _r\setminus \varOmega )\cup (\varOmega \setminus \varOmega _r). \end{aligned}$$

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\le c_{S,1}\, \varUpsilon$ with a constant $c_{S,1}\ge 0$. Following conditions hold:

• 1.a) On each rectangle R, the boundary lines $\partial \varOmega \cap R$ and $\partial \varOmega _r\cap R$ allow for unique parametrizations $g^R_\varOmega (t)$ and $g^R_{\varOmega _r}(t)$ over the base t, or $(t^1,t^2)$ for $d=3$, respectively.

• 1.b) The area of the cover is bounded by the area of the remainder S, i.e.

  $$\begin{aligned} \big | \bigcup _{i=1}^{n(S)} R_i\cap S\big | \le c_{S,2} |S|, \end{aligned}$$

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

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

Fig. 1

The domain $\varOmega$ (bold line) and its reconstruction $\varOmega _r$ (dashed). The cover of the domain remainders $S=(\varOmega \setminus \varOmega _r)\cup (\varOmega _r{\setminus } \varOmega )$ by a set of rectangles. The left configuration fulfils Assumption 1. The height $h_i$ of each rectangle $R_i$ is bounded by $c_{S,1}\varUpsilon$ and the intersections of each rectangle with the domain remainder $R_i\cap S$ do not overlap excessively. The shaded areas show the overlap. The right configuration is excluded by Assumption 1

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 Fig. 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 \varOmega _r$ and $\partial \varOmega$. Such approaches play an important role in isogeometric analysis. We refer to [34, 35] for examples on the construction of such maps.

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

Assumption 2

(Right hand side) Let $f\in H^m_\text {loc}(\mathbb {R}^d)$, i.e. $f\in H^m(G)$ for each compact subset $G\subset \mathbb {R}^d$. In addition we assume that the right hand side on $\varOmega _r$ can be bounded by the right hand side on $\varOmega$, i.e.

$$\begin{aligned} \Vert f\Vert _{H^m(\varOmega _r)}\le c \Vert f\Vert _{H^m(\varOmega )}. \end{aligned}$$

(4)

An alternative would be to use Sobolev extension theorems to extend functions $f\in H^m(\varOmega )$ from $\varOmega$ to $\varOmega _r$, see  [8].

On $\varOmega _r$ we define the solution $u_r\in H^1_0(\varOmega _r)$ to the perturbed Laplace problem

$$\begin{aligned} (\nabla u_r,\nabla \phi _r)_{\varOmega _r} = (f,\phi _r)_{\varOmega _r} \quad \forall \phi _r\in H^1_0(\varOmega _r), \end{aligned}$$

(5)

The unique solution to (5) satisfies the bound

$$\begin{aligned} \Vert u_r\Vert _{H^{m+2}(\varOmega _r)} \le c \Vert f\Vert _{H^{m}(\varOmega _r)} \le c \Vert f\Vert _{H^{m}(\varOmega )}. \end{aligned}$$

(6)

Remark 1

(Extension of the solutions) A difficulty for deriving error estimates is that u is defined on $\varOmega$ and $u_r$ on $\varOmega _r\ne \varOmega$. Since the domains do not match, u may not be defined on all of $\varOmega _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 $\varOmega _r\setminus \varOmega$ and $u_r:=0$ on $\varOmega \setminus \varOmega _r$. Globally, both functions still have the regularity $u,u_r\in H^1(\varOmega \cup \varOmega _r)$. We will use the same notation for discrete functions $u_h\in V_h$ defined on a mesh $\varOmega _h$ and extend them by zero to $\mathbb {R}^d$.

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 \mathbb {R}, \gamma >0$, $V\subset \mathbb {R}^d$ and $W\subset \mathbb {R}^d$ for $d\in \{2,3\}$ be two domains with boundaries $\partial V$ and $\partial W$ that satisfy Assumption 1 with distance

$$\begin{aligned} \gamma :={\text {dist}}(\partial V,\partial W). \end{aligned}$$

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

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

(7)

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 Fig. 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 1b). Integrating the function $\psi$ along this line gives

$$\begin{aligned} \big |\psi (x_{\partial W})\big |^2 \le 2\big |\psi (x_{\partial V})\big |^2 + 2\big | \int _{x_{\partial V}}^{x_{\partial W}}\psi '(s)\,\text {d}s\big |^2. \end{aligned}$$

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

$$\begin{aligned} |\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. \end{aligned}$$

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

$$\begin{aligned} \int |\psi (g^R_{W}(t))|^2\,\text {d}t \le \int |\psi (g^R_{V}(t))|^2\,\text {d}t +2c_S\gamma \int \int _{g^R_{V}(t)}^{g^R_{W}(t)}|\nabla \psi (s)|^2\,\text {d}s\,\text {d}t. \end{aligned}$$

(8)

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

$$\begin{aligned}&\frac{1}{\max _t \{1+|\nabla g^R_{W}(t)|^2\}} \int _{(\partial W\cap V)\cap R} |\psi |^2\,\text {d}s \nonumber \\&\quad \le \frac{1}{\min _t \{1+|\nabla g^R_{V}(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{aligned}$$

(9)

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

$$\begin{aligned} \Vert \psi \Vert ^2_{(\partial W\cap V)\cap R}\le c(\partial V,\partial W) c_S \Big (\Vert \psi \Vert ^2_{R\cap \partial V} + \gamma \Vert \nabla \psi \Vert ^2_{R\cap S}\Big ). \end{aligned}$$

(10)

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

$$\begin{aligned} \Vert \psi \Vert ^2_{\partial W\cap V}\le c(\partial V,\partial W) c_S \Big (\Vert \psi \Vert ^2_{\partial V} + \gamma \Vert \nabla \psi \Vert ^2_{V\setminus W}\Big ). \end{aligned}$$

For $\psi \in H^1_0(V)$, the term on $\partial V$ vanishes.

Fig. 2

Illustration of the proof of Lemma 3. Points in each rectangle R can be represented by local coordinates (t, s), where $0\le s\le c_S\gamma$ and the range of t depends on the size of the rectangle. The two boundary segments are given by the (smooth) parametrizations $g^R_V(t)$ and $g^R_{W}(t)$. Each line in the direction of s cuts both boundaries exactly once. In the 3d setting, the base is represented by two coordinates $\mathbf {t}=(t_1,t_2)$

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

$$\begin{aligned} |\psi (x(t,s))|^2 \le 2|\psi (g^R_V(t))|^2 + 2c_S \gamma \int \int _{x(t)}^{g^R_V(t)}|\nabla \psi |^2\,\text {d}s\text {d}t. \end{aligned}$$

We integrate over s and t to obtain

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

Summing over all rectangles gives the desired result.$\square$

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

We continue by estimating the difference between the solutions of the Laplace equations on $\varOmega$ and on $\varOmega _r$.

Lemma 4

Let $\varOmega , \varOmega _r \in \mathbb {R}^d$ with $\partial \varOmega ,\partial \varOmega _r \in C^{m+2}$ satisfying ${\text {dist}}(\partial \varOmega ,\partial \varOmega _r) < \varUpsilon$ as well as Assumption 1. Furthermore, let $f\in L^2_\text {loc}(\mathbb {R}^d)$ satisfy Assumption 2 and let $f_r:=f|_{\varOmega _r}$. For the solutions $u\in H^1_0(\varOmega )\cap H^2(\varOmega )$ and $u_r\in H^1_0(\varOmega _r)\cap H^2(\varOmega _r)$ to (2) and (5) respectively, it holds that

$$\begin{aligned} \Vert u-u_r\Vert _\varOmega + \varUpsilon ^\frac{1}{2}\Vert \nabla (u-u_r)\Vert _{\varOmega } \le c \varUpsilon \Vert f\Vert _{\varOmega \cup \varOmega _r}. \end{aligned}$$

Proof

(i) We continuously extend u and $u_r$ by zero to $\mathbb {R}^d$, c.f. Remark 1, such that $u-u_r\in H^1(\varOmega \cup \varOmega _r)$ is well defined.

We separate the domains of integration and integrate by parts

$$\begin{aligned} \begin{aligned} \Vert \nabla (u-u_r)\Vert _\varOmega ^2&= \big (\nabla (u-u_r),\nabla (u-u_r)\big )_{\varOmega \cap \varOmega _r} + \big (\nabla (u-u_r),\nabla (u-u_r)\big )_{\varOmega \setminus \varOmega _r}\\&=-\big (\Delta (u-u_r),u-u_r\big )_\varOmega + \langle \partial _n(u-u_r),u-u_r\rangle _{\partial (\varOmega \cap \varOmega _r)} \\&\quad +\,\langle \partial _n(u-u_r),u-u_r\rangle _{\partial (\varOmega \setminus \varOmega _r)}. \end{aligned} \end{aligned}$$

(11)

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

$$\begin{aligned} \begin{aligned} \Vert \nabla (u-u_r)\Vert _\varOmega ^2&= -\big (\Delta (u-u_r),u-u_r\big )_\varOmega + \langle \partial _n(u-u_r),u-u_r\rangle _{\partial \varOmega } \\&\quad +\,\langle [\partial _n (u-u_r)],u-u_r\rangle _{\partial \varOmega _r\cap \varOmega }. \end{aligned} \end{aligned}$$

(12)

In $\varOmega \cap \varOmega _r$ it holds $f=f_r$ and hence (weakly) $-\Delta (u-u_r)=0$, such that

$$\begin{aligned} -(\Delta (u-u_r),u-u_r)_\varOmega = -(\Delta (u-u_r),u-u_r)_{\varOmega \cap \varOmega _r} -(\Delta u,u)_{\varOmega \setminus \varOmega _r} = (f,u)_{\varOmega \setminus \varOmega _r}. \end{aligned}$$

(13)

On $\partial \varOmega$ it holds $u=0$ and on $\partial \varOmega _r\cap \varOmega$ it holds $u_r=0$. Further, since $u\in H^2(\varOmega )$ it holds that $[\partial _n u]=0$ on $\partial \varOmega _r\cap \varOmega$. Finally, $u_r=0$ on $\varOmega \setminus \varOmega _r$, such that the boundary terms reduce to

$$\begin{aligned}&\langle \partial _n(u-u_r),u-u_r\rangle _{\partial \varOmega } + \langle [\partial _n(u-u_r)],u-u_r\rangle _{\partial \varOmega _r\cap \varOmega } \nonumber \\&\quad =-\langle \partial _n(u-u_r),u_r\rangle _{\partial \varOmega \cap \varOmega _r} -\langle \partial _n u_r,u\rangle _{\partial \varOmega _r\cap \varOmega }. \end{aligned}$$

(14)

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

$$\begin{aligned} \Vert \nabla (u-u_r)\Vert _\varOmega ^2\le & {} \Vert f\Vert _{\varOmega \setminus \varOmega _r} \Vert u\Vert _{\varOmega \setminus \varOmega _r} \nonumber \\&+\Vert \partial _n(u-u_r)\Vert _{\partial \varOmega \cap \varOmega _r} \Vert u_r\Vert _{\partial \varOmega \cap \varOmega _r} +\Vert \partial _nu\Vert _{\partial \varOmega _r\cap \varOmega } \Vert u\Vert _{\partial \varOmega _r\cap \varOmega }. \end{aligned}$$

(15)

Since $u,u_r\in H^2(\varOmega \cap \varOmega _r)$, the trace inequality gives

$$\begin{aligned} \Vert \nabla (u-u_r)\Vert _\varOmega ^2\le & {} \Vert f\Vert _{\varOmega \setminus \varOmega _r} \Vert u\Vert _{\varOmega \setminus \varOmega _r} \nonumber \\&+\, c\big (\Vert u\Vert _{H^2(\varOmega )}+\Vert u_r\Vert _{H^2(\varOmega _r)}\big ) \big ( \Vert u_r\Vert _{\partial \varOmega \cap \varOmega _r} + \Vert u\Vert _{\partial \varOmega _r\cap \varOmega }\big ). \end{aligned}$$

(16)

Applying Lemma 3 twice: to $\psi =u$ and to $\psi =\nabla u$ (same for $u_r)$, and extending the norms from $\varOmega \setminus \varOmega _r$ to $\varOmega$ and from $\varOmega _r\setminus \varOmega$ to $\varOmega _r$ give the bounds

$$\begin{aligned} \begin{aligned} \Vert u\Vert _{\partial \varOmega _r\cap \varOmega }&\le c \varUpsilon ^\frac{1}{2} \Vert \nabla u\Vert _{\varOmega \setminus \varOmega _r} \le c\varUpsilon \Big ( \Vert \nabla u\Vert _{\partial \varOmega } + \varUpsilon ^\frac{1}{2} \Vert u\Vert _{H^2(\varOmega )}\Big ), \\ \Vert u_r\Vert _{\partial \varOmega \cap \varOmega _r}&\le c\varUpsilon ^\frac{1}{2} \Vert \nabla u_r\Vert _{\varOmega _r\setminus \varOmega }\le c \varUpsilon \Big ( \Vert \nabla u_r\Vert _{\partial \varOmega _r} + \varUpsilon ^\frac{1}{2} \Vert u_r\Vert _{H^2(\varOmega _r)}\Big ). \end{aligned} \end{aligned}$$

(17)

With the trace inequality and the a priori estimates $\Vert u\Vert _{H^2(\varOmega )}\le c \Vert f\Vert _{\varOmega }$ and $\Vert u_r\Vert _{H^2(\varOmega _r)} \le c \Vert f_r\Vert _{\varOmega _r} \le c \Vert f\Vert _{\varOmega }$ we obtain the bounds

$$\begin{aligned} \Vert u\Vert _{\partial \varOmega _r\cap \varOmega }\le c \varUpsilon \Vert f\Vert _{\varOmega },\quad \Vert u_r\Vert _{\partial \varOmega \cap \varOmega _r}\le c \varUpsilon \Vert f\Vert _{\varOmega }. \end{aligned}$$

(18)

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

$$\begin{aligned} \Vert u\Vert _{\varOmega \setminus \varOmega _r}\le c \varUpsilon \Vert \nabla u\Vert _{\varOmega \setminus \varOmega _r} \le c \varUpsilon ^\frac{3}{2}\Big (\Vert u\Vert _{H^2(\varOmega )} + \varUpsilon ^\frac{1}{2} \Vert u\Vert _{H^2(\varOmega )}\Big )\le c \varUpsilon ^\frac{3}{2} \Vert f\Vert _{\varOmega }. \end{aligned}$$

(19)

We can then estimate $\Vert f\Vert _{\varOmega \setminus \varOmega _r}\le \Vert f\Vert _{\varOmega }$ by extending to the complete domain. Combining (16) with (18) and (19) we obtain the estimate

$$\begin{aligned} \Vert \nabla (u-u_r)\Vert ^2_{\varOmega }\le c\Big (\varUpsilon ^\frac{3}{2}+ \varUpsilon \Big ) \Vert f\Vert ^2_{\varOmega }, \end{aligned}$$

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

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

$$\begin{aligned} z\in H^1_0(\varOmega ):\quad -\Delta z = \frac{u-u_r}{\Vert u-u_r\Vert _\varOmega }\text { in }\varOmega , \end{aligned}$$

which allows for a unique solution satisfying the a-priori bound $\Vert z\Vert _{H^2(\varOmega )}\le c_s$ with the stability constant $c_s<\infty$. Testing with $u-u_r$ and integrating by parts twice gives

$$\begin{aligned} \begin{aligned} \Vert u-u_r\Vert _\varOmega&= - (z,\Delta (u-u_r))_{\varOmega } + \langle z,\partial _n(u-u_r)\rangle _{\partial \varOmega } \\&\quad +\, \langle z,[\partial _n(u-u_r)]\rangle _{\partial \varOmega _r\cap \varOmega } -\langle \partial _n z,u-u_r\rangle _{\partial \varOmega }. \end{aligned} \end{aligned}$$

It holds $z=0$ and $u=0$ on $\partial \varOmega$, $[\partial _n u]=0$ on $\partial \varOmega _r\cap \varOmega$ and $-\Delta (u-u_r)=0$ in $\varOmega \cap \varOmega _r$ such that we get

$$\begin{aligned} \begin{aligned} \Vert u-u_r\Vert _\varOmega&= (z,f)_{\varOmega \setminus \varOmega _r} -\langle z,\partial _n u_r\rangle _{\partial \varOmega _r\cap \varOmega } +\langle \partial _n z,u_r\rangle _{\partial \varOmega }\\&\le \Vert z\Vert _{\varOmega \setminus \varOmega _r} \Vert f\Vert _{\varOmega \setminus \varOmega _r} + \Vert z\Vert _{\partial \varOmega _r\cap \varOmega } \Vert \partial _n u_r\Vert _{\partial \varOmega _r\cap \varOmega } + \Vert \partial _n z\Vert _{\partial \varOmega } \Vert u_r\Vert _{\partial \varOmega }. \end{aligned} \end{aligned}$$

The boundary terms $\Vert z\Vert _{\partial \varOmega _r\cap \varOmega }$ and $\Vert u_r\Vert _{\partial \varOmega }$ are estimated with Lemma 3, the normal derivatives by the trace inequality and the terms on $\varOmega \setminus \varOmega _r$ by (7)

$$\begin{aligned} \Vert u-u_r\Vert _\varOmega \le c \varUpsilon ^\frac{3}{2} \Vert z\Vert _{H^2(\varOmega )} \Vert f\Vert _\varOmega + c\varUpsilon \Vert z\Vert _{H^2(\varOmega )} \Vert u_r\Vert _{H^2(\varOmega _r)} + c\varUpsilon \Vert z\Vert _{H^2(\varOmega )}\Vert u_r\Vert _{H^2(\varOmega _r)}. \end{aligned}$$

The $L^2$-norm estimate follows by using the bounds $\Vert u\Vert _{H^2(\varOmega )}\le c\Vert f\Vert _{\varOmega }$, $\Vert u_r\Vert _{H^2(\varOmega _r)}\le c \Vert f\Vert _\varOmega$ and $\Vert z\Vert _{H^2(\varOmega )}\le c$.$\square$

Remark 2

The estimate $\Vert f\Vert _{\varOmega \setminus \varOmega _r}\le c \Vert f\Vert _{\varOmega }$ is not optimal. Further powers of $\varUpsilon$ are easily generated at the cost of a higher right hand side regularity. Also, the estimate $\Vert \partial _n (u-u_r)\Vert \le c (\Vert u\Vert _{H^2(\varOmega )}+\Vert u_r\Vert _{H^2(\varOmega _r)})$ by Cauchy Schwarz and the trace inequality could be enhanced to produce powers of $\varUpsilon$. The limiting term in (12) however is the boundary integral $|\langle \partial _n u_r,u\rangle _{\partial \varOmega _r\cap \varOmega }|=\mathcal {O}(\varUpsilon ^\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 $\varOmega \cap \varOmega _r$ only and that allows us to improve the order to $\mathcal {O}(\varUpsilon )$ in the $H^1$-case by avoiding exactly this boundary integral.

3 Discretization

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

We partition $\varOmega _r$ into a parametric triangulation $\varOmega _h$, consisting of open elements $T\subset \mathbb {R}^d$. Each element $T\in \varOmega _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 Sect. 4 are based on quadrilateral meshes. The map $T_T:\hat{T}\rightarrow T$ is a polynomial of degree $r\in \mathbb {N}$. We will consider iso-parametric finite element spaces, that are based on polynomials of the same degree $r\ge 1$. In the following we assume structural and shape regularity of the mesh such that standard interpolation estimates

$$\begin{aligned} \begin{aligned} \Vert \nabla ^k (u_r-I_h u_r)\Vert _T&\le c h^{r+1-k}\Vert u_r\Vert _{H^{r+1}(T)},\quad k=0,\dots ,r\le m,\\ \Vert \nabla ^k (u_r-I_h u_r)\Vert _{\partial T}&\le c h^{r+\frac{1}{2}-k}\Vert u_r\Vert _{H^{r+1}(T)},\quad k=0,\dots ,r\le m, \end{aligned} \end{aligned}$$

(20)

will hold for all elements $T\in \varOmega _h$, c.f.  [7]. The discretization parameter h represents the size of the largest element in the mesh. See  [31, 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.

$$\begin{aligned} \hat{P} \;\widehat{=}\; Q^r := {\text {span}}\{x_1^{\alpha _1}\cdots x_d^{\alpha _d}\;:\; 0\le \alpha _1,\ldots ,\alpha _d\le r\} \end{aligned}$$

on quadrilateral and hexahedral meshes. Then, the finite element space $V_h^r$ on the mesh $\varOmega _h$ is defined as

$$\begin{aligned} V_h^r=\{\phi _h\in C(\bar{\varOmega }_h)\;:\; \phi _h\circ T_T\in \hat{P}\text { on every }T\in \varOmega _h\}. \end{aligned}$$

This parametric finite element space does not exactly match the domain $\varOmega _r$. Given an iso-parametric mapping of degree r it holds ${\text {dist}}(\partial \varOmega _r,\partial \varOmega _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  [14], optimal interpolation and finite element error estimates have been presented in  [11, Section 4.4]. The case of higher order elements with optimal order energy norm estimates is covered in  [21].

From  [31, 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 1

Let $m\in \mathbb {N}_{0}$ and let $\varOmega _r$ be a domain with a boundary that allows for a parametrization of degree $m+2$. Let $f_r\in H^{m}(\varOmega _r)$ and $u_h\in V_h^r\cap H^1_0(\varOmega _h)$ be the iso-parametric finite element discretization of degree $1\le r\le m+1$

$$\begin{aligned} (\nabla u_h,\nabla \phi _h)_{\varOmega _h}= (f_r,\phi _h)_{\varOmega _h}\quad \forall \phi _h\in V_h^{r}. \end{aligned}$$

It holds

$$\begin{aligned} \Vert u_r-u_h\Vert _{H^1(\varOmega _r)} \le ch^r \Vert f_r\Vert _{H^{r-1}(\varOmega _r)},\quad \Vert u_r-u_h\Vert _{\varOmega _r} \le ch^{r+1} \Vert f_r\Vert _{H^{r-1}(\varOmega _r)}. \end{aligned}$$

We formulated the error estimate on the domain $\varOmega _r$ although the finite element functions are given on $\varOmega _h$ only. To give Theorem 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 2

Let $m\in \mathbb {N}_{0}$, $\varOmega$ and $\varOmega _r$ be domains with $C^{m+2}$ boundary, distance $\varUpsilon$ and that satisfy Assumption 1. Let $\varOmega _h$ be the iso-parametric mesh of $\varOmega _r$ with degree $1\le r\le m+1$ and let $f\in H^{r-1}_{loc}(\mathbb {R}^d)$ satisfy Assumption 2. For the finite element error between the fully discrete solution $u_h\in V_h^{r}$

$$\begin{aligned} (\nabla u_h,\nabla \phi _h)_{\varOmega _h}= (f,\phi _h)_{\varOmega _h} \quad \forall \phi _h\in V_h^{r} \end{aligned}$$

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

$$\begin{aligned} \Vert u-u_h\Vert _{H^1(\varOmega )} \le c \big (\varUpsilon ^\frac{1}{2} + h^r \big )\Vert f\Vert _{H^{r-1}(\varOmega )}, \end{aligned}$$

as well as

$$\begin{aligned} \Vert u-u_h\Vert _{\varOmega } \le c (\varUpsilon + h^{r+1}) \Vert f\Vert _{H^{r-1}(\varOmega )}. \end{aligned}$$

Proof

(i) We start with the $H^1$ error. Inserting $\pm u_r$ and extending the finite element error $u_r-u_h$ from $\varOmega$ to $\varOmega _r$, where a small remainder appears, we have

$$\begin{aligned} \Vert \nabla (u-u_h)\Vert _\varOmega ^2 \le 2\Big ( \Vert \nabla (u-u_r)\Vert _\varOmega ^2 + \Vert \nabla (u_r-u_h)\Vert _{\varOmega _r}^2 + \Vert \nabla (u_r-u_h)\Vert ^2_{\varOmega \setminus \varOmega _r}\Big ). \end{aligned}$$

(21)

The first and the second term on the right hand side are estimated by Lemma 4 and Theorem 1 and, since $u_r=0$ on $\varOmega \setminus \varOmega _r$, we obtain

$$\begin{aligned} \Vert \nabla (u-u_h)\Vert _\varOmega ^2 \le c \Big (\varUpsilon + h^{2r}\Big ) \Vert f\Vert _{H^{r-1}(\varOmega )}^2 + 2 \Vert \nabla u_h\Vert ^2_{\varOmega \setminus \varOmega _r}. \end{aligned}$$

(22)

We continue with the remainder $\nabla u_h$ on $\varOmega \setminus \varOmega _r$, which is non-zero on $\varOmega _h$ only

$$\begin{aligned} \Vert \nabla u_h\Vert _{\varOmega \setminus \varOmega _r}^2 = \Vert \nabla u_h\Vert ^2_{(\varOmega _h\setminus \varOmega _r)\cap (\varOmega \setminus \varOmega _r)}. \end{aligned}$$

This remaining stripe has the width

$$\begin{aligned} \gamma _{h,\varUpsilon }:={\mathscr {O}}(\min \{h^{r+1},\varUpsilon \}), \end{aligned}$$

and we apply Lemma 3 to get

$$\begin{aligned} \Vert \nabla u_h\Vert _{(\varOmega \setminus \varOmega _r)\cap (\varOmega _h\setminus \varOmega _r)}^2\le c \gamma _{h,\varUpsilon } \Vert \nabla u_h\Vert _{\partial \varOmega _r}^2 +c \gamma _{h,\varUpsilon }^2 \Vert \nabla ^2 u_h\Vert _{\varOmega _h\setminus \varOmega _r}^2, \end{aligned}$$

(23)

where the second derivative $\nabla ^2 u_h$ is understood element wise. This term is extended to $\varOmega _h$ and with the inverse estimate and the a priori estimate for the discrete solution we obtain with $\gamma _{h,\varUpsilon }^2={\mathscr {O}}(h^{2r+2})$ that

$$\begin{aligned} \gamma _{h,\varUpsilon }^2 \Vert \nabla ^2 u_h\Vert _{\varOmega _h\setminus \varOmega _r}^2 \le c_{inv} \gamma _{h,\varUpsilon }^2\ h^{-2} \Vert \nabla u_h\Vert _{\varOmega _h}^2 \le c_{inv} h^{2r} \Vert f\Vert _{\varOmega _h}^2\le c h^{2r} \Vert f\Vert ^2_{\varOmega }. \end{aligned}$$

(24)

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

$$\begin{aligned} \gamma _{h,\varUpsilon } \Vert \nabla u_h\Vert ^2_{\partial \varOmega _r}\le c\gamma _{h,\varUpsilon }\big ( \Vert \nabla u_r\Vert ^2_{\partial \varOmega _r} +\Vert \nabla (u_r-I_h u_r)\Vert ^2_{\partial \varOmega _r} +\Vert \nabla (u_h-I_h u_r)\Vert ^2_{\partial \varOmega _r}\big ). \end{aligned}$$

(25)

Here, the first and last terms are estimated with the trace inequalities and, in the case of the discrete term with the inverse inequality¹, followed by adding $\pm u_r$ we get

$$\begin{aligned}&\gamma _{h,\varUpsilon } \Vert \nabla u_h\Vert ^2_{\partial \varOmega _r}\le c\varUpsilon \Vert f\Vert _{L^2(\varOmega )}^2 \nonumber \\&\quad + ch^{r+1} \Vert \nabla (u_r-I_h u_r)\Vert ^2_{\partial \varOmega _r}\nonumber \\ +ch^r \Vert \nabla (u_r-I_h u_r)\Vert ^2_{\varOmega _r} + h^r \Vert \nabla (u_h-u_r)\Vert ^2_{\varOmega _r}. \end{aligned}$$

(26)

We used both $\gamma _{h,\varUpsilon }={\mathscr {O}}(h^{r+1})$ and $\gamma _{h,\varUpsilon }={\mathscr {O}}(\varUpsilon )$.

Then, collecting all terms in (22)–(26) and using the interpolation estimates as well as Theorem 1 we finally get

$$\begin{aligned} \Vert \nabla (u-u_h)\Vert ^2_\varOmega \le c \big (\varUpsilon + h^{2r}\big ) \Vert f\Vert _{H^{r-1}(\varOmega )}^2 +ch^{3r-1}\Vert f\Vert _{H^{r-1}(\varOmega )}^2, \end{aligned}$$

(27)

which shows the a priori estimate since $3r-1\le 2r$ for all $r\ge 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 \varOmega _h$, where $u_h=0$, can be applied, i.e.

$$\begin{aligned} \Vert u-u_h\Vert ^2_\varOmega \le c \Big ( \Vert u-u_r\Vert ^2_\varOmega + \Vert u_r-u_h\Vert _{\varOmega _r}^2 +\varUpsilon ^2 \Vert \nabla u_h\Vert ^2_{\varOmega \setminus \varOmega _r} \Big ). \end{aligned}$$

The $L^2$-estimate directly follows with Lemma 4, Theorem 1 and by the a priori estimate $\Vert \nabla u_h\Vert ^2_{\varOmega \setminus \varOmega _r}\le c \Vert \nabla u_h\Vert ^2_{\varOmega _h}\le c\Vert f\Vert ^2_{\varOmega _r}$.$\square$

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(\varOmega )$ and $u_r\in H^2(\varOmega _r)$, is relatively straightforward. In this case, $\varOmega _h$ fits $\varOmega _r$ such that the finite element error $u_r-u_h$ can be estimated with the standard a priori result $\Vert u_r-u_h\Vert +h\Vert \nabla (u_r-u_h)\Vert \le c \Vert f\Vert$. 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 Fig. 1: All rectangles are split in such a way that the corners of $\partial \varOmega$ and $\varOmega _r$ are cut by the edges of rectangles. This allows to derive the optimal error estimates $\Vert u-u_h\Vert _{H^1(\varOmega )}= {\mathscr {O}}(\varUpsilon ^\frac{1}{2}+h)$ and $\Vert u-u_h\Vert _{H^1(\varOmega )}= {\mathscr {O}}(\varUpsilon +h^2)$. In three dimensions, such a simple refinement of the cover is not possible and the extension to polygonal domains is more involved.

Remark 4

(Optimality of the estimates) Two ingredients govern the error estimates:

• 1. A geometrical error of order $\mathcal {O}(\varUpsilon ^\frac{1}{2})$ and $\mathcal {O}(\varUpsilon )$, that describes the discrepancy between $\varOmega$ and $\varOmega _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 Fig. 3, namely $-\Delta u=4$ on the unit disc $\varOmega =B_1(0)$ and $-\Delta u_r=4$ on the shifted domain $\varOmega _r = B_1(\varUpsilon )$. The errors in $H^1$ norm and $L^2$ norms expressed on the complete domain $\varOmega$ are estimated by

  $$\begin{aligned} \Vert u-u_r\Vert _\varOmega =\sqrt{\pi }\varUpsilon + \mathcal {O}(\varUpsilon ^3),\quad \Vert \nabla (u-u_r)\Vert _\varOmega =\sqrt{8\varUpsilon } + \mathcal {O}(\varUpsilon ). \end{aligned}$$

  A closer analysis shows that the main error – in the $H^1$-case – occurs on the small shaded stripe $\varOmega \setminus \varOmega _r$ such that

  $$\begin{aligned} \Vert \nabla (u-u_r)\Vert _{\varOmega \setminus \varOmega _r} = \mathcal {O}(\varUpsilon ^\frac{1}{2}),\quad \Vert \nabla (u-u_r)\Vert _{\varOmega \cap \varOmega _r} = \mathcal {O}(\varUpsilon ), \end{aligned}$$

  while the $L^2$-error in $\varOmega \cap \varOmega _r$ is optimal

  $$\begin{aligned} \Vert u-u_r\Vert _{\varOmega \setminus \varOmega _r} = \mathcal {O}(\varUpsilon ^\frac{3}{2}),\quad \Vert u-u_r\Vert _{\varOmega \cap \varOmega _r} = \mathcal {O}(\varUpsilon ). \end{aligned}$$

• 2. The usual Galerkin error $\Vert u_r-u_h\Vert _{\varOmega _r} + h\Vert \nabla (u-u_r)\Vert _{\varOmega _r} = \mathcal {O}(h^{r+1})$ of iso-parametric finite element approximations contributes to the overall error. For $\varOmega =\varOmega _r$, i.e. $\varUpsilon =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.

Fig. 3

Illustration concerning Remark 4. The error estimates for $u-u_h$ are optimal, if the error is evaluated on $\varOmega$. The lowest order terms $\mathcal {O}(\varUpsilon ^\frac{1}{2})$ appear in the shaded area $\varOmega \setminus \varOmega _r$ where $u_r$ and (most of) $u_h$ are zero

In Sect. 4 we discuss the difficulty of measuring errors on an unknown domain $\varOmega$. 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 $\varOmega \setminus \varOmega _r$, where no finite element mesh is given. These remainders contribute the lowest order parts $\varUpsilon ^\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 $\varUpsilon$ in the $H^1$-norm error. In addition to the previous setting we require a regular map $T_r:\varOmega \rightarrow \varOmega _r$ between the two domains. By pulling back $\varOmega _r$ to $\varOmega$ via this map a Jacobian arises that controls the geometrical error and that hence has to be controllable by $\varUpsilon$.

Corollary 1

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

$$\begin{aligned} T_r:\varOmega \rightarrow \varOmega _r \end{aligned}$$

satisfying

$$\begin{aligned} \Vert I-{\text {det}}(\nabla T_r) \nabla T_r^{-1}\nabla T_r^{-T}\Vert _{L^\infty (\varOmega )} =\mathcal {O}(\varUpsilon ). \end{aligned}$$

(28)

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

$$\begin{aligned} f\in W^{1,\infty }_{loc}(\mathbb {R}^d)\cap H^{r-1}_{loc}(\mathbb {R}^d) \end{aligned}$$

(29)

and let the solution satisfy

$$\begin{aligned} \Vert u\Vert _{W^{2,\infty }(\varOmega )} + \Vert u_r\Vert _{W^{2,\infty }(\varOmega _r)} \le c. \end{aligned}$$

(30)

Then, it holds

$$\begin{aligned} \Vert \nabla (u-u_h)\Vert _{\varOmega \cap \varOmega _r\cap \varOmega _h}\le c \big ( \varUpsilon + h^r \big ). \end{aligned}$$

Proof

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

$$\begin{aligned} \Vert \nabla (u-u_h)\Vert _{\varOmega \cap \varOmega _r\cap \varOmega _h}\le \Vert \nabla (u-u_r)\Vert _{\varOmega \cap \varOmega _r}+ \Vert \nabla (u_r-u_h)\Vert _{\varOmega _r\cap \varOmega _h}. \end{aligned}$$

(31)

An optimal order estimate of the finite element error

$$\begin{aligned} \Vert \nabla (u_r-u_h)\Vert _{\varOmega _r\cap \varOmega _h}\le \Vert \nabla (u_r-u_h)\Vert _{\varOmega _r}=\mathcal {O}(h^r) \end{aligned}$$

(32)

is given in Theorem 1.

To estimate the first term of the right hand side of (31) we introduce the function

$$\begin{aligned} \hat{u}_r(x):=u_r(T_r(x)), \end{aligned}$$

which satisfies $\hat{u}_r\in H^1_0(\varOmega )$ and solves the problem

$$\begin{aligned} (J_r F_r^{-1}F_r^{-T}\nabla \hat{u}_r,\nabla \hat{\phi }_r)_{\varOmega } = (\hat{f}_r,\hat{\phi }_r)\quad \forall \hat{\phi }_r\in H^1_0(\varOmega ), \end{aligned}$$

where $\hat{f}_r(x):=f(T_r(x))$ and where $F_r:=\nabla T_r$ and $J_r:={\text {det}}(F_r)$. See  [31, 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

$$\begin{aligned} \Vert \nabla (u-u_r)\Vert _{\varOmega \cap \varOmega _r}\le \Vert \nabla (u-\hat{u}_r)\Vert _{\varOmega \cap \varOmega _r}+ \Vert \nabla (\hat{u}_r-u_r)\Vert _{\varOmega \cap \varOmega _r}. \end{aligned}$$

(33)

We introduce the notation $e_r:=u-\hat{u}_r$, extend the first term from ${\varOmega \cap \varOmega _r}$ to $\varOmega$ and insert $\pm J_rF_r^{-1}F_r^{-T}\nabla \hat{u}_r$ which gives

$$\begin{aligned}&\Vert \nabla (u-\hat{u}_r)\Vert _{\varOmega \cap \varOmega _r}^2\le \Vert \nabla (u-\hat{u}_r)\Vert _{\varOmega }^2 \nonumber \\&=(\nabla u,\nabla e_r)_{\varOmega } - (J_rF_r^{-1}F_r^{-T}\nabla \hat{u}_r,\nabla e_r)_{\varOmega } +(J_rF_r^{-1}F_r^{-T}\nabla \hat{u}_r,\nabla e_r)_{\varOmega } -(\nabla \hat{u}_r,\nabla e_r)_{\varOmega }\nonumber \\&=(f-\hat{f}_r,e_r)_{\varOmega } + ( [J_rF_r^{-1}F_r^{-T}-I]\nabla \hat{u}_r,\nabla e_r)_{\varOmega }\nonumber \\&\quad \le \Vert f-\hat{f}_r\Vert _{\varOmega } c \Vert \nabla e_r\Vert _\varOmega + \Vert [J_rF_r^{-1}F_r^{-T}-I]\Vert _{L^\infty (\varOmega )} \Vert \nabla e_r\Vert _\varOmega , \end{aligned}$$

(34)

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

$$\begin{aligned} |f(x)-\hat{f}_r(x)| = |f(x)-f(T_r(x))|= |\nabla f(\xi )\cdot (T_r(x)-x)| \le \varUpsilon |\nabla f(\xi )|, \end{aligned}$$

(35)

where $\xi \in \varOmega$ is some point on the line from x to $T_r(x)$. We take the square and integrate over $\varOmega$ to get the estimate

$$\begin{aligned} \Vert f-\hat{f}_r\Vert _\varOmega \le c \varUpsilon \Vert f\Vert _{W^{1,\infty }(\varOmega _\varUpsilon )}, \end{aligned}$$

(36)

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

$$\begin{aligned} \Vert \nabla (\hat{u}_r-u_r)\Vert _M \le c\varUpsilon \Vert u_r\Vert _{W^{2,\infty }(\varOmega \cap \varOmega _r)} \le c \varUpsilon . \end{aligned}$$

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

Unfortunately this corollary can not be applied universally as the existence of a suitable map $T_r:\varOmega \rightarrow \varOmega _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}:\varOmega \rightarrow \mathbb {R}^2$

$$\begin{aligned} T_r(x) = x+\hat{d}(x),\quad F_r(x) = I + \nabla \hat{d}(x). \end{aligned}$$

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

$$\begin{aligned} \Vert J_r\Vert _{L^\infty (\varOmega )} = 1+\mathcal {O}(\varUpsilon ),\quad \Vert I-J_rF_r^{-1}F_r^{-T}\Vert _{L^\infty (\varOmega )} = \mathcal {O}(\varUpsilon ). \end{aligned}$$

While the assumption $|\hat{d}|=\mathcal {O}(\varUpsilon )$ is easy to satisfy since ${\text {dist}}(\partial \varOmega ,\partial \varOmega _r)\le \varUpsilon$, the condition $|\nabla \hat{d}|=\mathcal {O}(\varUpsilon )$ 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 $\varOmega$ be the unit sphere, $\varOmega _r$ be an ellipse

$$\begin{aligned} \varOmega = \{x\in \mathbb {R}^2\,:\, x_1^2+x_2^2<1\},\quad \varOmega _r = \{x\in \mathbb {R}^2\,:\, (1+\varUpsilon )^2x_1^2+(1+\varUpsilon )^{-2}x_2^2<1\}. \end{aligned}$$

It holds ${\text {dist}}(\partial \varOmega ,\partial \varOmega _r)\le \varUpsilon$ and we define the map $T_r:\varOmega \rightarrow \varOmega _r$ by

$$\begin{aligned} T_r(x) = \begin{pmatrix} (1+\varUpsilon )^{-1} x_1 \\ (1+\varUpsilon ) x_2, \end{pmatrix},\quad F_r=\nabla T_r = \begin{pmatrix} (1+\varUpsilon )^{-1}&{}0\\ 0&{}(1+\varUpsilon ) \end{pmatrix},\quad J_r=1. \end{aligned}$$

This map satisfies the assumptions of the corollary

$$\begin{aligned} I-J_r F_r^{-1}F_r^{-T} =\varUpsilon (\varUpsilon +2) \begin{pmatrix} -1&{}0\\ 0&{} (1+\varUpsilon )^{-2} \end{pmatrix},\quad \Vert I-J_r F_r^{-1}F_r^{-T}\Vert _{\infty }=2\varUpsilon +\varUpsilon ^2. \end{aligned}$$

Fig. 4

Illustration of an example for the application of Corollary 1

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 $\varUpsilon$. Moreover, we numerically extend the analytical predictions and show that a similar behavior holds for the Stokes system.

We consider $\varOmega$ to be a unit ball in two and three dimensions and define a family of perturbed domains $\varOmega _\varUpsilon$, with the amplitude of the perturbation being dependent on the coefficient $\varUpsilon$, cf. Fig. 5.

Fig. 5

Sketch of the computational domains w.r.t. the parameter $\varUpsilon$ in two dimensions (left) and for $\varUpsilon =0.1$ in three dimensions (right)

In two dimensions, the boundary of the domain $\varOmega _\varUpsilon$ is given in polar coordinates $(\rho ,\varphi )$ by

$$\begin{aligned} \partial \varOmega _\varUpsilon = \{(1 - \varUpsilon /5+\varUpsilon \sin (8\varphi ),\varphi ) \text { for } \varphi \in [0,2\pi )\}, \end{aligned}$$

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

$$\begin{aligned} \partial \varOmega _\varUpsilon = \{(1 - \varUpsilon /5+\varUpsilon \sin (3\varphi )\sin (3\theta ),\theta ,\varphi ) \text { for } \theta \in [0,\pi ), \varphi \in [0,2\pi )\}. \end{aligned}$$

For computations we take

$$\begin{aligned} \varUpsilon \in \{0,0.0125,0.025,0.05,0.1 \}. \end{aligned}$$

In order to illustrate the convergence result from Theorem 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 $\varOmega _\varUpsilon$, with a mesh size h, by $\varOmega _{h,\varUpsilon }$.

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

4.1 Laplace Equation in Two and Three Dimensions

We consider the following problem

$$\begin{aligned} -\Delta u = f \text { in }\varOmega ,\quad u = 0 \text { on }\partial \varOmega , \end{aligned}$$

(37)

where $\varOmega$ 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

$$\begin{aligned} u(r)= - \cos \left( \frac{\pi }{2} r\right) \end{aligned}$$

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

$$\begin{aligned} f_{2d}(r) = \frac{\pi }{2r} \sin \left( \frac{\pi }{2} r\right) + \frac{\pi ^2}{4} \cos \left( \frac{\pi }{2} r\right) ,\quad f_{3d}(r) = \frac{\pi }{r} \sin \left( \frac{\pi }{2} r\right) + \frac{\pi ^2}{4} \cos \left( \frac{\pi }{2} r\right) . \end{aligned}$$

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

$$\begin{aligned} \begin{aligned} \varOmega '_{2d}&=\{(\varphi ,\rho ) \text { for } \varphi \in [0,2\pi ) \text { and } \rho \in (0,0.88) \},\\ \varOmega '_{3d}&= \{(\varphi ,\theta , \rho ) \text { for }\theta \in [0,\pi ), \varphi \in [0,2\pi ) \text { and } \rho \in (0,0.88) \}, \end{aligned} \end{aligned}$$

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

In Figs. 7 and 6 we see the resulting $L^2$- and $H^1$-errors. We observe that for finer meshes, $\varUpsilon$ 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. 6. The result is consistent with Corollary 1. As soon as the FE error is smaller than the geometry perturbation $\varUpsilon$, 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 $\varUpsilon$. 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.

Fig. 6

$L^2$- and $H^1$-errors w.r.t. mesh-size $h_{max}$ for varying parameter $\varUpsilon$ computed for the Laplace problem in two-dimensions with FE. Left: linear finite elements. Right: quadratic finite elements

Fig. 7

$L^2$- and $H^1$-errors w.r.t. mesh-size $h_{max}$ for varying parameter $\varUpsilon$ computed for the Laplace problem in three-dimensions with linear finite elements

Fig. 8

$L^2$- and $H^1$-errors w.r.t. parameter $\varUpsilon$ computed for the Laplace problem in two and three-dimensions with linear and quadratic finite elements

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 $\mathbf{u}$ and the pressure p such that

$$\begin{aligned} {\text {div}}\mathbf{u} = 0,\quad -\Delta \mathbf{u} +\nabla p = \mathbf{f} \text { in }\varOmega , \end{aligned}$$

(38)

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

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

$$\begin{aligned} \mathbf{u} (x,y) = \cos \left( \frac{\pi }{2}(x^2+y^2)\right) \begin{pmatrix}y \\ -x \end{pmatrix}, \end{aligned}$$

where the corresponding forcing term reads

$$\begin{aligned} \mathbf{f} (x,y) =\pi \cos \left( \frac{\pi }{2}(x^2+y^2)\right) \begin{pmatrix} yr^2 \pi +4(y-x)\tan \left( \frac{\pi }{2}(x^2+y^2)\right) \\ -xr^2 \pi - 4(x+y)\tan \left( \frac{\pi }{2}(x^2+y^2)\right) \end{pmatrix}. \end{aligned}$$

In Fig. 9 we see the resulting $L^2$- and $H^1$-errors. Again we observe that $\varUpsilon$ 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.

Fig. 9

$L^2$- and $H^1$-errors w.r.t. mesh-size $h_{max}$ for varying parameter $\varUpsilon$ computed for the Stokes problem in two dimensions with linear finite elements

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 $\varUpsilon$ between the real and the approximated domains, cf. Theorem 2. We have illustrated the sharp nature of this bound in the computations performed in Sect. 4.

Particularly, in the case of first and second order approximation we observe how the relation between the mesh size h and aforementioned $\varUpsilon$ 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 [22]. 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  [30] and the difficulty of nonlinearities introduced by the convective term, and thus the non-uniqueness of solutions  [20].