Finite Element Error Estimates on Geometrically Perturbed Domains

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 H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H^1$$\end{document}- and L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document}-error estimates for the Laplace problem. Theoretical considerations are validated by a computational example.


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 B Piotr Minakowski piotr.minakowski@ovgu.de Thomas Richter thomas.richter@ovgu.de 1 Institute of Analysis and Numerics, Otto-von-Guericke University Magdeburg, Universitätsplatz 2, 39106 Magdeburg, Germany 2 Interdisciplinary Center for Scientific Computing, Heidelberg University, INF 205, 69120 Heidelberg, Germany 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.

Notation
Let Ω ⊂ R d be a domain with dimension d ∈ {2, 3}. By L 2 (Ω) we denote the Lebesgue space of square integrable functions equipped with the norm · Ω . By H 1 (Ω) we denote the space of L 2 (Ω) functions with first weak derivative in L 2 (Ω) and by H m (Ω) for m ∈ N 0 we denote the corresponding generalizations with weak derivatives up to degree m ∈ N 0 . The norms in H m (Ω) are denoted by · H m (Ω) . For convenience we use the notation H 0 (Ω) := L 2 (Ω). By H 1 0 (Ω) we denote the space of those H 1 (Ω) functions that have vanishing trace on the domain's boundary ∂Ω and we use the notation H 1 0 (Ω; Γ ) if the trace only vanishes on a part of the boundary, Γ ⊂ ∂Ω. Further, by (·, ·) Ω we denote the L 2 (Ω)-scalar product and ·, · Γ the L 2 -scalar product on a d − 1 dimensional manifold Γ , e.g. Γ = ∂Ω. Moreover, [∂ n ψ] is the jump of the normal derivative of ψ, i.e. for x ∈ Γ with normal n (that is normal to

Laplace Equation and Domain Perturbation
On Ω ⊂ R d let f ∈ L 2 (Ω) be the given right hand side. We consider the Laplace problem with homogeneous Dirichlet boundary conditions, The variational formulation of this problem is given by: The boundary ∂Ω is supposed to have a parametrization in C m+2 , where m ∈ N 0 . Given the additional regularity f ∈ H m (Ω), H 0 (Ω) := L 2 (Ω), there exists a unique solution satisfying the a-priori estimate see e.g. [16]. In the following we assume that the real domain Ω is not exactly known but only given up to an uncertainty. We hence define a second domain, the reconstructed domain Ω r with a boundary that allows for C m+2 parametrization. The Hausdorff distance between both domains is then denoted by Υ ∈ R, This distance Υ is not necessarily small. When it comes to spatial discretization we will be interested in both cases, h Υ as well as Υ 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 Ω and Ω r . Assumption 1 (Domains) Let Ω and Ω r be two domains with Ω∩Ω r = ∅ and with Hausdorff distance Υ ∈ R. Both boundaries allow for a local C m+2 parametrization, m ∈ N 0 . Let We assume that there exists a cover of S by a finite number of open rectangles (or rectangular cuboids) {R 1 , . . . , R n(S) }. Each rectangle R i is given as translation and rotation of (0, 1 Υ with a constant c S,1 ≥ 0. Following conditions hold: 1.a) On each rectangle R, the boundary lines ∂Ω ∩ R and ∂Ω r ∩ R allow for unique parametrizations g R Ω (t) and g R Ω 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.
where c S,2 > 0 is a constant.
For the following we set c S := max{c S,1 , c S,2 }.  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 ∂Ω r and ∂Ω. 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 Ω r we must face the technical difficulty that the right hand side f ∈ H m (Ω) is not necessarily defined on Ω r . We therefore weaken the assumptions on the right hand side.
In addition we assume that the right hand side on Ω r can be bounded by the right hand side on Ω, i.e.
An alternative would be to use Sobolev extension theorems to extend functions f ∈ H m (Ω) from Ω to Ω r , see [8].
On Ω r we define the solution u r ∈ H 1 0 (Ω r ) to the perturbed Laplace problem The unique solution to (5) satisfies the bound Remark 1 (Extension of the solutions) A difficulty for deriving error estimates is that u is defined on Ω and u r on Ω r = Ω. Since the domains do not match, u may not be defined on all of Ω 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 Ω r \ Ω and u r := 0 on Ω \ Ω r . Globally, both functions still have the regularity u, u r ∈ H 1 (Ω ∪ Ω r ). We will use the same notation for discrete functions u h ∈ V h defined on a mesh Ω h and extend them by zero to 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.
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 as, if the line would leave this remainder, it would cut each line more than once which opposes Assumption 1b). Integrating the function ψ along this line gives Applying Hölder's inequality to the second term on the right hand side, with the length of the line bounded by c S γ , we obtain Using the parametrizations The volume integral on the right hand side is exactly the integral over R ∩ (V \ W ). The boundary integrals can be interpreted as path integrals and therefore be estimated by As the boundaries allow for a C 2 parametrization, we estimate Summation over all rectangles and estimation of all overlaps by means of Assumption 1 gives For ψ ∈ H 1 0 (V ), the term on ∂ V vanishes. To show the second estimate on V \ W we again pick one rectangle R and consider a point By the same arguments as above it holds We integrate over s and t to obtain 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 Ω and Ω r , specifically to the case of use.
We continue by estimating the difference between the solutions of the Laplace equations on Ω and on Ω r . (2) and (5) respectively, it holds that We separate the domains of integration and integrate by parts Combining the boundary terms on the right-hand side of (11), into an integral over ∂Ω and a jump term over ∂Ω r ∩ Ω, we obtain In Finally, u r = 0 on Ω \ Ω r , such that the boundary terms reduce to Combining (12)- (14) and using the Cauchy-Schwarz inequality, we estimate Since u, u r ∈ H 2 (Ω ∩ Ω r ), the trace inequality gives Applying Lemma 3 twice: to ψ = u and to ψ = ∇u (same for u r ), and extending the norms from Ω \ Ω r to Ω and from Ω r \ Ω to Ω r give the bounds With the trace inequality and the a priori estimates u H 2 (Ω) ≤ c f Ω and u r H 2 (Ω r ) ≤ c f r Ω r ≤ c f Ω we obtain the bounds Using the fact that u = 0 on ∂Ω we apply (7) twice and use the trace inequality to get the estimate We can then estimate f Ω\Ω r ≤ f Ω by extending to the complete domain. Combining (16) with (18) and (19) we obtain the estimate which concludes the H 1 -norm bound.
(ii) For the L 2 -estimate we introduce the adjoint problem which allows for a unique solution satisfying the a-priori bound z H 2 (Ω) ≤ c s with the stability constant c s < ∞. Testing with u − u r and integrating by parts twice gives It holds z = 0 and u = 0 on ∂Ω, The boundary terms z ∂Ω r ∩Ω and u r ∂Ω are estimated with Lemma 3, the normal derivatives by the trace inequality and the terms on Ω \ Ω r by (7) u − u r Ω ≤ cΥ The L 2 -norm estimate follows by using the bounds u H 2 (Ω) ≤ c f Ω , u r H 2 (Ω r ) ≤ c f Ω and z H 2 (Ω) ≤ c.

