Multilevel quadrature for elliptic problems on random domains by the coupling of FEM and BEM

Elliptic boundary value problems which are posed on a random domain can be mapped to a fixed, nominal domain. The randomness is thus transferred to the diffusion matrix and the loading. While this domain mapping method is quite efficient for theory and practice, since only a single domain discretisation is needed, it also requires the knowledge of the domain mapping. However, in certain applications, the random domain is only described by its random boundary, while the quantity of interest is defined on a fixed, deterministic subdomain. In this setting, it thus becomes necessary to compute a random domain mapping on the whole domain, such that the domain mapping is the identity on the fixed subdomain and maps the boundary of the chosen fixed, nominal domain on to the random boundary. To overcome the necessity of computing such a mapping, we therefore couple the finite element method on the fixed subdomain with the boundary element method on the random boundary. We verify the required regularity of the solution with respect to the random domain mapping for the use of multilevel quadrature, derive the coupling formulation, and show by numerical results that the approach is feasible.


Introduction
Many practical problems in science and engineering lead to elliptic boundary value problems for an unknown function. Their numerical treatment by e.g. finite difference or finite element methods is in general well understood provided that the input parameters are given exactly. This, however, is often not the case in practical applications.
If a statistical description of the input data is available, one can mathematically describe data and solutions as random fields and aim at the computation of corresponding deterministic statistics of the unknown random solution. The present article is dedicated to the treatment of uncertainties in the description of the computational domain. Applications are, besides traditional engineering, for example uncertain domains which are derived from inverse methods such as tomography. In recent years, this situation has become of growing interest, see e.g. [5,6,22,24,27,29,30] and the references therein.
In this article, we are first going to focus on the so-called domain mapping method, which has been introduced in [30] and rigorously analysed in [6,22], where analytic dependency of the solution on the random domain mapping with regard to the energy norm has been verified. Given enough spatial regularity of the random domain mapping, we first prove that the solution is analytically dependent on the random domain mapping also in the H s (D)-norm. The key idea of the method is to map the boundary value problem While the random domain mapping approach is mathematically natural, it is not neccessarily the setting, which is directly encountered in practical applications. This mainly stems from the fact that the random domain mapping does not only describe the random domains themselves but also includes a specific point correspondence between the domain realisations. In applications often only a description of the random boundary might be known, however in such cases the quantity of interest (5) QoI is generally sought on a deterministic subdomain, B, which almost surely is a subset of the domain realisations. Therefore, it is then necessary to be able to transform the description of the random domains given by a description of the random boundary and the specification of the subdomain into the form of a random domain mapping. In [30], the authors consider by using the vector-valued Laplace equation to compute such a random domain mapping. If more structure is given, for example when the random domains are described by star-shaped boundaries or more generally when they are directly given by a boundary mapping from a nominal boundary, one may also consider other approaches, such as transfinite interpolation techniques, see e.g. [12,13,14], to extend the mapping onto the whole reference domain. However, to overcome the necessity of computing such a random domain mapping in this setting, we propose to compute the quantity of interest by performing the calculations on the realisations of the random domains. Moreover, we can also sidestep the generation of a mesh on the random part of the domain D[ω] \ B, by coupling finite element methods with boundary element methods for the spatial approximation as follows: we apply finite elements on the subdomain B and treat the rest of the domain by a boundary element method. This is advantageous, since large domain deformations on coarse discretisations can be handled easily, as we do not need to mesh the random part of the domain but only its boundary. We present the resulting coupling formulation and then discuss the efficient solution by multilevel quadrature methods. Especially, since we verify the required regularity with respect to the random perturbation field, we also know that we have the required regularity on the deterministic subdomain B, under the assumption that there exists a transform from the random boundary description to the random domain mapping, which has sufficient regularity.
The rest of this article is organized as follows. Section 2 is dedicated to the mathematical formulation of the problem under consideration. The problem's regularity is studied in Section 3. Here, we provide estimates in stronger spatial norms which are needed for multilevel accelerated quadrature methods. The coupling of finite elements and boundary elements is the topic of Section 4. The multilevel quadrature method for the solution of the random boundary value problem is then introduced in Section 5. Numerical experiments are carried out in Section 6. Finally, we state concluding remarks in Section 7.

