Efficient PDE constrained shape optimization based on Steklov-Poincaré type metrics

# Efficient PDE constrained shape optimization based on Steklov-Poincaré type metrics

## Abstract

Recent progress in PDE constrained optimization on shape manifolds is based on the Hadamard form of shape derivatives, i.e., in the form of integrals at the boundary of the shape under investigation, as well as on intrinsic shape metrics. From a numerical point of view, domain integral forms of shape derivatives seem promising, which rather require an outer metric on the domain surrounding the shape boundary. This paper tries to harmonize both points of view by employing a Steklov-Poincaré type intrinsic metric, which is derived from an outer metric. Based on this metric, efficient shape optimization algorithms are proposed, which also reduce the analytical labor, so far involved in the derivation of shape derivatives.

P

DE constrained shape optimization, optimization on shape manifolds.

## 1 Introduction

Shape optimization is of interest in many fields of application – in particular in the context of partial differential equations (PDE). As examples, we mention aerodynamic shape optimization [22], acoustic shape optimization [30] or optimization of interfaces in transmission problems [10, 18, 20] and in electrostatics [4]. In industry, shapes are often represented within a finite dimensional design space. However, often this reduction is felt as being too restrictive [27], which motivates shape optimization based on shape calculus. Major effort in shape calculus [7, 26] has been devoted towards expressions for shape derivatives in so-called Hadamard-form, i.e., in boundary integral form. It is known that the second order shape derivative, formerly coined as shape Hessian, is nonsymmetric in general, which for a long time has been an obstacle for algorithmic developments in shape optimization in the fashion of nonlinear programming. Recently [23, 24, 25], shape optimization has been considered as optimization on Riemannian shape manifolds, which enables design and analysis of NLP-like algorithms including one-shot sequential quadratic programming and theoretical insights into the structure of the second order shape derivative in comparison to the Riemannian shape Hessian. Coercivity results for shape Hessians for elliptic problems can be found in [8]. The scalar product used in this work is in line with these results.

On the other hand, it is often a very tedious, not to say painful, process to derive the boundary formulation of the shape derivative. Along the way, there frequently appears a domain formulation in the form of an integral over the whole domain as an intermediate step. Recently, it has been shown that this intermediate formulation has numerical advantages [5, 10, 12, 20]. In [14], also practical advantages of the domain shape formulation have been demonstrated, since it requires less smoothness assumptions. Furthermore, the derivation as well as the implementation of the domain integral formulation requires less manual and programming work. Thus, there arises the natural goal of combining the favorable domain integral formulation of shape derivatives with the favorable NLP-type optimization strategies on shape manifolds, which seem so far tightly coupled with boundary integral formulations of shape derivatives. This publication aims at demonstrating that this coupling is indeed possible and that it naturally leads to a novel family of Poincaré-Steklov type metrics on shape manifolds. In contrast to [24] this work consciously avoids surface formulations of shape derivatives in order to provide more handy optimization algorithms.

The paper is organized in the following way. First, in section 2, we set up notation and terminology and formulate the model problem. In section 3, we discuss generalized Poincaré-Steklov operators as the basis for Riemannian metrics for shape manifolds. Section 6 is devoted to the set of all shapes in the context of the novel metric introduced in section 3. Section 4 rephrases NLP-like optimization algorithms on shape manifolds within the framework of domain integral formulations of shape derivatives. Finally, section 5 discusses algorithmic and implementation details, as well as, numerical results for a parabolic transmission shape optimization problems.

## 2 Problem Formulation

We first set up notation and terminology in shape calculus. Then we recall the model problem in [24], which is motivated by electrical impedance tomography and given by a parabolic interface shape optimization problem.

### 2.1 Notations in shape calculus

Let and . We denote by a bounded domain with Lipschitz boundary and by a real-valued functional depending on it. Moreover, let be a family of bijective mappings such that . This family transforms the domain into new perturbed domains with and the boundary into new perturbed boundaries with . If you consider the domain as a collection of material particles, which are changing their position in the time-interval , then the family describes the motion of each particle, i.e., at the time a material particle has the new position with . The motion of each such particle could be described by the velocity method, i.e., as the flow determined by the initial value problem

 dξ(t,x)dt=V(ξ(t,x))ξ(0,x)=x (1)

