Faculty of Electrical Engineering and Computer Science Department of Applied Mathematics
Ph.D. Thesis
The Boundary Element Method for Shape Optimization in 3D
Author:
Ing. Jan Zapletal
Supervisor:
Doc. RNDr. Jiří Bouchala, Ph.D.
2017
or other qualification. Where other sources of information have been used, they have been acknowledged.
In Ostrava on January 18, 2017 . . . .
adventurous world of mathematics and the boundary element method and provided me with his guidance and valuable advice whenever needed. For numerous fruitful discussions I would like to thank doc. Ing. Dalibor Lukáš, Ph.D. It has always been pleasure to work in the friendly atmosphere of the Department of Applied Mathematics.
The topic of the thesis was formed during my internship at the Institut für Numerische Mathematik at TU Graz. My thanks thus also go to Univ.–Prof. Dr. Olaf Steinbach and Ass. Prof. Günther Of for their leadership during the work on shape optimization problems.
The stay was an enjoyable one thanks to Günther who shared his office with me and other members of the department including Arno Kimeswenger, Martin Neumüller, Elias Karabelas, Dominik Amann, Marco Zank, Gerhard Unger, Barbara Pöltl, and others.
I would like to thank prof. Ing. Tomáš Kozubek, Ph.D., the Head of the Research Pro- gramme 3 at the IT4Innovations National Supercomputing Center, for his support of the devel- opment of the boundary element library described in the thesis. Special thanks go to Michal Merta for his long-term friendship and fruitful cooperation on the BEM4I project and other activities at IT4Innovations and the Department of Applied Mathematics.
I am particularly grateful for the support given by my family not only throughout my studies.
Last but not least, I would like to thank Monika for making the life outside the office enjoyable.
has allowed it to be considered in the shape optimization context, where it is necessary to solve a given state problem many times. We present a BEM-based shape optimization concept, which can also be used for the solution of inverse problems including the presented Bernoulli free-surface problem. To separate the computational and optimization meshes we use the hierarchy of control meshes constructed by means of subdivision surfaces known from the computer graphics. In the thesis we also address the important topic of efficient implementation of BEM on modern hardware architectures and accelerators. The theory is supported by a series of numerical experiments validating the proposed approach.
Keywords: boundary integral equations, boundary element method, Bernoulli problem, shape optimization, subdivision surfaces, vectorization, manycore and multicore acceleration
1 Introduction 1
2 The Bernoulli free boundary problem 5
2.1 Function spaces . . . 5
2.1.1 Continuous functions . . . 5
2.1.2 Lebesgue and Sobolev spaces . . . 7
2.1.3 Continuous functions on manifolds . . . 8
2.1.4 Lebesgue and Sobolev spaces on manifolds . . . 10
2.2 Continuous optimization problem . . . 14
2.2.1 Existence of optimal domains . . . 18
2.2.2 Shape derivative of the cost . . . 29
2.3 Boundary integral equations . . . 37
2.3.1 Potentials and boundary integral operators . . . 38
2.3.2 Boundary integral equations . . . 41
3 Representation of smooth curves and surfaces 43 3.1 Bézier curves . . . 43
3.2 B-spline curves . . . 44
3.3 NURBS curves and surfaces . . . 46
3.4 Box splines . . . 47
3.5 Subdivision curves . . . 49
3.6 Loop’s subdivision . . . 54
4 Efficient implementation of BEM and shape optimization problems 59 4.1 Boundary element method . . . 59
4.1.1 Ansatz and test spaces . . . 59
4.1.2 Dirichlet problem . . . 61
4.1.3 Neumann problem . . . 62
4.2 Efficient implementation of BEM . . . 64
4.2.1 Numerical scheme . . . 67
4.2.2 Numerical line quadrature . . . 69
4.2.3 Vectorization of the numerical scheme . . . 70
4.2.4 Many core acceleration of the numerical scheme - offload mode . . . 79
4.2.5 Many core acceleration of the numerical scheme - native mode . . . 86
4.2.6 Semi-analytic scheme . . . 90
4.2.7 Analytic quadrature of collocation integrals . . . 91
4.2.8 Vectorized evaluation of the semi-analytic scheme . . . 94
4.3 Discrete optimization problem . . . 96
4.3.1 BEM4I interface for shape optimization problems . . . 102
4.4 Numerical experiments . . . 103
4.4.1 Ball . . . 103
4.4.2 Toy duck . . . 106
4.4.3 Industrial example . . . 109
A Acknowledgements 115
B Summary of publications 117
2.1 Local parametrization of∂Ω. . . . 8
2.2 Exterior Bernoulli problem configuration. . . . 14
2.3 Space-time cylinder. . . . 30
3.1 Bézier curve. . . . 44
3.2 B-spline. . . . 45
3.3 Refinement of a uniform cubic B-spline.. . . 46
3.4 Spline surface as a tensor product of one dimensional splines. . . . 47
3.5 NΞ based on the triangulation ofR2. . . . 48
3.6 The Chaikin corner cutting algorithm. . . . 49
3.7 Subdivison of the control polygon. . . . 50
3.8 1D subdivision and the invariant neighbourhood for the interval[−1,1]. . . . 51
3.9 Smoothing property of the operator SR=S(STS)−1ST. . . . 53
3.10 Triangle quadrisection. . . . 54
3.11 Loop subdivison rules.. . . 54
3.12 Multiresolution editing of a subdivision surface. . . . 57
4.1 Triangular mesh elements parametrized byR: ˆτ →τ. . . . 60
4.2 Discrete spaces on the boundary. . . . 60
4.3 Scalar and vectorized addition of two vectorsc:=a+b. . . . 65
4.4 Unaligned memory access, loop peeling.. . . 72
4.5 Array of structures vs. structure of arrays. . . . 75
4.6 OpenMP parallelized assembly of the BEM matrices. . . . 78
4.7 OpenMP vectorized assembly of the BEM matrices. . . . 79
4.8 Asymmetric offload programming model. . . . 80
4.9 Distribution of the matrix assembly. . . . 81
4.10 Accelerated assembly of the BEM matrices, double precision. . . . 83
4.11 Accelerated assembly of the BEM matrices, single precision. . . . 83
4.12 OpenMP parallelized assembly of the BEM matrices, double precision. . . . 87
4.13 OpenMP parallelized assembly of the BEM matrices, single precision. . . . 87
4.14 OpenMP vectorized assembly of the BEM matrices, double precision. . . . 88
4.15 OpenMP vectorized assembly of the BEM matrices, single precision. . . . 88
4.16 Analytic integration over triangles. . . . 92
4.17 Masked evaluationc(d>0)=a+b. . . . 97
4.18 One-level optimization loop. . . . 101
4.19 Optimization interface in BEM4I. . . . 102
4.20 Multiresolution optimization with known analytic result (ball). . . . 104
4.21 Convergence graph (ball). . . . 105
4.22 Single level optimization (ball). . . . 105
4.23 Multiresolution optimization withg= 7(duck). . . . 106
4.24 Convergence graphs forg= 2andg= 7(duck).. . . 107
4.25 Bernoulli surfaces forg∈ {2,3,4,5,6,7}(duck). . . . 108
4.26 Setting of the industrial optimization problem.. . . 110
4.27 Optimization of a high-voltage circuit disconnector. . . . 110
4.1 Element-based assembly. . . . 71
4.2 Simple quadrature. . . . 71
4.3 Vectorized quadrature. . . . 72
4.4 Naive Sauter-Schwab quadrature. . . . 73
4.5 Computation of reference quadrature nodes.. . . 74
4.6 Computation of quadrature nodes. . . . 75
4.7 Evaluation of the double-layer kernel. . . . 76
4.8 Vectorized Sauter-Schwab quadrature. . . . 77
4.9 Vectorized quadrature on MIC.. . . 81
4.10 Transfer of raw data to Intel Xeon Phi. . . . 82
4.11 Matrix assembly accelerated by the Intel Xeon Phi coprocessors. . . . 85
4.12 Scalar semi-analytic quadrature. . . . 95
4.13 Scalar evaluation of the collocation integral.. . . 95
4.14 Vectorized semi-analytic quadrature. . . . 96
4.15 Vectorized evaluation of the collocation integral. . . . 97
4.1 Convergence of the solution to the Dirichlet and Neumann problems. . . . 76
4.2 Speedup of OpenMP parallelized assembly for Vh(2×HSW, AVX2).. . . 78
4.3 Speedup of OpenMP parallelized assembly for Kh(2×HSW, AVX2).. . . 78
4.4 Speedup of OpenMP vectorized assembly for Vhand Kh. . . . 78
4.5 Speedup of CPU+MIC vs. CPU assembly for Vh. . . . 82
4.6 Speedup of CPU+MIC vs. CPU assembly for Kh. . . . 83
4.7 Speedup of OpenMP parallelized assembly for Vh(2×HSW, AVX2).. . . 86
4.8 Speedup of OpenMP parallelized assembly for Kh(2×HSW, AVX2).. . . 87
4.9 Speedup of OpenMP parallelized assembly for Vh(KNC, IMCI).. . . 87
4.10 Speedup of OpenMP parallelized assembly for Kh(KNC, IMCI). . . . 88
4.11 Speedup of OpenMP parallelized assembly for Vh(KNL, AVX-512). . . . 88
4.12 Speedup of OpenMP parallelized assembly for Kh(KNL, AVX-512). . . . 88
4.13 Speedup of OpenMP vectorized assembly for Vh. . . . 89
4.14 Speedup of OpenMP vectorized assembly for Kh. . . . 89
4.15 History of the normalized cost and gradient norm (ball). . . . 103
4.16 History of the normalized cost and gradient norm forg= 2(duck). . . . 107
4.17 History of the normalized cost and gradient norm forg= 7(duck). . . . 107
4.18 History of normalizedL2 andC cost functionals (high-voltage circuit disconnector).. . . 111
1 Introduction
Problems arising in the design and manufacture of components in many industrial areas can often be stated as minimization problems over a class of admissible shapes. Differently from optimal control problems with the state satisfying a certain partial differential equation on a fixed domain, the state variable in shape optimization problems is controlled by the shape of the computational domain. Thus, the set of admissible controls does not directly build the structure of a vector space. The pioneering attempts to mathematically describe such problems are due to Hadamard [54], where the author optimizes the shape of a clamped plate and describes the admissible domains by normal perturbations of a smooth reference boundary.
This work gave birth to shape calculus, a research field further pursued in the now classical monographs [34,114,125]. In this thesis we use the tools provided by shape calculus to solve free boundary problems of the Bernoulli type in combination with the boundary element method (BEM) for the numerical analysis of the state and adjoint problems and the subdivision surfaces representing the admissible domains.
There are several possible approaches to solving the Bernoulli type problems. First of all, the existence analysis and regularity of solutions has been studied in, e.g., [2, 5, 47, 49, 50]. The trial methods as in [18, 59, 60, 130, 131] relax the boundary value problem overdetermined by both Dirichlet and Neumann boundary condition by iteratively solving the problem with one condition only and subsequently moving the free boundary in such a way that the second condition is satisfied exactly. The approach adopted in the following text is based on the reformulation of the problem as a shape optimization problem with one of the boundary conditions enforced by minimizing a least-squares-type tracking functional. This strategy has been followed together with the fictitious domain method for solving the state and adjoint boundary value problems in [45, 62, 63, 64], wavelet-based BEM has been used in [39, 41, 42, 43, 57].
To find a minimizer of the studied cost functional one can rely on gradient-free methods.
Despite the fact that the aim of such methods is to find the global minimizer, they usually require a high number of cost evaluations which is equal to the number of possibly costly numerical analyses. The gradient-based methods, on the other hand, search for a local minimizer but provide an improved speed of convergence. The gradient information can be provided either by the first-discretize-then-optimize approach, where the state problem is first discretized and the sensitivity analysis is performed on the discrete level, or the first-optimize-then-discretize approach relying on the speed method or the perturbation of identity as described in the classical monographs [34, 65, 125] and followed for the solution of the Bernoulli problem in this thesis and in [61, 77]. For completeness, let us also mention the possibility of automatic differentiation techniques [53], which consider the whole computer program solving the state problem as a series of elementary arithmetic operations and apply the chain rule for automatized differentiation of the program with respect to the design variables.
The boundary element method employed in the thesis provides an efficient tool for solving the underlying boundary value problems if it is possible to reformulate them as boundary integral equations. Although this is especially the case for problems with (piecewise) constant material coefficients, it also allows for the natural solution of exterior problems. The boundary element approach seems natural for solving shape optimization problems since the shape of the volume is fully given by its boundary. In the context of volume finite element techniques, the boundary
perturbations have to be transferred to the volume mesh. While this can be easy for a number of academic examples, where the unknown part of the boundary can be represented, e.g., as a graph, for general surfaces this is much more difficult. One possible approach is to solve an auxiliary linear elasticity problem with boundary conditions given by the current perturbation, however, strong deformations of the boundary would still require remeshing of the volume mesh after every couple of iterations. Another option is to embed diffeomorphic shape perturbations in the variational formulation of the state problem as in [67, 110] and keep the mesh is constant during the whole optimization process.
The novelty of this thesis lies in the combination of BEM for the solution of state and adjoint boundary value problems and the subdivision surfaces used for the discretization of boundary perturbations as already presented in [10]. The idea of subdivision as a tool for constructing smooth curves as a limit of a sequence of control polygons is well-known in computer graphics and dates back to the corner cutting algorithm of Chaikin [20]. In this work we concentrate on subdivision surfaces based on triangular meshes defining quartic splines in the regular case. This subdivision technique has been originally generalized for triangular meshes of arbitrary topology by Loop [92]. In the proposed shape optimization algorithm the coarse control polygon defines the set of design parameters, while the mesh obtained afternsubdivision steps serves for the boundary element analysis. Since the subdivision basis functions define smooth perturbations, no mesh smoothing steps are necessary. The nature of the subdivision process also allows us to add design parameters when the optimum in the current design space is found and thus also serves as a globalization strategy. Further works on the topic of Loop’s subdivision include, e.g., [123, 126, 141, 142]. The combination of finite element approximation and subdivision have been discussed in [95], the isogeometric approach with the subdivision functions used for approximating both the geometry and the state variable has been suggested in [22,23]. The regularity of subdivision basis functions and the necessary approximation properties are further discussed in [8]. The isogeometric approach has also been used for shape optimization in [9, 11].
The second part of the thesis is devoted to the efficient implementation of BEM not only on modern PCs but also in the environment of High Performance Computing (HPC) centres. The drawback of the classical BEM is its quadratic complexity both in terms of the computational time and memory requirements restricting its applicability to moderate problem dimensions.
To overcome this issue, several fast BEM methods can be employed to lower the complexity to almost linear. This includes the fast multipole method [52,107,120] based on the approximation of the kernel by a truncated series, or the adaptive cross approximation (ACA) [12,119] building low-rank blocks based on an algebraic point of view. Although these sparsification methods are inevitable for large-scale engineering problems, it is still crucial to efficiently assemble the so-called non-admissible full blocks. Moreover, in the case of ACA, the low-rank approximation requires evaluation of several rows and columns of every admissible block, which relies on the standard full assembly. In the thesis we thus concentrate on the acceleration of the standard BEM assembly in shared memory.
The distribution of the workload among available CPU cores by, e.g., OpenMP pragmas has become more or less standard in similar scientific codes. Although such parallelization techniques lead to a significant speedup, overlooking further acceleration by utilizing the Single Instruction Multiple Data (SIMD) concept can hardly exploit the full potential of modern CPUs. The first
Intel’s effort to employ vector instructions on multiple operands dates back to 1997 with the introduction of the MMX instructions set reusing the existing scalar registers for integer opera- tions only. The SSE-SSE4 instruction sets introduced in 1999-2006 added eight 128-bit registers designed to work in parallel with four single-precision or two double-precision floating-point operands. The AVX-AVX2 sets supported by Intel processors since 2011-2013 further increase the size of the vector register to 256 bit and incrementally add more instructions including the three-operand Fused Multiply Add (FMA) combining addition and multiplication in one step.
The need for code vectorization is even more apparent in connection with the Many Integrated Core (MIC) architecture represented by the Intel Xeon Phi (co)processors. Currently, the Xeon Phi coprocessors support the Initial Many Core Instructions (IMCI) set doubling the size of the AVX2 registers. Although programming for an accelerator may seem a daunting task, the MIC architecture allows to reuse the majority of the already existing CPU code, which is in contrast with the GP-GPU acceleration on graphics processing units. Moreover, the current AVX-512 instruction set providing 512-bit registers is supported (at least partially) both by the Skylake CPU architecture and the Knights Landing (KNL) MIC architecture. Contrary to the Knights Corner (KNC) MIC copocessors, the KNL version is also available as a standalone host processor, possibly showing the future of high performance computing.
The ideas of vectorization in connection with BEM have already been discussed in several publications. The automatic vectorization based on the capabilities of modern compilers is discussed in [29], the direct use of vector intrinsic functions for the BEM matrix assembly, however, without a detailed look on the integration of singular kernels, is presented in [69], and the vectorization of the evaluation of the representation formula is shown in [93]. In the thesis we present two different strategies for code vectorization, namely by using OpenMP 4.0 pragmas [101, 109] and the Vc library wrapping the vector intrinsic functions [84, 98]. Moreover, the computationally most intensive parts allow the offload to MIC coprocessors further speeding up the assembly [101]. In contrast to the previously mentioned references we also discuss the treatment of singularities in the underlying surface integrals both by the semi-analytic [119] and fully numeric quadrature schemes [121]. The presented approaches build the core of the BEM4I library [97] developed at the IT4Innovations Czech National Supercomputing Center.
The structure of the thesis follows. Section 2 first reviews basic properties of function spaces appearing throughout the thesis and presents the shape optimization problem under consideration together with its solvability result under reasonable regularity assumptions and the derivation of the shape derivative. The final part of the section summarizes the boundary integral approach to solving boundary value problems. Section 3 is devoted to the representation of smooth surfaces generated by Loop’s subdivision of quartic box splines. These surfaces are used to represent the unknown shape and to discretize the admissible boundary perturbations.
In Section 4 we first review the discretization of boundary integral equations, i.e., the boundary element method. We present two approaches to the assembly of system matrices, namely the semi-analytic and fully numerical quadrature, and describe its efficient realization on modern hardware including the acceleration by the Intel Xeon Phi coprocessors. With all the necessary tools available we present the discrete version of the optimization problem validated by solving two academic examples and a real-life experiment originating in the ABB Corporate Research Center Switzerland.
2 The Bernoulli free boundary problem
Before we concentrate on the Bernoulli problem itself, we introduce the necessary function spaces and review the theory of boundary integral equations.
2.1 Function spaces
In this section we introduce function spaces necessary for the weak formulation of elliptic boundary value problems, the associated boundary integral equations, and for the boundary transformations used in the shape optimization procedure. In the following we denote by Ω ⊂ Rd a bounded domain. For a more detailed treatment of the topic the reader is referred to [3, 85, 86, 105].
2.1.1 Continuous functions
We denote by C0(Ω) or C(Ω) the set of functions u:Ω → R continuous in Ω. By Ck(Ω) for k∈N0 we denote the space of k-times continuously differentiable functions inΩ, i.e.,
Ck(Ω) :={u∈C(Ω) :Dαu∈C(Ω) for all α,|α| ≤k} with the multiindex
α= [α1, . . . , αd]∈Nd0 of the length |α|:=α1+· · ·+αd and the partial derivative operator
Dαu:= ∂|α|
∂xα11· · ·∂xαddu, ∂0
∂x0iu:=u.
For k∈N0, the set of k-times continuously differentiable functions with compact support in Ω is denoted as
C0k(Ω) :={u∈Ck(Ω) : suppu:={x∈Ω:u(x)̸= 0} ⊂Ω is compact}. We define
C∞(Ω) := ⋂
k∈N
Ck(Ω), C0∞(Ω) := ⋂
k∈N
C0k(Ω).
Definition 2.1. We say that a function u is called uniformly continuous on a set M ⊂ Rd (and denote u∈C(M)) if for every ε > 0 there existsδ >0 such that for all points x,y ∈M satisfying ∥x−y∥< δ it holds that|u(x)−u(y)|< ε.
By C0(Ω) or C(Ω) we understand functions bounded and uniformly continuous in Ω, and thus continuously extensible to∂Ω. Fork∈N0 we introduce the spaces
Ck(Ω) :={u∈C(Ω) :Dαu∈C(Ω) for allα,|α| ≤k}.
Equipped with the norm
∥u∥Ck(Ω):= ∑
α,|α|≤k
max
x∈Ω
|Dαu(x)|,
where we identifyDαu with its extension to Ω, these spaces are complete, i.e., Banach spaces.
Again,
C∞(Ω) := ⋂
k∈N
Ck(Ω) C0∞(Ω) := ⋂
k∈N
C0k(Ω).
Definition 2.2. A set K ⊂ C(Ω) is called equicontinuous if the parameter δ = δ(ε) from Definition 2.1 can be chosen independently ofu∈K.
Definition 2.3. Let X denote a metric space. The set Y ⊂ X is relatively compact in X if for every sequence (yn)⊂Y there exists a subsequence (ynk) and an element y ∈X such that ynk →y in X.
Theorem 2.4(Arzelà-Ascoli). A setK ⊂C(Ω)is relatively compact if and only if it is bounded and equicontinuous.
Proof. The proof of this standard result can be found, e.g., in [86].
Corollary 2.5. A set K ⊂ Ck(Ω) is relatively compact if and only if it is bounded in Ck(Ω) and the sets
Ks:={Dαu:u∈K,|α|=s} are equicontinuous for all s≤k.
Proof. Let (un) be a sequence in K. From boundedness in Ck(Ω) it follows, that it is also bounded in C(Ω). Since the set K0 is equicontinuous, it follows from Theorem 2.4 that there exists a subsequence (unk) convergent inC(Ω). Using the same arguments, i.e., boundedness of K1 inC(Ω) and equicontinuity of
{Dαunk:k∈N,|α|= 1}⊂K1,
leads to a subsequence (unkℓ)⊂(unk) for which it holds that (Dαunkℓ) is convergent inC(Ω) for any α, |α|= 1. In other words, (unkℓ) is convergent in C1(Ω). Repeating this process finitely many times leads to the assertion.
Finally, we define the Banach space of Hölder continuous functions for k∈N0,λ∈(0,1] as Ck,λ(Ω) :={u∈Ck(Ω) :Hα,λ(u)<∞for all α,|α|=k}
with
Hα,λ(u) := sup
x,y∈Ω x̸=y
|Dαu(x)−Dαu(y)|
∥x−y∥λ and the norm
∥u∥Ck,λ(Ω):=∥u∥Ck(Ω)+ ∑
α,|α|=k
Hα,λ(u).
The setC0,1(Ω) is the space of Lipschitz continuous functions.
2.1.2 Lebesgue and Sobolev spaces
Forp∈[1,∞)∪ {∞} we define the Lebesgue spaces
Lp(Ω) :={u measurable on Ω:∥u∥Lp(Ω)<∞} with
∥u∥Lp(Ω):=
(∫
Ω
|u(x)|pdx )1/p
, ∥u∥L∞ := ess sup
Ω
|u|. (2.1)
Equipped with the respective norms (2.1) these spaces are Banach spaces. With the inner product
⟨u, v⟩L2(Ω) :=
∫
Ω
u(x)v(x) dx
inducing the ∥ · ∥L2(Ω) norm, the spaceL2(Ω) is also a Hilbert space. Additionally, we define the space L1loc(Ω) of functionsf ∈L1(O) for any compact subsetO ⊂Ω.
To define weak solutions to boundary value problems it is necessary to introduce appropriate Sobolev spaces. On C∞(Ω),C0∞(Ω) we define for p∈[1,∞) and s∈Nthe Sobolev norms
∥u∥Ws,p(Ω):=
⎛
⎝
∑
α,|α|≤s
Dαu
p Lp(Ω)
⎞
⎠
1/p
, ∥u∥Ws,∞(Ω) := max
α,|α|≤s
Dαu
L∞(Ω) (2.2) and the spacesWs,p(Ω),W0s,p(Ω) withp∈[1,∞)∪ {∞} as the completion ofC∞(Ω),C0∞(Ω), respectively, with respect to the norms (2.2).
Remark 2.6. For Lipschitz domains (see Definition 2.7), the previous definitions can be shown to be equivalent to
Ws,p(Ω) :={u∈Lp(Ω) : Dαu∈Lp(Ω) for allα,|α| ≤s}. (2.3) The partial derivatives from the definition (2.3) are understood in the distributional sense, for more details see, e.g., [3, 96].
For p ∈ [1,∞) and s = k+κ with k ∈ N0, κ ∈ (0,1), the definition of Ws,p(Ω) can be extended by considering the Sobolev–Slobodeckij norms
∥u∥Ws,p(Ω):=(∥u∥pWk,p(Ω)+|u|pWs,p(Ω)
)1/p
(2.4) with the semi-norm
|u|Ws,p(Ω):=
⎛
⎝
∑
α,|α|=k
∫
Ω
∫
Ω
|Dαu(x)−Dαu(y)|p
∥x−y∥d+pκ dydx
⎞
⎠
1/p
. For convenience, we denote the L2(Ω) based Hilbert spaces
Hs(Ω) :=Ws,2(Ω), H0s(Ω) :=W0s,2(Ω)
Figure 2.1: Local parametrization of∂Ω.
with the inner product
⟨u, v⟩Hs(Ω):= ∑
α,|α|≤k
⟨Dαu, Dαv⟩L2(Ω)
fors∈N and
⟨u, v⟩Hs(Ω):= ∑
α,|α|≤k
⟨Dαu, Dαv⟩L2(Ω)
+ ∑
α,|α|=k
∫
Ω
∫
Ω
(Dαu(x)−Dαu(y))(Dαv(x)−Dαv(y))
∥x−y∥d+2κ dydx
fors=k+κ with k∈ N0, κ ∈(0,1). The above inner products induce the norms (2.2), (2.4) forp= 2.
2.1.3 Continuous functions on manifolds
To introduce relevant function spaces on the boundary of a domain we first need to address the topic of its regularity. Since we deal with bounded domainsΩ ⊂Rd, we can assume that the boundary∂Ω is a compact set.
Definition 2.7. A domainΩ⊂Rd is a Ck,λ domain, orΩ∈Ck,λ, if there exists a finite open cover of∂Ωdenoted by{Oi}ni=1 such that for everyi∈ {1, . . . , n}there exist a Cartesian system of coordinates
[yi1, . . . , yd−1i , ydi] = [yi, ydi], where yi := [y1i, . . . , yid−1], numbersεi, δi∈R+, and a function
ai:Rd−1 →R, ai∈Ck,λ({yi:∥yi∥ ≤δi})
satisfying
Oi∩∂Ω={(yi, ydi) :∥yi∥< δi, ydi =ai(yi)},
Ω⊃ {(yi, ydi) :∥yi∥< δi, ai(yi)< yid< ai(yi) +εi}, (Rd\Ω)⊃ {(yi, ydi) :∥yi∥< δi, ai(yi)−εi< yid< ai(yi)}.
A C0,1 domain is called a Lipschitz domain.
Remark 2.8. In the following we shall assume that all domains are at least Lipschitz if not stated otherwise.
Definition 2.9. We say that a domain Ω ⊂ Rn is L-Lipschitz with L ∈ R+ if it is Lipschitz and there exists a finite open cover{Oi}ni=1 such that for the functionsai it holds
|ai(x)−ai(y)| ≤L∥x−y∥ for all x,y∈ {yi:∥yi∥ ≤δi} with the notation as in Definition 2.7.
Remark 2.10. It can be shown that a domain Ω ⊂ Rn is L-Lipschitz if and only if it satisfies the so-called ε-cone property, see [21, 114].
With each cover {Oi}ni=1 of∂Ω we associate the partition of unity, i.e., the functions
{λi:Rd→R}ni=1 (2.5)
such that
(∀i∈ {1, . . . , n})(∀x∈Rd) : 0≤λi(x)≤1,
∀i∈ {1, . . . , n}:λi ∈C0∞(Rn),suppλi⊂Oi,
∀x∈∂Ω:
n
∑
i=1
λi(x) = 1.
Remark 2.11. In the following we also deal with vector fields defining function spaces [X]d. These can be understood componentwise as functions in the space X. IfX is a normed space, the norm in [X]d can be defined, e.g., as
∥V∥[X]d =∥(V1, . . . , Vd)∥[X]d:=
( d
∑
i=1
∥Vi∥2X )1/2
.
Where these spaces can be distinguished from context, we use the standard scalar notation.
Remark 2.12. It can be shown that for every Ck,λ domainΩ there exist a finite family of open sets{Oi}ni=1 covering∂Ω and mappings
ψi∈Ck,λ(Bd(0,1), Oi)
, such that [ψi]−1∈Ck,λ(Oi, Bd(0,1))
and
[ψi]−1(∂Ω∩Oi) =Bd(0,1)∩[Rd−1× {0}] =:B0d(0,1)≡Bd−1(0,1), [ψi]−1(Ω∩Oi) =Bd(0,1)∩[Rd−1×(−∞,0)],
[ψi]−1(Rd\Ω∩Oi) =Bd(0,1)∩[Rd−1×(0,∞)].
The pairs (Oi,[ψi]−1) are called charts, the collection of all charts defines an atlas. For an illustration see Figure 2.1, for more details and examples we also refer to [34].
Taking into account the partition of unity, a function u:∂Ω →Rcan be decomposed as u(x) =
n
∑
i=1
λi(x)u(x) =
n
∑
i=1
ui(x) with ui:=λiu.
Transferringui:∂Ω→R to the parameter domain byψi leads to
uˆi:Bd−1(0,1)→R, uˆi(y) :=ui(ψi(y)) =λi(ψi(y))u(ψi(y)). (2.6) We say thatu is continuous on∂Ω, oru∈C(∂Ω) if and only if
uˆi∈C(Bd−1(0,1)) for all i∈ {1, . . . , n}.
Since the chain rule for partial derivatives ofuˆi yields
∂uˆi
∂yj
(y) =
d
∑
k=1
∂ui
∂xk
(ψi(y))∂ψki
∂yj
(y),
in general it makes sense to define forΩ∈Ck,λ the spaces of maximal smoothness
Ck,λ(∂Ω) :={u∈C(∂Ω) :uˆi ∈Ck,λ(Bd−1(0,1)) for alli∈ {1, . . . , n}} (2.7) with the norm
∥u∥Ck,λ(∂Ω) :=
n
∑
i=1
∥uˆi∥
Ck,λ(Bd−1(0,1)). (2.8)
2.1.4 Lebesgue and Sobolev spaces on manifolds
To define the Lebesgue spaces on boundaries of domains we follow the definition of continuous functions. We say thatu∈Lp(∂Ω),p∈[1,∞)∪ {∞} if for every i∈ {1, . . . , n}the compound functionuˆi from (2.6) belongs toLp(Bd−1(0,1)). Forp∈[1,∞) the norm is defined as
|||u|||Lp(∂Ω):=
( n
∑
i=1
∫
Bd−1(0,1)
|u(ψi(y))|pdy )1/p
. (2.9)
Definition 2.13. Let Ω ∈ C0,1 and {λi}ni=1 denote a partition of unity corresponding to the open cover{Oi}ni=1and Rd∋ed:= [0, . . . ,1]T. Foru∈L1(∂Ω) we define the surface integral as
∫
∂Ω
u(x) dsx:=
n
∑
i=1
∫
Bd−1(0,1)
λi(ψi(y))u(ψi(y))|detDψi(y)|[Dψi(y)]−Teddy.
Remark 2.14. For Ω⊂R3 the previous definition is equivalent to
∫
∂Ω
u(x) dsx=
n
∑
i=1
∫
B2(0,1)
λi(ψi(y))u(ψi(y))
∂ψi
∂y1(y)×∂ψi
∂y2(y)
dy (2.10)
with the component-wise derivatives
∂ψi
∂yj :=
[∂ψi1
∂yj, ∂ψ2i
∂yj, ∂ψi3
∂yj ]T
.
Remark 2.15. The norms (2.9) dependent on the parametrization of the boundary are equivalent to the norm
∥u∥Lp(∂Ω) :=
(∫
∂Ω
|u(x)|pdsx
)1/p
. (2.11)
The space L2(∂Ω) is a Hilbert space with the inner product
⟨u, v⟩L2(∂Ω):=
∫
∂Ω
u(x)v(x) dsx. (2.12)
Similarly as in the case of continuous functions, for k∈Nand Ω∈Ck−1,1 we only consider the Sobolev spacesHs(∂Ω) of the order s∈R,|s| ≤k. For 0≤s≤kwe set
Hs(∂Ω) :={u:uˆi ∈Hs(Bd−1(0,1)) for alli∈ {1, . . . , n}}
with the norm
|||u|||Hs(∂Ω):=
( n
∑
i=1
∥ui◦ψi∥2Hs(Bd−1(0,1))
)1/2
. (2.13)
Remark 2.16. Fors= 0 the norm (2.13) is equivalent to theL2(∂Ω) norm (2.11). Fors∈(0,1) we have an equivalent Sobolev–Slobodeckij norm
∥u∥Hs(∂Ω):=(∥u∥2L2(∂Ω)+|u|2Hs(∂Ω)
)1/2
with the semi-norm
|u|Hs(∂Ω):=
( ∫
∂Ω
∫
∂Ω
|u(x)−u(y)|2
∥x−y∥d−1+2sdsydsx
)1/2
. For s <0 we introduce the dual spaces
Hs(∂Ω) := [H−s(∂Ω)]∗. (2.14)
Theorem 2.17. For k∈N,λ∈(0,1],0≤s≤k, and a Ck−1,λ domain Ω it holds Hs(∂Ω)↪→L2(∂Ω) = [L2(∂Ω)]∗ ↪→H−s(∂Ω),
where both embeddings are continuous and dense, and the equality of L2(∂Ω) and its dual is understood in the sense of Riesz.
Proof. For details on Gelfand triples see, e.g., [121, Section 2.1.2.4].
The preceding theorem ensures that the inner product (2.12) can be continuously extended to the duality pairing onH−s(∂Ω)×Hs(∂Ω), i.e.,
⟨u, v⟩∂Ω :=⟨u, v⟩H−s(∂Ω)×Hs(∂Ω)= lim
∫
∂Ω
un(x)v(x) dsx
for a sequence (un)⊂L2(∂Ω),un→u in the standard dual norm ofH−s,
∥u∥H−s(∂Ω) := sup
v∈Hs(∂Ω)\{0}
⟨u, v⟩∂Ω
∥v∥Hs(∂Ω)
. (2.15)
LetΓ0 ⊂∂Ω denote a patch of the boundary. Fors >0 we define the space Hs(Γ0) :={u= ˜u|Γ0: ˜u∈Hs(∂Ω)}
with the norm
∥u∥Hs(Γ0):= inf
u∈H˜ s(∂Ω) : ˜u|Γ
0=u
∥˜u∥Hs(∂Ω)
and the space
H˜s(Γ0) :={u= ˜u|Γ0: ˜u∈Hs(∂Ω),supp ˜u⊂Γ0}.
Fors <0 we define
Hs(Γ0) := [H˜−s(Γ0)]∗, H˜s(Γ0) := [H−s(Γ0)]∗. Remark 2.18. In case thatΓ0 is a closed surface we trivially have
Hs(Γ0) =H˜s(Γ0).
For a partitioning of the boundary
∂Ω =
n
⋃
i=1
Γi
with disjoint and smooth open panelsΓi we define fors >0 the spaces Hpws (∂Ω) :={u∈L2(∂Ω) :u|Γi ∈Hs(Γi) for i= 1, . . . , n}
with the norm
∥u∥Hs
pw(∂Ω):=
( n
∑
i=1
∥u|Γi∥2Hs(Γi)
)1/2
. Fors <0 we introduce
Hpws (∂Ω) :=
n
∏
i=1
H˜s(Γi) and the norm
∥u∥Hs
pw(∂Ω):=
( n
∑
i=1
∥u∥2
H˜s(Γi)
)1/2
.
Proposition 2.19. Let Ω ⊂Rd denote a Lipschitz domain. Then there exists a unique linear continuous mapping
γ0:H1(Ω)→H1/2(∂Ω) satisfying
γ0u=u|∂Ω for all u∈C(Ω).
Remark 2.20. The function γ0uis called the (Dirichlet) trace of the function u∈H1(Ω).
Proof. See, e.g., [24, 96]
Due to the result of Rademacher, for Lipschitz domain the unit outward normal vector n can be defined almost everywhere on ∂Ω and n∈L∞(∂Ω). For u ∈H2(Ω) we can define the normal derivative as
γ1u:=⟨γ0∇u,n⟩ ∈L2(∂Ω) = [L2(∂Ω)]∗ ↪→H−1/2(∂Ω),
with the trace operator γ0 understood as a component-wise application of γ0. To formulate boundary integral equations for the Laplace equation it is necessary to generalize the concept, since in general the solution is sought in
H∆1(Ω) :={u∈H1(Ω) : ∆u∈L2(Ω)}
with the Laplace operator understood in the distributional sense. It can be shown that the following generalized Green’s first identity holds, see, e.g., [121, 137].
Theorem 2.21 (Green’s first formula). Let Ω denote a Lipschitz domain and u ∈ H∆1(Ω).
Then there exists a unique linear continuous mapping γ1:H∆1(Ω)→H−1/2(∂Ω) satisfying
γ1u= ∂u
∂n ∈L2(∂Ω) = [L2(∂Ω)]∗ ↪→H−1/2(∂Ω) for all u∈C1(Ω) and it holds
∫
Ω
∆u(x)v(x) dx=⟨γ1u, γ0v⟩∂Ω−
∫
Ω
⟨∇u(x),∇v(x)⟩dx for all v∈H1(Ω). (2.16) Corollary 2.22(Green’s second formula). Let Ωdenote a Lipschitz domain andu, v∈H∆1(Ω).
Then it holds
∫
Ω
∆u(x)v(x)−u(x)∆v(x) dx=⟨γ1u, γ0v⟩∂Ω− ⟨γ1v, γ0u⟩∂Ω. (2.17) Remark 2.23. The functionγ1uis called the generalized normal derivative or the Neumann trace of the function u∈H∆1(Ω).
2.2 Continuous optimization problem
We are now in position to define the optimization problem treated in the following sections.
Firstly, we define the corresponding forward problem and later concentrate on the solvability of the inverse problem reformulated as a shape optimization problem.
The classical formulation of the direct problem to determine the electrostatic potential uΩ
in absence of volume charges reads
find uΩ ∈C2(Ω)∩C(Ω), such that
⎧
⎪⎪
⎨
⎪⎪
⎩
−∆uΩ = 0 inΩ, uΩ =h on Γ0, uΩ = 0 on Γf,
(2.18)
with an annular domainΩ ⊂R3 with its boundary ∂Ω decomposed into separated subsets Γ0 andΓf. We denote byΩ0 the domain enclosed by the componentΓ0 and byD⊂R3 a bounded hold-all domain satisfyingΩ⊂Dand dist(Γf, ∂D∪Γ0)≥δ >0. See Figure 2.2 for the described configuration. Additionally, we assumeh∈H1/2(Γ0).
The weak formulation of (2.18) reads
finduΩ∈ {u∈H1(Ω) :γ0u|Γ0 =h, γ0u|Γf = 0}, such that
∫
Ω
⟨∇uΩ(x),∇v(x)⟩dx= 0 for all v∈H01(Ω). (P(Ω)) Due to standard arguments the problem (P(Ω)) is uniquely solvable and the solution depends continuously on the Dirichlet data.
The classical formulation of the inverse problem, more specifically the exterior Bernoulli free boundary problem, reads
find Ω∈ O, uΩ ∈C2(Ω)∩C1(Ω), such that
⎧
⎪⎪
⎪⎪
⎪⎪
⎨
⎪⎪
⎪⎪
⎪⎪
⎩
−∆uΩ = 0 inΩ, uΩ =h on Γ0, uΩ = 0 on Γf,
∂uΩ
∂n =−g on Γf,
(2.19)
Figure 2.2: Exterior Bernoulli problem configuration.