Notation and model problem
Before we complete the mathematical setting of our model problem, we will introduce the notations used throughout the rest of the article. Especially, for the regularity considerations in Section 3 some of the notation -and the choice of a certain weighting in the Sobolev-Bochner norms -helps keep formulas somewhat more concise and compact.
2.1. Notation and precursory remarks. We use N to denote the natural numbers including 0 and N * when excluding 0.
For a sequence of natural numbers, α = {α n } n∈N * ∈ N N * , we define the support of the sequence as supp α = {n ∈ N * | α n = 0} and say that α is finitely supported, if supp α is of finite cardinality, Then, N N * f denotes the set of finitely supported sequences of natural numbers and we refer to its elements as multi-indices. Furthermore, for all m ∈ N * we will identify the elements α = (α 1 , . . . , α m ) ∈ N m with their extension by zero into N N * f , that is α = (α 1 , . . . , α m , 0, . . .). Thus, by this identification, all notations defined for elements of N N * f also carry over to the elements of N m and we also refer to elements of N m as multi-indices. For and a sequence of real numbers γ = {γ n } n∈N * ∈ R N * , we use the following common notations: Furthermore, we say that α ≤ β holds, when α j ≤ β j holds for all j ∈ supp α ∪ supp β, and α < β, when α ≤ β and α = β hold. Subsequently, we will always equip R m with the norm · 2 induced by the canonical inner product ·, · and R m×m with the induced norm · 2 . Moreover, when considering R m itself or an open domain D ⊂ R m as a measure space we always equip it with the Lebesgue measure. Similarly, we always equip N and N * with the counting measure, when considering them as measure spaces.
Let X , X 1 , . . . , X r and Y be Banach spaces, then we denote the Banach space of bounded, linear maps from X to Y as B(X ; Y); furthermore, we recursively define B(X 1 , . . . , X r ; Y) := B X 1 ; B(X 2 , . . . , X r ; Y) and the special case For T ∈ B(X 1 , . . . , X r ; Y) and v j ∈ X j we use the shorthand notation Tv 1 · · · v r : For a given Banach space X and a complete measure space M with measure µ the space L p µ (M; X ) for 1 ≤ p ≤ ∞ denotes the Bochner space, see [26], which contains all equivalence classes of strongly measurable functions v : A function v : M → X is strongly measurable if there exists a sequence of countably-valued measurable functions v n : M → X , such that for almost every m ∈ M we have lim n→∞ v n (m) = v(m). Note that, for finite measures µ, we also have the usual inclusion L p µ (M; X ) ⊃ L q µ (M; X ) for 1 ≤ p < q ≤ ∞.
For a given Banach space X and an open domain D ⊂ R d , with d ∈ N * , the space W η,p (D; X ) for η ∈ N and 1 ≤ p ≤ ∞ denotes the Sobolev-Bochner space, which contains all equivalence classes of strongly measurable functions v : D → X , such that the function itself and all weak derivatives up to total order η are in L p (D; X ) with the norm v η,p,D;X := v W η,p (D;X ) := Moreover, W η,p 0 (D; X ) denotes the closure of the linear subspace of smooth functions with compact support, C ∞ c (D; X ), in W η,p (D; X ) and we set H η (D; X ) := W η,2 (D; X ) and H η 0 (D; X ) := W η,2 0 (D; X ). As usual, we use C ω (D; X ) to denote the real analytic functions from D to X and C k,s (D; X ) to denote the Hölder spaces. For a bi-Lipschitz function v : D → X we denote its bi-Lipschitz constants by In the notation for the Bochner, Sobolev-Bochner and Hölder spaces, we may omit specifying the Banach space X when X = R. Especially, H −η (D) denotes the topological dual space of H η 0 (D). Moreover, if the X we are considering is itself a Bochner or Sobolev-Bochner space, then we replace the X in the subscript of the norm with the subscripts of its norm, for example . Lastly, to avoid the use of generic but unspecified constants in certain formulas, we use c d to mean that c can be bounded by a multiple of d, independently of parameters which c and d may depend on. Obviously, c d is defined as d c and we write c d if c d and c d.