or by the perturbation of identity, which is defined by where denotes a sufficiently smooth vector field. We will use the perturbation of identity throughout the paper. The Eulerian derivative of at in direction is defined by

 DJ(Ω)[V]:=limt→0+J(Ωt)−J(Ω)t. (2)

The expression is called the shape derivative of at in direction and shape differentiable at if for all directions the Eulerian derivative (2) exists and the mapping is linear and continuous. For a thorough introduction into shape calculus, we refer to the monographs [7, 26]. In particular, [31] states that shape derivatives can always be expressed as boundary integrals due to the Hadamard structure theorem. The shape derivative arises in two equivalent notational forms:

 DJΩ[V] :=∫ΩF(x)V(x)dx (domain formulation) (3) DJΓ[V] :=∫Γf(s)V(s)⊤n(s)ds (boundary formulation) (4)

where is a (differential) operator acting linearly on the perturbation vector field and with

 DJΩ[V]=DJ(Ω)[V]=DJΓ[V]. (5)

The boundary formulation (4), , acting on the normal component of has led to the interpretation as tangential vector of a corresponding shape manifold in [23].

### 2.2 PDE model definition

We use the same model problem as in [24], which is briefly recalled. Let this domain be an open subset of and split into the two disjoint subdomains such that , and where the interior boundary is assumed to be smooth and variable and the outer boundary Lipschitz and fixed. An example of such a domain is illustrated in figure 1.

###### Remark 1

In the shape optimization method proposed in this work the topology of the domain is fixed. This means we do not consider topology optimization.

The parabolic PDE constrained shape optimization problem is given in strong form by

 minJ(Ω)=j(Ω)+ jreg(Ω):=∫T0∫Ω(y−¯y)2dxdt+μ∫Γint1ds (6) s.t. ∂y∂t−div(k∇y) =fmodelin Ω×(0,T] (7) y =1on Γtop×(0,T] (8) ∂y∂n =0on (Γbottom∪Γleft∪Γright)×(0,T] (9) y =y0in Ω×{0} (10)

