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

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


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.


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


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


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:

(domain formulation) (3)
(boundary formulation) (4)

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


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.



Figure 1: Example of a domain where and denotes the unit outer normal to at

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



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


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


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


The linear form in (12) is given by




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


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


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


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


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


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 :




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


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.,


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


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


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


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




with defined as the solution of the variational problem


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


with defined as the solution of the variational problem


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


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


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

where , and thus the projected operators


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


where solves the Neumann problem


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


4 Shape quasi-Newton methods based on the metric

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


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


which is equivalent to


Thus, , where solves


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


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.,


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


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

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.

Our investigations focus on the comparison between two l-BFGS optimization approaches: The first approach is based on the surface expression of the shape derivative, as intensively described in [24]. Here, a representation of the shape gradient at with respect to the Sobolev metric (33) is computed and applied as a Dirichlet boundary condition in the linear elasticity mesh deformation. This involves two operations, which are non-standard in finite element tools and thus leads to additional coding effort. Since we are dealing with linear finite elements the gradient expressions of state and adjoint in (28) are piecewise constant and can not be applied directly to the mesh as deformations. We thus have to implement a kind of -projection on (cf. [24]) bringing back the sensitivity information into the space of continuous, linear functions. The next additional piece of code is a discrete version of the Laplace-Beltrami operator for the Sobolev metric (33). The essential part of this is the solution of a tangential Laplace equation on the surface . Therefore, we follow the presentations [15] and artificially extend our 2D grid in the third coordinate direction. The second approach, discussed in sections 3 and 4, involves the volume formulation of the shape derivative and a corresponding metric, which is very attractive from a computational point of view. The computation of a representation of the shape gradient with respect to the chosen inner product of the tangent space is now moved into the mesh deformation itself. The elliptic operator (cf. (37)) – here the linear elasticity – is both used as inner product and mesh deformation leading to only one linear system, which has to be solved. Besides saving brain work in the calculation of the shape derivative, a lot of coding work is obsolete using surface formulation of shape derivatives. Moreover, it is not always clear how the surface formulation looks like and which additional assumptions have to be made in its derivation. A discussion of the l-BFGS algorithm used within this algorithm can be seen in section 4.

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


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


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


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


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


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


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


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.

Figure 2: Convergence of the optimization iteration measured as an approximation to the geodesic distance in the shape space on a grid with approx. 100,000 cells



BFGS iterates with unmodified approximation of volume shape derivative


Magnitude of unmodified volumic source term

Figure 3: Wrong mesh deformations and source term due to discretization errors in the unmodified right hand side of (61)



Smooth deformations and convergence to optimal shape due to modified source term {subfigure}1.0


The domain-based optimization approach also enables the use of much coarser spatial discretizations (approx. 1000 cells)

Figure 4: BFGS iterates with corrected source term (27) indicating mesh independent convergence
Figure 5: Smooth mesh deformations even with kinks in the initial configuration due to regularization

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