2.2.
Model problem. Let τ ∈ N and d ∈ N * ; D ⊂ R d denote the reference domain with boundary ∂D that is of class C τ,1 -when τ = 1 then we also consider the case where D is a bounded and convex domain with Lipschitz continuous boundary -and (Ω, F, P) be a separable, complete probability space with σ-field F ⊂ 2 Ω and probability measure P. Furthermore, let V ∈ L ∞ P Ω; C τ,1 (D; R d ) be the random domain mapping. Moreover, we require that, for P-almost any ω, is bi-Lipschitz and fulfils the uniformity condition Finally, we require that the we have a hold-all domain D that satisfies D[ω] ⊂ D for P-almost any ω ∈ Ω and consider f ∈ C ω (D).
Note that while we restrict ourselves to the Poisson equation here to simplify the analysis, the extension of the regularity result to an operator div x A ∇ x , with an A ∈ C ω (D; R d×d ) and A fulfilling an ellipticity condition is straightforward.
While, by definition, we know that for Palmost any ω ∈ Ω, we also have the following stronger result.
Proof. The fact that V[ω] is a C τ -diffeomorphism follows directly from the inverse funtion theorem. Then, with the explicit formula for the τ -th derivative of V[ω] −1 from the inverse funtion theorem, one can bound we can use the one-to-one correspondence to pull back the model problem onto the reference domain D instead of considering it on the actual domain realisations D[ω]. According to the chain rule, we then have Now, with (3) this leads us to the following formulation of our model problem (2) on the reference domain, cf. [22]: for P-almost every ω ∈ Ω and allv ∈ H 1 0 (D). Note, especially, that by the uniformity condition we have that Without loss of generality, we assume σ ≤ 1 ≤ σ.
From here on, we assume that the spatial variable x and the stochastic parameter ω of the random field have been separated by the Karhunen-Loève expansion of V coming from the mean field E[V] and the covariance Cov[V] yielding a parametrised expansion is a sequence of uncorrelated random variables, see e.g. [22]; we denote the pushforward measure of P onto as P y . Thus, we then also view all randomness as being parametrised by y, i.e. ω, Ω and P are replaced by y, and P y .
We now impose some common assumptions, which make the Karhunen-Loève expansion computationally feasible.

Regularity
To prove the analyticity of the mappingû : → H τ +1 (D), we first investigate the analyticity of the mappingsÂ : in a first subsection. Based on that analyticity we then prove the analyticity forû in the second subsection.
To make the notation less cumbersome, since we are considering the norm of spaces of the form L ∞ Py ( ; X ), we introduce the shorthand notatioñ v~X := v ∞, ;X .
We will especially make use it for spaces of the form L ∞
3.1. Parametric regularity of the diffusion coefficient and the right-hand side. To provide regularity estimates for the diffusion coefficientÂ and the right hand sidef , that are based on the decay of the expansion of V as per Assumption 2.2, we first note that we can write Therefore, we first discuss the regularity of the combined mapping , for which we have the following result.
Here, c τ denotes the constant coming from the embedding Proof. By definition we have that J[y] = D x V[y] and so it follows that From this we can derive that first order derivatives are given by and all higher derivatives vanish. Clearly, this affine dependence on y implies the bounds.
Next, we supply bounds on the derivatives of the mappings T and s.
We start with the mappings Then, using [23, Lemma 3] we see, that the mapping Next, we consider the mapping Clearly, the r-th Fréchet derivative of inv is given by where S r is the set of all bijections on the set {1, 2, · · · , r}. Thus, we have Therefore, we can use [23,Lemma 4] to see, that the mapping where M [i1,H1],..., [it,Ht] denotes the matrix M whose i k -th column is replaced by the i k -th column of the matrix H k for all k from 1 to r. Now, since we can bound the determinant of a matrix by the product of the norms of its columns, i.e.
and since we know that z j ≤ z 1 · · · z d .
it follows that, with k det = (2σ) d and c det = 1. As before, we can use [23,Lemma 4] to see, that the mapping Proof. We start with the mapping which is infinitely Fréchet differentiable with Moreover, as shown in the previous proof we also have that Now, these results enable us to show the following regularity estimates for the diffusion coefficientÂ and the right hand sidef .
as well as, for α = 0, where C(α, s) is the set of all compositions of the multi-index α into s non-vanishing multiindices β 1 , . . . , β s , see [23], and we make use of the combinatorial identity This proves the assertion forÂ, while the assertion forf follows analogously after remarking that