where

 k≡{k1=const. in Ω1×(0,T]k2=const. in Ω2×(0,T]

and denotes the unit outer normal to at . Of course, the formulation (7) of the differential equation is to be understood only formally because of the jumping coefficient . We observe that the unit outer normal to is equal to , which enables us to use only one normal for the subsequent discussions. Furthermore, we have interface conditions at the interface . We formulate explicitly the continuity of the state and of the flux at the boundary as

 Missing or unrecognized delimiter for \left (11)

where the jump symbol denotes the discontinuity across the interface and is defined by where and . The perimeter regularization, , with in the objective (6) is frequently used in this kind of problems. In [28] a weaker but more complicated regularization is instrumental in order to show existence of solutions.

We assume and . In our setting, the boundary value problem (7-11) is written in weak form

 amodel(y,p)=bmodel(p,p1p2), ∀p∈W(0,T;H1(Ω)) (12)

and for all , as in [24]. For properties of the function spaces, we refer the reader to the literature, e.g. [11, 29]. The bilinear form in (12) is achieved by applying integration by parts on and on . Thus, we get

 Missing or unrecognized delimiter for \left (13)

The linear form in (12) is given by

 bmodel(p,p1,p2):=bmodel1(p)+bmodel2(p1,p2) (14)

where

 bmodel1(p) :=∫T0∫Ωfmodelpdxdt, (15) bmodel2(p1,p2) :=∫T0∫Γtopp1(y−1)dsdt+∫T0∫Γout∖Γtopp2∂y∂ndsdt. (16)

In the following, we assume for the observation . The Lagrangian of (6-11) is defined as

 L(Ω,y,p):=J(Ω)+amodel(y,p)−bmodel(p,p1,p2) (17)

where is defined in (6), in (13) and in (14-16).

The adjoint problem, which we obtain from differentiating the Lagrangian with respect to , is given in strong form by

 −∂p∂t−div(k∇p) =−(y−¯¯¯y)in Ω×[0,T) (18) p =0in Ω×{T} (19) Missing or unrecognized delimiter for \left =0on Γint×[0,T) (20) Missing or unrecognized delimiter for \left =0on Γint×[0,T) (21) ∂p∂n =0on (Γbottom∪Γleft∪Γright)×[0,T) (22) p =0on Γtop×[0,T) (23) p1 =−k1pon (Γbottom∪Γleft∪Γright)×[0,T) (24) p2 =k1∂p∂non Γtop×[0,T) (25)

and the state equation, which we get by differentiating the Lagrangian with respect to , is given in strong form by

 ∂y∂t−div(k∇y)=fmodelin Ω×(0,T]. (26)

As mentioned earlier, in many cases, the shape derivative arises in two equivalent forms. If we consider the objective (6) without the perimeter regularization , the shape derivative can be expressed as an integral over the domain as well as an integral over the interface . Assume that a solution of the parabolic PDE problem (7-11) exists and is at least in . Moreover, assume that the adjoint equation (18-23) admits a solution . Then the shape derivative of the objective without perimeter regularization, i.e., the shape derivative of at in the direction expressed as an integral over the domain is given by

 DjΩ[V]=∫T0∫Ω−k∇y⊤(∇V+∇V⊤)∇p−p(∇fmodel)⊤V+div(V)(12(y−¯¯¯y)2+∂y∂tp+k∇y⊤∇p−fmodelp)dxdt. (27)

This domain integral allows us to calculate the boundary expression of the shape derivative, which is given by

 Missing or unrecognized delimiter for \left (28)

The derivations are very technical. Note that we need a higher regularity of and to provide the boundary shape derivative expression (28). More precisely, has to be an -function having weak first derivatives in and has to be an element of . Here denotes the dual space of . We achieve (27) by an application of the theorem of Correa and Seeger [6, theorem 2.1] and (28) by an application of integration by parts. We refer the reader for its derivations to [24]. By combining theorem 2.1 and 2.2 in [24] with proposition 5.1 in [19] we get the following two expressions for the shape derivative of the objective (with perimeter regularization) at in the direction :

 DjΩ[V]+Djreg(Ω)[V]=DJ(Ω)[V]=DjΓint[V]+Djreg(Ω)[V] (29)

with

 Djreg(Ω)[V]=∫Γint⟨V,n⟩μκds (30)

where denotes the curvature corresponding to the normal .

## 3 Steklov-Poincaré type metrics on shape manifolds

We first discuss the definition of shape manifolds and metrics. Then, we introduce novel metrics dovetailed to shape optimization based on domain formulations of shape derivatives.

### 3.1 Shape manifolds

As pointed out in [23], shape optimization can be viewed as optimization on Riemannian shape manifolds and resulting optimization methods can be constructed and analyzed within this framework, which combines algorithmic ideas from [1] with the differential geometric point of view established in [16]. As in [23], we first study connected and compact subsets with and boundary where denotes a bounded domain with Lipschitz-boundary (cf. figure 1). We now identify the variable boundary with a simple closed curve . Additionally, we need to describe a space including all feasible shapes and the corresponding tangent spaces. In [16], this set of smooth boundary curves is characterized by

 Be(S1,R2):=Emb(S1,R2)/Diff(S1), (31)

i.e., as the set of all equivalence classes of embeddings of into the plane (), where the equivalence relation is defined by the set of all re-parameterizations, i.e., diffeomorphisms of into itself (). A particular point on the manifold is represented by a curve . Because of the equivalence relation (), the tangent space is isomorphic to the set of all normal vector fields along , i.e.,

 TcBe≅{h:h=αn,α∈C∞(S1,R)} (32)

where is the unit exterior normal field of the shape defined by the boundary such that for all and denotes the circumferential derivative as in [16]. Several intrinsic metrics are discussed in [16], among which the following Sobolev metric seems the most natural intrinsic one from a numerical point of view. For , the Sobolev metric is induced by the scalar product

 g1:TcBe×TcBe→R,(h,k)↦((id−A△c)α,β)L2(c) (33)

where and denote two elements from the tangent space at and denotes the Laplace-Beltrami operator on the surface . In [16] it is shown that the condition guarantees that the scalar product defines a Riemannian metric on and thus, geodesics can be used to measure distances.

With the shape space and its tangent space in hand we can now form the Riemannian shape gradient corresponding to a shape derivative given in the form

 dJ(Ω)[V]=∫cψ⟨V,n⟩ds. (34)

In our model setting the objective function is given in (6) and its shape derivative in (29). Finally, the Riemannian shape gradient with respect to the Riemannian metric is obtained by

The metric , which is also used in [24], necessitates a shape derivative in Hadamard form as an efficient means to solve linear systems involving the Laplace Beltrami operator in surfaces. All of that is certainly not impossible but requires computational overhead which we can get rid of by usage of the metric discussed below. We compare the algorithmic aspects of both approaches below in section 5.

### 3.2 Steklov-Poincaré type Riemannian metrics

The ideal Riemannian metric for shape manifolds in the context of PDE constrained shape optimization problems is to be derived from a symmetric representation of the second shape derivative in the solution of the optimization problems. Often, this operator can be related to the Dirichlet to Neumann map, aka Steklov-Poincaré operator, or the Laplace-Beltrami operator [21]. If one aims at mesh independent convergence properties, one of these two will be appropriate in most cases. Since it can be observed that the Laplace-Beltrami operator is spectrally equivalent to the square of the Steklov-Poincaré operator, the latter operator seems to be more fundamental and we will focus on it as a basis for the scalar product on . Another advantage of this operator is that is blends well in with a corresponding mesh deformation strategy.

Most often, the Dirichlet to Neumann map is associated with the Laplace operator. However, as pointed out in [2, 13] more general elliptic operators can be involved. For the purpose of mesh deformation, an elasticity operator may be the ideal choice. In numerical computations, its inverse, the Neumann to Dirichlet map or Poincaré-Steklov is also of importance. Therefore, we first define these operators.

In the sequel, we use the continuous generalized trace map

 γ:H10(Ω,Rd)→H1/2(Γint,Rd)×H−1/2(Γint,Rd),U↦(γ0Uγ1U):=⎛⎜ ⎜⎝UΓint∂nUΓint⎞⎟ ⎟⎠. (36)

Analogously to [13], we define for vector fields with or , the Neumann solution operator for the inner boundary derived from a symmetric and coercive bilinear form

 a:H10(Ω,Rd)×H10(Ω,Rd)→R (37)

by

 Misplaced & (38)

with defined as the solution of the variational problem

 a(U,V)=∫Γintu⊤(γ0V) ds, ∀ V∈H10(Ω,Rd) (39)

where we note that the integral in the right hand side of equation (39) is to be understood as the duality pairing. Furthermore, we define the Dirichlet solution operator for the inner boundary by

 ED:H1/2(Γint,Rd)→H10(Ω,Rd),u↦U (40)

with defined as the solution of the variational problem

 Missing dimension or its units for \rule (41)

Now, we can define the Dirichlet to Neumann map and the Neumann to Dirichlet map as done in the following definition:

{definition}

In the setting above, the Dirichlet to Neumann map and the Neumann to Dirichlet map are defined by

 T :=γ1∘ED:H1/2(Γint,Rd)→H−1/2(Γint,Rd), (42) S :=γ0∘EN:H−1/2(Γint,Rd)→H1/2(Γint,Rd) (43)

where are given in (36).

In obvious generalization of theorem 2.3.1 in [13] from scalar fields to vector fields, we conclude that both operators are symmetric w.r.t. the standard dual pairing , coercive, continuous and that , an observation, for which [2] is cited in [13]. For the purpose of defining an appropriate scalar product on the tangent space of shape manifolds, we define the following mappings. {definition}In the setting above, we define

 η:H(Γint) →H(Γint,Rd) η⊤:H(Γint,Rd) →H(Γint) α ↦α⋅n U ↦n⊤U

where , and thus the projected operators

 Tp :=η⊤∘T∘η:H1/2(Γint)→H−1/2(Γint), (44) Sp :=η⊤∘S∘η:H−1/2(Γint)→H1/2(Γint). (45)

Both operators, and , inherit symmetry, coercivity, continuity and invertibility from the operators . However, we observe in general . Both operators can be used for the definition of a scalar product on the tangent space. In line with the discussion of Sobolev type metrics in [16], we would prefer a scalar product with a smoothing effect like the projected Dirichlet to Neumann map . However, we need its inverse in numerical computations, which is usually not , although spectrally equivalent. We can limit the computational burden, if we use directly as a metric on the tangent space, having a similar smoothing effect but also the advantage of the straight forward inverse . In order to summarize, let us explicitly formulate the operator

 Sp:H−1/2(Γint)→H1/2(Γint),α↦(γ0U)⊤n (46)

where solves the Neumann problem

 a(U,V)=∫Γintα⋅(γ0V)⊤n ds, ∀ V∈H10(Ω,Rd) (47)

which corresponds to an elliptic problem with fixed outer boundary and forces at the inner boundary . Thus, we propose to use the scalar product defined below. {definition} In the setting above, we define the scalar product on by

 gS:H1/2(Γint)×H1/2(Γint)→R,(α,β)↦⟨α,(Sp)−1β⟩=∫Γintα(s)⋅[(Sp)−1β](s) ds. (48)

## 4 Shape quasi-Newton methods based on the metric gS

As already mentioned in section 2 the shape derivative can always be expressed as boundary integral (cf. (4)) due to the Hadamard structure theorem. If , this can be written more concisely as

 DJΓint[V]=∫Γintαf ds. (49)

Due to the isomorphism (32) and the handy expression (49) we can state the connection of with shape calculus, i.e., we can determine a representation of the shape gradient in terms of defined in (48) by

 gS(ϕ,h)=(f,ϕ)L2(Γint),∀ϕ∈C∞(Γint,R), (50)

which is equivalent to

 ∫Γintϕ(s)⋅[(Sp)−1h](s) ds=∫Γintf(s)ϕ(s) ds,∀ϕ∈C∞(Γint,R). (51)

Thus, , where solves

 Extra open brace or missing close brace (52)

which means that the representation of the domain integral formulation in terms of the elliptic form as used in [10] can be – projected to the normal component on – interpreted as the representation of the boundary integral formulation in terms of . However, in both points of view, the information of the shape derivative is in physical terms used as a force (in the domain or on the boundary) and we obtain a vector field as an (intermediate) result, which can serve as a deformation of the computational mesh – identical to Dirichlet deformation.

###### Remark 2

In general, is not necessarily an element of because it is not ensured that is . Under special assumptions depending on the coefficients of a second-order partial differential operator, the right hand-side of a PDE and the domain on which the PDE is defined, a weak solution of a PDE is by the theorem of infinite differentiability up to the boundary [9, theorem 6, section 6.3].

Now, we rephrase the l-BFGS-quasi-Newton method for shape optimization from [24] in terms of the metric and in generalization to domain formulations of the shape derivative. We note that the complete deformation of a shape optimization algorithms is just the (linear) sum of all iterations, which means that the BFGS update formulas can be rephrased directly in terms of the deformation vector field, rather than only as boundary deformations to be transferred to the domain mesh in each iteration.

BFGS update formulas need the evaluation of scalar products, where at least one argument is a gradient-type vector. According to the metric introduced in section 3, we can assume that a gradient type vector can be written as

 u=(γ0U)⊤n (53)

for some vector field . The other argument is either of gradient-type or deformation-type, which can also be assumed of being of the form (53), i.e.,

 v=(γ0V)⊤n (54)

for some . If is a gradient of a shape objective , we observe

 gS(u,v)=DJΓint[V]=DJΩ[V]=a(U,V). (55)

This observation can be used to reformulate the scalar product on the boundary equivalently as for domain representations. In the sequel, we consider only domain representations of , mesh deformations and differences where denotes the vector transport as in [24].

With this notation we formulate the double-loop of an l-BFGS quasi-Newton method:

for  do

end for
for  do

end for
return

The resulting vector is simultaneously a shape deformation as well as a deformation of the domain mesh.

## 5 Numerical results and implementation details

We compare the limited memory BFGS shape optimization algorithms of [24] with the analogous algorithm based on the Riemannian metric , introduced above. We use a test case within the domain , which contains a compact and closed subset with smooth boundary. The parameter is valid in the exterior and the parameter is valid in the interior . First, we build artificial data , by solving the state equation for the setting with . Afterwards, we choose another initial domain and . Figure 1 illustrates the interior boundary around the initial domain and the target domain . The reason for this choice of artificial test data is that we obtain a representation of that can be evaluated at arbitrary points in space since it is represented in finite element basis. Moreover, this construction guarantees that the optimization converges to a reasonable shape that is within the boundaries and not too different to the initial shape such that the mesh remains feasible under deformations.

###### Remark 3

Choosing the measurements as the solution of the model equation (7-11) we obtain that as we assumed in section 2.1.

In the particular test case, which is studied in this section and can be seen in figure 4, the diffusion coefficients are chosen to be and . Further, the initial condition is for all , in and the final time of the simulation is . The results shown in this section are computed under a mild perimeter regularization with , where we did not notice any numerical difference with the case . Yet, for the non-smooth initial configuration shown in figure 5 a stronger regularization has to be chosen in the first iterations as . In this particular case the regularization is controlled by a decreasing sequence from to .

The numerical solution of the boundary value problem (7-10) is obtained by discretizing its weak formulation (12) with linear finite elements in space and an implicit Euler scheme in time. For the time discretization 30 time steps are chosen, which are equidistantly distributed. The diffusion parameter is discretized as a piecewise constant function in contrast to the continuous trial and test functions. This choice of function spaces ensures that the transmission conditions (11) are automatically fulfilled. The corresponding adjoint problem (18-25) can be discretized in the same way. More precisely, it is not necessary to assemble different linear operators, which is attractive in terms of computational effort. All arising linear systems are then solved using the preconditioned conjugate gradient method.

An essential part of a shape optimization algorithm is to update the finite element mesh after each iteration. For this purpose, we use a solution of the linear elasticity equation

 div(σ) =felasinΩ (56) U =0onΓout σ :=λTr(ϵ)I+2μϵ ϵ :=12(∇U+∇UT)

where and are the strain and stress tensor, respectively. Here and denote the Lamé parameters, which can be expressed in terms of Young’s modulus and Poisson’s ratio as

 λ=νE(1+ν)(1−2ν),μ=E2(1+ν). (57)

The solution is then added to the coordinates of the finite element nodes. The Lamé parameters do not need to have a physical meaning here. It is rather essential to understand their effect on the mesh deformation. states the stiffness of the material, which enables to control the step size for the shape update. gives the ratio how much the mesh expands in the remaining coordinate directions when compressed in one particular direction. The numerical results in this work are obtained using and .

Equation (56) is modified according to the optimization approach under consideration. In case we use the surface formulation of the shape derivative (28), the following Dirichlet condition is added on the variable boundary

 U=UsurfonΓint (58)

where is the representation of the shape gradient with respect to the Sobolev metric given in (33). The source term is then set to zero. Otherwise, when the mesh deformation operator is also used as shape metric, assembled according to (27) and there is no Dirichlet condition on . This only covers the portion of the shape derivative for which a volume formulation is available. Parts of the objective function leading only to surface expressions, such as, for instance, the perimeter regularization , are incorporated as Neumann boundary conditions given by

 ∂U∂n=fsurfonΓint. (59)

In the notation of section 3.2 we set as the weak form of the linear elasticity equation leading to

 a(U,V)=∫Ωσ(U):ϵ(V)dx. (60)

For our model problem given in section 2 we have to solve in the context of the domain formulation of the shape derivative and its representation in terms of

 a(U,V)=DjΩ[V]+Djreg(Ω)[V], ∀ V∈H10(Ω,Rd) (61)

where the right hand side is given by the left formulation in (29), in particular in (27) and in (30). Note that (61) is justified by the main result (52) of the previous sections stating that the connection between the volume formulation of shape derivatives and a bilinear form leads to a representation of the shape gradient with respect to .

Both approaches (A versus B below) follow roughly the same steps with a major difference in the way the shape sensitivity is incorporated into the mesh deformation. The appraoch A (domain formulation) is clearly to be preferred because of its implementational ease and less computational effort, if a technical detail discussed below is taken into account. One optimization iteration can be summarized as follows:

1. The measured data has to be interpolated to the current iterated mesh and the corresponding finite element space. Here, this consists of the interpolation between two finite element spaces on non-matching grids.

2. The state and the adjoint equation are solved.

3. Assembly of the linear elasticity equation.
A) Domain formulation:

• The volume form of the shape derivative is assembled into a source term for the linear elasticity mesh deformation. Only test functions whose support includes are considered, which is justified in the subsequent discussion. The behavior of the algorithm with full assembly for all test functions is illustrated in figure 3. Here, the magnitude of the unmodified discretization of the source term is visualized, which shows not only non-zero values outside of due to discretization errors, but leads also to detrimental mesh deformations.

• Shape derivative contributions, which are only available in surface formulation, such as the perimeter regularization, are assembled into the right hand side in form of Neumann boundary conditions.

B) Surface formulation:

• The preliminary gradient given in (28) is evaluated at .

• The -projection of into the space of piecewise linear, continuous functions is conducted. Let this be denoted by .

• Finally the contributions resulting from the regularization, which is here , is added to and we solve the Laplace-Beltrami equation to obtain a representation of the gradient with respect to the Sobolev metric as given in (33).

• then yields the Dirichlet boundary condition (58).

4. Solve linear elasticity equations, apply the resulting deformation to the current finite element mesh and go to the next iteration.

Assembling the right hand side of the discretized weak form (equation (61)) only for test functions whose support intersects with in the volume formulation of step 3 above is due to the following reasoning. In exact integration, the integral should be zero for all test functions which do not have within their support. Thus, nonzero integral contributions are caused by discretization noise, On the other hand, its effect on the optimization algorithm can be understood from a perturbation point of view. We may assume that the Riemannian shape Hessian (where means covariant derivative), whose action in the optimal solution coincides with the action of the shape Hessian, i.e.,