Remark 2
The estimate f Ω\Ω r ≤ c f Ω is not optimal. Further powers of Υ are easily generated at the cost of a higher right hand side regularity. Also, the estimate ∂ n (u − u r ) ≤ c( u H 2 (Ω) + u r H 2 (Ω r ) ) by Cauchy Schwarz and the trace inequality could be enhanced to produce powers of Υ . The limiting term in (12) however is the boundary integral | ∂ n u r , u ∂Ω r ∩Ω | = O(Υ 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 Ω ∩ Ω r only and that allows us to improve the order to O(Υ ) in the H 1 -case by avoiding exactly this boundary integral.

Discretization
The starting point of a finite element discretization is the mesh of the domain Ω. In our setting we do not mesh Ω directly, because the domain Ω is not exactly known. Instead, we consider a mesh of the reconstructed domain Ω r .
We partition Ω r into a parametric triangulation Ω h , consisting of open elements T ⊂ R d . Each element T ∈ Ω h stems from a unique reference elementT 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 :T → T is a polynomial of degree r ∈ N. We will consider iso-parametric finite element spaces, that are based on polynomials of the same degree r ≥ 1. In the following we assume structural and shape regularity of the mesh such that standard interpolation estimates will hold for all elements T ∈ Ω 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 elementT letP be a polynomial space of degree r , e.g.
on quadrilateral and hexahedral meshes. Then, the finite element space V r h on the mesh Ω h is defined as This parametric finite element space does not exactly match the domain Ω r . Given an isoparametric mapping of degree r it holds dist(∂Ω r , ∂Ω h ) = 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 ∈ N 0 and let Ω r be a domain with a boundary that allows for a parametrization of degree m + 2. Let f r ∈ H m (Ω r ) and u h ∈ V r h ∩ H 1 0 (Ω h ) be the isoparametric finite element discretization of degree 1 ≤ r ≤ m + 1 It holds We formulated the error estimate on the domain Ω r although the finite element functions are given on Ω 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.
as well as Proof (i) We start with the H 1 error. Inserting ±u r and extending the finite element error u r − u h from Ω to Ω r , where a small remainder appears, we have 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 Ω \ Ω r , we obtain We continue with the remainder ∇u h on Ω \ Ω r , which is non-zero on Ω h only . This remaining stripe has the width and we apply Lemma 3 to get where the second derivative ∇ 2 u h is understood element wise. This term is extended to Ω h and with the inverse estimate and the a priori estimate for the discrete solution we obtain with To the first term on the right hand side of (23) we add ±u r and ±I h u r , the nodal interpolation of u r into the finite element space Here, the first and last terms are estimated with the trace inequalities and, in the case of the discrete term with the inverse inequality 1 , followed by adding ±u r we get We used both γ h, Then, collecting all terms in (22)- (26) and using the interpolation estimates as well as Theorem 1 we finally get which shows the a priori estimate since 3r − 1 ≤ 2r for all r ≥ 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 ∂Ω h , where u h = 0, can be applied, i.e.
The L 2 -estimate directly follows with Lemma 4, Theorem 1 and by the a priori estimate

Remark 3 (Polygonal domains)
In two dimensions, the extension of the error estimates to the case of convex polygonal domains, where u ∈ H 2 (Ω) and u r ∈ H 2 (Ω r ), is relatively straightforward. In this case, Ω h fits Ω 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 ∇(u r −u h ) ≤ c f . The extension of Lemma 3, which locally requires smoothness of the parametrizations g R W (·) and g R V (·), 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 ∂Ω and Ω r are cut by the edges of rectangles. This allows to derive the optimal error estimates . 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 O(Υ 1 2 ) and O(Υ ), that describes the discrepancy between Ω and Ω 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 − u = 4 on the unit disc Ω = B 1 (0) and − u r = 4 on the shifted domain Ω r = B 1 (Υ ). The errors in H 1 norm and L 2 norms expressed on the complete domain Ω are estimated by A closer analysis shows that the main error -in the H 1 -case -occurs on the small shaded stripe Ω \ Ω r such that

The usual Galerkin error u r
) of iso-parametric finite element approximations contributes to the overall error. For Ω = Ω r , i.e. Υ = 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.
In Sect. 4 we discuss the difficulty of measuring errors on an unknown domain Ω. 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 Ω \ Ω r , where no finite element mesh is given. These remainders contribute the lowest order parts Υ 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 Υ in the H 1 -norm error. In addition to the previous setting we require a regular map T r : Ω → Ω r between the two domains. By pulling back Ω r to Ω via this map a Jacobian arises that controls the geometrical error and that hence has to be controllable by Υ .