Parametric regularity of the solution.
For this subsection, we require an elliptic regularity result, which we state as an assumption: Assumption 3.5. Let D be a sufficiently smooth such that, for all B ∈ C τ −1,1 (D; R d×d symm ) that fulfil (7), we have that the problem of solving for any h ∈ H τ −1 (D) has a unique solution u ∈ H 1 0 (D), which also lies in H τ +1 (D), with u τ +1,2,D ≤ C er h τ −1,2,D , where C er only depends on D, σ, σ, τ and a bound on B C τ −1,1 (D;R d×d symm ) . Such an elliptic regularity estimate for example is known for τ = 1, when the domain is convex and bounded, see [16, Propositions 3.2.1.2 and 3.1.3.1]. The elliptic regularity estimate is also known to hold for τ ≥ 1 and d = 2, when the domain's boundary is smooth, see [4].
This obviously directly implies the following result.
Proof. By differentiation of the variational formulation (6) with respect to y we arrive, for .
Applying the Leibniz rule on the left-hand side yields .
Then, by rearranging and using the linearity of the gradient, we find .
Using Green's identity, we can then write .

The coupling of FEM and BEM
While we have considered general random domain mappings in the previous subsections, we will now restrict them according to the remarks made in the introduction. That is, we assume for the rest of the article that we are given a random boundary description, Γ[ω], and a fixed, deterministic subdomain B, which describe our random domain, compare Figure 1 when Γ = Γ[ω]. Moreover, we assume that we are interested the some quantity of interest that is based on the knowledge of u| B as in (5).
We will assume that there is a random domain mapping V which fulfils the Assumption 2.2 as well as fulfilling V[ω]| B = Id B and V[ω](∂D) = Γ[ω] for almost any ω. Then, we know from the previous section thatû : → H τ +1 (D) is analytic which also implies that u| B : → H τ +1 (B) is analytic.
So, to be able to use multilevel quadrature to compute the quantity of interest efficiently, we consider a formulation here, that enables us to compute the Galerkin solution u h [ω] ∈ H 1 (B) with a mesh on B but without needing a mesh on D[ω] \ B or needing the knowledge of the random domain mapping. One arrives at such a formulation by reformulate the boundary value problem as two coupled problems involving only boundary integral equations on the random boundary Γ[ω].
4.1. Newton potential. For sake of simplicity in representation, we shall restrict ourselves the deterministic boundary value problem (12) − ∆ u = f in D, u = 0 on Γ := ∂D, i.e., the domain D is assumed to be fixed. Of course, when applying a sampling method for (1), the underlying domains are always different. In order to resolve the inhomogeneity in (12), we introduce a Newton potential N f which satisfies Here, D is a sufficiently large domain containing D[ω] almost surely. The Newton potential is supposed to be explicitly known like in our numerical example (see Section 6) or computed with sufficiently high accuracy. Especially, since the domain D can be chosen fairly simple, one can apply finite elements based on tensor products of higher order spline functions (in [−R, R] d ) or dual reciprocity methods. Notice that the Newton potential has to be computed only once in advance.
By making the ansatz (14) u = N f +ũ and settingg := g − N f , we arrive at the problem of seeking a harmonic functionũ which solves the following Dirichlet problem for the Laplacian Now, we are able to apply the coupling of finite elements and boundary elements.