is coercive on the boundary, i.e., for projections , which guarantees a well-posed problem. However, the Hessian operator approximated in the BFGS update strategy described in section 4 deals with a Hessian defined on the whole mesh, which posseses a huge kernel, determined by all vector fields with zero normal component on the boundary. Thus, the space of all admissible deformations has a decomposition

 H10(Ω,Rd)=HΓint⊕H⊥Γint (63)

where and denotes its orthogonal complement in the bilinear form . Shape gradients and increments in lie in only. It is abvious that l-BFGS update formulas produce steps which lie again in only, which means that the optimization algorithm in function spaces acts always on the coercive shape Hessian only. However, the discretized version is a perturbation of the infinite Hessian. Thus, perturbed coercive operators stay coercive, if the perturbation is not too large. But, positive semidefinite operators with a nontrivial kernel, almost inevitably will get directions of negative curvature, when perturbed. These directions of negative curvature will be chosen, if we allow nonzero components in the right hand side of the discretized mesh deformation equation (61) in the interior of the domain. On the other hand, if we do not allow zero components there, the algorithm only acts in the subspace of the discretization of where the projected Hessian is a perturbation of the shape Hessian and thus coercive, if the perturbation is not too large.

We now conclude this section with a discussion of the numerical results. The figures 3 to 5 show the initial configuration and the iterations 2, 4, 20 of the full BFGS algorithm as described in section 4. In figure 3 the algorithm is shown for the unmodified assembly of the right hand side in (61) leading to divergence. Whereas, figure 4 shows a selection of BFGS iterates for the modified source term and with . It is also demonstrated here that the domain-based shape optimization algorithm can be applied to very coarse meshes. This is due to the fact, that there is no dependence on normal vectors like in the case of surface shape derivatives. Finally, figure 2 shows the convergence of l-BFGS with three gradients in memory, full BFGS and the pure gradient method for suface and volume shape derivative formulation, respectively.