Corollary 1 In addition to the assumptions of Theorem
Further, let the following regularity of problem data hold in addition to Assumption 2 and let the solution satisfy Then, it holds Proof We start by splitting the error into domain approximation and finite element approximation errors An optimal order estimate of the finite element error is given in Theorem 1.
To estimate the first term of the right hand side of (31) we introduce the function which satisfiesû r ∈ H 1 0 (Ω) and solves the problem We introduce the notation e r := u −û r , extend the first term from Ω ∩ Ω r to Ω and insert where we also used Poincaré's estimate. For bounding f −f r we consider a point x ∈ Ω ∩Ω r , use the higher regularity of the right hand side (29) to estimate by a Taylor expansion where ξ ∈ Ω is some point on the line from x to T r (x). We take the square and integrate over Ω to get the estimate where Ω Υ is a enlargement of Ω by at most O(Υ ), since intermediate values ξ used in (35) are not necessarily part of Ω ∪ Ω r . This argument is also applicable to the second term on the right hand side of (33) such that it holds 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 : Ω → Ω r depends on the given application. Here a construction, corresponding to the ALE map, can be realised by means of a domain deformationd : Such a construction is common in fluid-structure interactions, see [31,Section 2.5.2]. Given that |d|, |∇d| = O(Υ ) it holds While the assumption |d| = O(Υ ) is easy to satisfy since dist(∂Ω, ∂Ω r ) ≤ Υ , the condition |∇d| = O(Υ ) 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 Ω be the unit sphere, Ω r be an ellipse It holds dist(∂Ω, ∂Ω r ) ≤ Υ and we define the map T r : Ω → Ω r by This map satisfies the assumptions of the corollary

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 Υ . Moreover, we numerically extend the analytical predictions and show that a similar behavior holds for the Stokes system. We consider Ω to be a unit ball in two and three dimensions and define a family of perturbed domains Ω Υ , with the amplitude of the perturbation being dependent on the coefficient Υ , cf. Fig. 5.
In two dimensions, the boundary of the domain Ω Υ is given in polar coordinates (ρ, ϕ) by 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 Ω Υ , with a mesh size h, by Ω h,Υ .
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].

Laplace Equation in Two and Three Dimensions
We consider the following problem where Ω 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 with r = x 2 + y 2 in two and r = x 2 + y 2 + z 2 in three dimensions, respectively, which results in the right hand sides For the ease of evaluations the errors, the H 1 -and L 2 -norms will be computed on the truncated domains Ω 2d = {(ϕ, ρ) for ϕ ∈ [0, 2π) and ρ ∈ (0, 0.88)}, Ω 3d = {(ϕ, θ, ρ) for θ ∈ [0, π), ϕ ∈ [0, 2π) and ρ ∈ (0, 0.88)}, see also Remark 4. We hence do not compute the errors ∇(u − u h ) and u − u h on the remainders Ω \ Ω r . Therefore we expect optimal order convergence in the spirit of Corollary 1. The restriction of the domain to an area within Ω 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, Υ 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 Υ , 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 Υ . 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.

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 div u = 0, − u + ∇ p = f in Ω,   (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 where the corresponding forcing term reads f(x, y) = π cos π 2 (x 2 + y 2 ) yr 2 π + 4(y − x) tan π 2 (x 2 + y 2 ) −xr 2 π − 4(x + y) tan π 2 (x 2 + y 2 ) .
In Fig. 9 we see the resulting L 2 -and H 1 -errors. Again we observe that Υ 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.

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 Υ 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 Υ 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].