4.2.
Reformulation as a coupled problem. For the subdomain B ⊂ D, we set Σ := ∂B, see Figure 1 for an illustration. The normal vectors n at Γ and Σ are assumed to point into D \ B. We shall split (15) in two coupled boundary value problems in accordance with In order to derive suitable boundary integral equations for the problem in D \ B, we define the single layer operator V ΦΨ , the double layer operator K ΦΨ and its adjoint K ΨΦ , and the hypersingular operator W ΦΨ with respect to the boundaries Φ, Ψ ∈ {Γ, Σ} by Figure 1. The domain D, the subdomain B, and the boundaries Γ = ∂D and Σ = ∂B.

4.3.
Variational formulation. We next introduce the product space H := H 1 (B)×H −1/2 (Σ)× H −1/2 (Γ), equipped by the product norm Further, let a : H × H → R, be the bilinear form defined by For sake of simplicity in representation, we omitted the trace operator in expressions like (w, W ΣΣ v) L 2 (Σ) etc.
Introducing the linear functional F : H → R, , the variational formulation is given by: for all (w, λ Σ , λ Γ ) ∈ H. In accordance with [11], the variational formulation (19) admits a unique solution (ũ, σ Σ , σ Γ ) ∈ H for all F ∈ H , provided that D has a conformal radius which is smaller than one if d = 2.

Galerkin discretization.
Since the variational formulation is stable without further restrictions, the discretization is along the lines of [19]. We first introduce a uniform triangulation of B which in turn induces a uniform triangulation of Σ. Moreover, we introduce a uniform triangulation of the free boundary Γ, which we suppose to have the same mesh size as the triangulation of the domain B. For the FEM part, we consider continuous, piecewise linear ansatz functions {ϕ B k : k ∈ ∆ B } with respect to the given domain mesh. For the BEM part, we employ piecewise constant ansatz functions {ψ Φ k : k ∈ ∇ Φ } on the respective triangulations of the boundaries Φ ∈ {Σ, Γ}.
For sake of simplicity in representation, we set ϕ Σ k := ϕ B k | Σ for all k ∈ ∆ B . Note that most of these functions vanish except for those with nonzero trace which coincide with continuous, piecewise linear ansatz functions on Σ. Finally, we shall introduce the set of continuous, piecewise linear ansatz functions on the triangulation of Γ, which we denote by Then, introducing the system matrices where again Φ, Ψ ∈ {Σ, Γ}, and the data vector we obtain the following linear system of equations We mention that G −1 Γ g corresponds to the L 2 (Γ)-orthogonal projection of the given Dirichlet datag ∈ H 1/2 (Γ) onto the space of the continuous, piecewise linear ansatz functions on Γ. That way, we can also apply fast boundary element techniques to the boundary integral operators on the right hand side of the system (21) of linear equations.

4.5.
Multilevel based solution of the coupling formulation. We shall encounter some issues on the efficient multilevel based solution of the system (21) of linear equations. The complexity is governed by the BEM part since the boundary element matrices are densely populated. Following [19,20], we apply wavelet matrix compression to reduce this complexity such that the over-all complexity is governed by the FEM part. On the other hand, according to [20,25], the Bramble-Pasciak-CG (see [2]) provides an efficient and robust iterative solver for the above saddle point system. Combining a nested iteration with the BPX preconditioner (see [3]) for the FEM part and a wavelet preconditioning (see [9,28]) for the BEM part, we derive an asymptotical optimal solver for the above system, see [20] for the details. We refer the reader to [20] for the details of the implementation of a similar coupling formulation.