In our tests, the convergence with the Laplace-Beltrami representation of the shape gradient seems to require a bit less iterations compared to the domain-based formulation. Yet, the domain-based form is computationally more attractive since it also works for much coarser discretizations. This can be seen in figure 4 where 4 shows the necessary fineness of the mesh for the surface derivative to lead to a reasonable convergence. The coarse grid in 4, however, only works for the domain-based formulation.

Since the volume term in approach A only has to be computed for discretization elements adjacent to it is computationally not more expensive than the surface formulation in approach B. Moreover, the computing time for -projection and solution of tangential Laplace equation is saved. Yet, the BFGS update algorithm in approach A is more expensive then the one in B due to the higher dimension of the involved matrix, which does not play a decisive role since we only store a few gradients in memory. These differences should yet not be overrated. The most expensive operation is the computation of the mesh deformation involving the solution of the linear elasticity equation in both approaches A and B making them comparable in terms of computational costs.

This changes for highly parallel application on supercomputers as investigated in [17]. Operations, which are only performed on surfaces, can drastically affect the scalability of the overall algorithm if the computational load is not balanced also with respect to surface elements. The higher demand for memory of the domain-based formulation seems also not be dramatic since the numerical tests suggest that very few gradients in memory are sufficient for good performance of the l-BFGS method (see figure 2).

## 6 Towards a novel shape space

The scalar product introduced above in section 3 connects shape gradients with deformations. These deformations evaluated at a prior shape give deformed shapes of class , if the deformations are injective and continuous. In the following, it is clarified what we mean by -shapes. The investigations done in the previous section are not limited to shapes, i.e., elements of the shape space . Therefore, this section is devoted to an extension of , i.e., to a novel shape space definition, and its connection to shape calculus.

We would like to recall once again that a shape in the sense of the shape space of Peter W. Michor and David Mumford introduced in [16] is given by the image of an embedding from the unit sphere into the Euclidean space . In view of our generalization, it has technical advantages to consider a prior shape as the boundary of a connected and compact subset with