Multilevel quadrature method
The crucial idea of the multilevel quadrature to compute the quantity of interest (5) is to combine an appropriate sequence of quadrature rules for the stochastic variable with the multilevel discretization in the spatial variable. To that end, we first parametrize the quantity of interest by using (8) over the cube = [−1, 1] N * and compute (22) QoI where u −1 := 0. Herein, for the spatial approximation, we shall use the multilevel representation from Subsection 4.5 to compute the Galerkin solution u ∈ H 1 (B) on level that corresponds to the step size h = 2 − . For the approximation in the stochastic variable y, we shall thus provide a sequence of quadrature formulae {Q } for the integral v(x, y) dP y of the form For our purposes, we assume that the number of points N of the quadrature formula Q is chosen such that the corresponding accuracy is (23) ε = 2 − , = 0, 1, . . . , L.
Since the multilevel quadrature can be interpreted as a sparse-grid approximation, cf. [21], it is known that mixed regularity results of the integrand have to be provided as derived in Section 3, compare [10,15,21] for example. Since the mapping u : → H τ +1 (B) is analytic, we can especially apply the quasi-Monte Carlo method, the Gaussian quadrature, or the sparse grid quadrature, see e.g. [15]. Especially, in case of H 2 -regularity (τ = 1) and F = u| B , i.e., QoI(u) = E(u| B ), we obtain then the error estimate Notice that the computational complexity of the multilevel quadrature (22) is considerably reduced compared to a standard single-level quadrature method which has the same accuracy, see e.g. [1,7,21].

Numerical results
In our numerical example, we consider the reference domain D to be the ellipse with semiaxis 0.6 and 0.4. We represent its boundary by γ ref : [0, 2π) → ∂D in polar coordinates and perturb this parametrization in accordance with where y k ∈ (−0.5, 0.5) for all k ∈ Z and ε = 0.05. The weights w k are chosen as w k = 1 for all |k| ≤ 5 and w k = (k − 5) −6 for all |k| > 5. Hence, we have the decay γ k ∼ k −4 for the choice τ = 1, which is sufficient for applying the quasi Monte Carlo method based on the Halton sequence. In practice, we set all w k to zero if |k| > 64 which corresponds to a dimension truncation after 129 dimensions. The random parametrization γ[ω] induces the random domain D[ω]. The fixed subset B ⊂ D is given as the ball of radius 0.2, centered in the origin. For an illustration of four random draws, see  1  37  32  2  129  64  3  481  128  4  1857  256  5  7297  512  6  28929  1024  7  115201  2048  8 459777 4096 Table 1. Number of the degrees of freedom of the finite element method and the boundary element method.  be given. A suitable Newton potential is analytically given by N f = −(x 2 1 + x 2 2 )/4. We consider the L 2 -tracking type functional as quantity of interest, where u is a given function. The coarse triangulation of B, based on Zlámal's curved finite elements [31], consists of 14 curved triangles on the coarse grid, which are then uniformly refined to get the triangulation on the finer grids. The 14 triangles correspond to eight piecewise linear and constant boundary elements each on the boundary ∂B. At the boundary ∂D, we likewise consider eight piecewise linear and constant boundary elements each on level 0. When applying uniform refinement, we arrive at the numbers of degrees of freedom in the finite and boundary element spaces found in Table 1.
In order to compute the quantity of interest, we will employ the quasi-Monte Carlo method based on the Halton sequence, see [17] for example. Since the exact solution is unknown, we compute first the quantity of interest on the spatial discretization level 8 by using 100 000 Halton points. Next, we compute the solution by the multilevel quasi-Monte Carlo method. Namely, for the multilevel quasi-Monte Carlo method on level L, we apply N = 2 L− N L and N = 4 L− N L Halton points, respectively, on the coarser levels 0 ≤ ≤ L, where we choose N L = 10, 20, 40 fine grid samples.
As it is seen in Figure 3, we always observe the same linear and quadratic convergence rate, respectively, but with different constants involved. Notice that linear convergence is in accordance with (24) while quadratic convergence is in accordance with (25). Only for the level = 8 and desired quadratic convergence, we observe a stagnation of the convergence. This issues from the fact that the reference solution computed by a single-level quadrature method is not accurate enough.

Conclusion
We provided regularity estimates of the solution to elliptic problems on random domains which allow for the application of multilevel quadrature methods. In order to avoid the need to compute either a random domain mapping or to generate meshes for every domain sample, we couple finite elements with boundary elements. It has been shown by numerical experiments that this approach is indeed able to exploit the additional regularity we have in the underlying problem without causing numerical problems on too coarse grids.