A shape optimization algorithm for cellular composites
AA shape optimization algorithm for cellular composites
Martin Siebenborn ∗ Andreas Vogel † Abstract
We propose and investigate a mesh deformation technique for PDE constrainedshape optimization. Introducing a gradient penalization to the inner product for lin-earized shape spaces, mesh degeneration can be prevented within the optimizationiteration allowing for the scalability of employed solvers. We illustrate the approachby a shape optimization for cellular composites with respect to linear elastic energyunder tension. The influence of the gradient penalization is evaluated and the parallelscalability of the approach demonstrated employing a geometric multigrid solver onhierarchically distributed meshes.
Keywords: shape optimization, multigrid methods, algorithmic scalability
Mathematics Subject Classification (2010):
Partial differential equation (PDE) constrained shape optimization (e.g. [10, 7, 30]) aims at findingoptimal structural designs for models described by PDEs. In contrast to classical optimal control,where model parameters are represented by distributed functions, the underlying assumption ofshape optimization is that there is a fixed number of available parameters (e.g. materials) whichare spatially separated by interfaces. The variable in this optimization technique is then the shapeof these contours while the topology is assumed to be known and kept fixed.Nowadays, shape optimization has reached a high level of maturity, which can be seen in thefield of applications ranging from, e.g., aerodynamics [19, 13, 9] and acoustics [25] over electrostatics[8] to solid mechanics [3]. The advantage over classical optimal control is a significant reduction ofdegrees of freedom allowing for high resolutions in regions of interest such as the interfaces betweenmaterials. Moreover, the choice of function spaces for an optimal control typically introduces acertain level of smoothness to the parameter, which prevents desired sharp interface solutions.The price one has to pay is to provide initial knowledge about the number of parameters, theirapproximate location, and topology. Typically, shapes are numerically encoded by edges (2d) orfaces (3d) in a finite element mesh. This enables true integer solutions with sharp interfaces, butthe initial topology can not be changed such that it is not straightforward to join or split shapes.For the special class of interface identification problems there are two prominent, alternativeapproaches to shape optimization: first, interfaces can be implicitly resolved by level-set (cf. [28])or phase-field functions (see for instance [22], where a similar problem to the one in this paper isstudied). Second, there is the topology optimization approach, where the binary variable describingmaterial and void is relaxed and a combination of continuous optimization and filter techniques areapplied [5]. Similarly to our considerations in this paper, there is also ongoing effort in maintainingmaterial thicknesses in the context of level-set methods (see for instance [4]).The focus of this article are applications where the shape topology should intentionally beretained. A representative example is a biological cell composite where cells may deform and moverelative to each other, but the topology of cells remains unchanged. E.g., the cell layout of theepidermis (the outermost layer of the human skin) consists of different (and differently shaped)cells surrounded by thin lipid layers. We refer the reader to, e.g., [14, 23, 20, 16] for a reviewof modeling aspects. For such a system, a reasonable question is to ask for an optimized cell ∗ Universit¨at Hamburg, Bundesstraße 55, 20146 Hamburg, [email protected] † Ruhr-Universit¨at Bochum, Universit¨atsstraße 150, 44801 Bochum, [email protected] a r X i v : . [ m a t h . O C ] A p r top Γ bottom Γ side Ω out Ω cell , Ω cell ,N Γ side Γ int Figure 1: Schematic view of the domain Ω consisting of an bulk material Ω out and cellinclusions Ω cell , , . . . , Ω cell ,N with different material properties arrangement with respect to elastic energy under outer force application. Using shape optimizationto analyze such objectives, a certain numerical difficulty is due to the thin layers between cells thatmust be treated appropriately and we develop mathematical tools to address such settings.In order to optimize shapes it is obligatory to define a set of admissible shapes and a corre-sponding tangent space for the descent directions. The probably most prominent approach is todescribe shapes as points on a Riemannian manifold [17, 18]. A lot of attention has been paidto finite element approximations of this manifold and the corresponding tangent spaces (see forinstance [8, 27, 26, 12]).In this paper, we study a modification of inner products in the linearized FEM approximationof shape spaces. We propose a formulation that adapts the modeling of biological tissues to someextend. During shape updates in the optimization process, descent direction with small gradientsshould be promoted, since in growth processes of biological cellular structures the material canonly be compressed to some extend and not completely vanish. This can be modeled by limitinggradients of shape updates. A resulting feature of limited deformation gradients is the preservationof the FEM mesh quality, which is demonstrated within this work. In previous work [29] it has beenobserved that the mesh quality may reduce dramatically during shape updates with the consequenceof unreliable finite element approximations. We also demonstrate that the inner product in shapespace can be implemented using standard finite element techniques, which is our foundation forscalable algorithms on parallel supercomputers.The main contributions of this work thus are: • Optimization algorithm with gradient penalization for cellular composites • Numerical study of the influence of the gradient penalization on optimization results, meshquality, and solver behavior for a cell composite optimization • Parallel scalability analysis for the presented implementationThis paper is structured as follows: In Section 2 the PDE model is introduced and the shapeoptimization problem is formulated. Section 3 describes a gradient based optimization algorithmand the novel inner product is proposed. Numerical experiments and parallel performance of ourimplementation are then investigated in Section 4. Finally, Section 5 discusses the achieved resultsand provides a brief outlook on future research.
The general idea of our model can be illustrated by a geometrical setting as sketched in Figure 1: weconsider a composite material that is fixed at its bottom surface and forces act on the top surface. e assume that the domain of interest consists of several material inclusions with different elasticbehavior. This setting is inspired by a composite of biological cells with some surrounding bulkmaterial. The goal is then to improve the stiffness of the entire object by changing its interiorgeometrical configuration, i.e. relocating and reshaping the inclusions. As an important constraint,the outer shape of the block has to remain unchanged in order to prevent degenerated solutionssuch as removing the force by a vanishing upper surface. If the bulk material and each of theinclusions consist of a homogeneous material, the stiffness of the composite block is a function ofthe shapes of the inclusions, and our aim is to devise an algorithm for the efficient computation ofan optimized shape of these cells. In the following, we provide a precise mathematical formulationfor this setting. Let Ω ⊂ R d and Ω out , Ω cell , , . . . , Ω cell ,N ⊂ Ω ( N ∈ N ) be bounded Lipschitz domains in a d -dimensional space for d ∈ { , } . We assume that the entire domain is made up by the bulkmaterial Ω out and the cell inclusions Ω cell ,i , 1 ≤ i ≤ N , and that cell inclusion do not touch eachother, i.e. we assumeΩ out = Ω \ (cid:32) N (cid:91) i =1 cl(Ω cell ,i ) (cid:33) and cl(Ω cell ,i ) ∩ cl(Ω cell ,j ) = ∅ for i (cid:54) = j, (1)where cl(Ω cell ,i ) denotes the closure of a set. This implies that there is always material betweentwo cells Ω cell ,i and Ω cell ,j . The boundary of the cells is denoted by Γ int := ∪ Ni =1 ∂ Ω cell ,i . Theouter boundary of Ω is subdivided into the disjoint parts Γ top , Γ bottom , Γ side , each with a positiveLebesgue measure. An example of this situation is depicted in Figure 1.The mathematical model is given by the equations of linear elasticity for a displacement field u : R d → R d as div( σ ( u )) = 0 , in Ω ,σ ( u ) · n = f, on Γ top ,σ ( u ) · n = 0 , on Γ side ,u = 0 , on Γ bottom , (2)where σ ( u ) = λ tr( ∇ u ) I + µ ( ∇ u + ∇ u T ) is the stress tensor in terms of the Lam´e parameter λ and µ , and f : R d → R d is a surface force.Defining a function space V := (cid:8) u ∈ H (Ω , R d ) : u = 0 a.e. on Γ bottom (cid:9) on Ω we obtain theweak formulation of problem (2) as follows: Find u ∈ V such that (cid:90) Ω σ ( u ) : ∇ w dx = (cid:90) Γ top f · w ds ∀ w ∈ V . (3)The purpose of this model is the simulation of the elastic behavior of a composite material. Weassume that Ω out is homogeneously filled with one material. In addition, each of the N inclusionsΩ cell ,N is filled with a homogeneous material, but the material may differ between inclusions. Thissituation is reflected in the mathematical model by piece-wise constant Lam´e parameter λ, µ ∈ L ∞ (Ω). For the above settings, we model the stiffness of the material as a function of the inner geometricalconfiguration. We thus consider a minimization problem asmin Γ int J elast ( u, Ω) := 12 (cid:90) Ω σ ( u ) : ∇ u dx s.t. (cid:90) Ω σ ( u ) : ∇ w dx = (cid:90) Γ top f · w ds ∀ w ∈ V . (4)This minimization problem is to be understood in a formal way, i.e. we optimize over the inte-rior shape Γ int of the domain Ω. Note that a change of the location and shape of the inclusionsΩ cell , , . . . , Ω cell ,N changes the spatial distribution of the Lam´e parameter and thereby affects theobjective function in (4). n order to represent shape modifications we concentrate on the so-called perturbation of iden-tity. This means that, in terms of the initial shape Ω = Ω , we considered for t ≥ t = {F t ( x ) : x ∈ Ω } with F t : Ω → Ω t , F t ( x ) = x + tv ( x ) (5)for a given displacement vector field v ∈ W , ∞ (Ω , R d ). Note that this choice guarantees that theouter shape of Ω is retained. Moreover, for t sufficiently small, Ω t remains admissible in the sensethat there are no overlappings of the cell inclusions.In order to investigate sensitivities of the objective J with respect to the shape we mainly followthe concepts introduced in [30] and [7]. For a vector field v as defined above, the directional shapederivative of the functional J evaluated at Ω in direction v is defined by dJ (Ω)[ v ] = lim t (cid:38) J (Ω t ) − J (Ω) t . (6)Note that the general definition above implicitly takes into account the dependence of J on thestate variable u . Since u is the solution of the elliptic equation (3) it depends on the interior shapeof the domain Ω t , in particular on the distribution of Lam´e parameters. Therefore, J ( u t , Ω t ) hasto be considered in (6) and derivatives with respect to t have to be computed.We consider the associated Lagrangian to problem (4) given by L ( u, p, Ω) = 12 (cid:90) Ω σ ( u ) : ∇ u dx − (cid:90) Ω σ ( u ) : ∇ p dx + (cid:90) Γ top fp dx. (7)The state variable u and the corresponding adjoint p can then be characterized as the saddle pointof (7) as J ( u, Ω) = min u ∈ H (Ω , R d ) max p ∈ H (Ω , R d ) L ( u, p, Ω) . (8)In order to obtain the shape derivative dJ (Ω)[ v ], the equation above is reformulated in terms of Ω t , u t and p t . Following the presentations in [7, Chapter 10], where an application of the Correa andSeeger theorem [6, Theorem 2.1] to shape sensitivity analysis is described, the differentiability ofthe right hand side in (8) with respect to t is obtained. From these results it follows that dJ (Ω)[ v ] = ddt L ( u t , p t , Ω t ) | t =0 . (9)First we observe that in a saddle point for t = 0 it holds L u ( u, p, Ω) = 0. This yields thecondition p = u for the adjoint variable, which we can therefore substitute in the following. Let u t := u t ◦ F t be the pull-back of u t onto the domain Ω and observe that ∇ u t = ( D F t ) − T ∇ u t . Inaddition, due to the special choice of F t it holds D F t = I + t ∇ v, ddt det( D F t ) (cid:12)(cid:12)(cid:12)(cid:12) t =0 = div( v ) , ddt ( D F t ) − (cid:12)(cid:12)(cid:12)(cid:12) t =0 = − Dv. (10)With these identities we can compute ddt L ( u t , p t , Ω t ) | t =0 = ddt (cid:90) Ω t σ ( u t ) : ∇ u t dx (cid:12)(cid:12)(cid:12)(cid:12) t =0 = ddt (cid:90) Ω t (cid:16)(cid:104) ( λ ◦ F t − ) tr( ∇ u t ) I + ( µ ◦ F t − ) (cid:16) ∇ u t + ∇ u Tt (cid:17)(cid:105) : ∇ u t (cid:17) dx (cid:12)(cid:12)(cid:12)(cid:12) t =0 = ddt (cid:90) Ω (cid:104) λ tr(( D F t ) − T ∇ u t ) I + µ (cid:16) ( D F t ) − T ∇ u t + (( D F t ) − T ∇ u t ) T (cid:17)(cid:105) : (( D F t ) − T ∇ u t ) det( D F t ) dx (cid:12)(cid:12)(cid:12) t =0 = 12 (cid:90) Ω (cid:104) λ tr(( − Dv ) T ∇ u ) I + µ (cid:16) ( − Dv ) T ∇ u + (( − Dv ) T ∇ u ) T (cid:17)(cid:105) : ∇ u + σ ( u ) : (cid:16) ( − Dv ) T ∇ u (cid:17) + div( v ) ( σ ( u ) : ∇ u ) dx e finally obtain the directional shape derivative dJ elast ( u, Ω)[ v ] = (cid:90) Ω
12 div( v ) ( σ ( u ) : ∇ u ) − σ ( u ) : ( ∇ v ∇ u ) dx = (cid:90) Ω (cid:18)
12 ( σ ( u ) : ∇ u ) I − ( ∇ u ) T σ ( u ) (cid:19) : ∇ v dx, (11)where we exploited the special structure of the stress tensor for isotropic linear elastic material.This expression can also be found in [21] with a slightly different derivation. Note that this formulais different to the expression used in the previous work [29] where we found numerical artifactsinfluencing the optimization. In this paper, we use the above expression and modify the algorithmin order to address these issues.Maintaining the outer shape of Ω is not yet incorporated into the model. Quite the contrary: thesensitivities dJ (Ω)[ v ] might show largest absolute values close to the boundaries due to the fact thatshrinking Γ top or extending Γ bottom typically leads to a steep descent of the objective. To avoid thisunintended behavior, we thus define a domain D ⊂ Ω, slightly smaller than Ω, which leaves out theouter boundaries. The directional derivatives dJ (Ω)[ v ] are then only computed within D . Note thatthis procedure is equivalent to considering deformation fields v with support in D only. In addition,from a practical aspect, one typically encounters nonzero values in the discrete approximation of dJ (Ω)[ v ] even for finite element shape functions that do not intersect Γ int . However, a FE shapefunction φ : D → R d can not contribute to a descent direction v if supp( φ ) ∩ Γ int = ∅ . Since Lam´eparameter are constant for each cell and the bulk material, only deformations of Γ int can affect theobjective functional. We thus add an additional step to the optimization algorithm (cf. Alg. 1),which sets these shape derivative values to zero in order to avoid unintended deformations. Our problem setting is motivated by an elasticity model for biological cell composites such as thehuman skin. For such a system, also geometrical optimization aspects play a role. We thereforemodify the objective function by two weighted additional terms using ν elast > ν vol , ν peri ≥ Γ int J ( u, Ω) := ν elast J elast ( u, Ω) + ν vol J vol (Ω) + ν peri J peri (Ω)= ν elast (cid:90) Ω σ ( u ) : ∇ u dx + ν vol (cid:90) Ω out dx + ν peri (cid:90) Γ int ds s.t. (cid:90) Ω σ ( u ) : ∇ w dx = (cid:90) Γ top f · w ds ∀ w ∈ V . (12)The additional terms represent the minimization of the volume of the bulk material in Ω out and theminimization of the perimeter (2d) or surface (3d) of the inclusions Γ int . The motivation for theseterms is to find a space-filling design with minimal cell perimeter/surfaces. For a derivation of theshape derivatives in direction v of the geometrical quantities, given by dJ vol (Ω)[ v ] = (cid:90) Γ int v · n ds = − (cid:90) Ω out div( v ) dx, (13) dJ peri (Ω)[ v ] = (cid:90) Γ int div( v ) − ∂v∂n · n ds, (14)we refer the reader to [30]. In the equation above, n is assumed to be the normal vector field atΓ int , which points outwards of the inclusions Ω cell , , . . . , Ω cell ,N . The equations (11), (13), and (14) state the sensitivities of the objective in terms of variationsaccording to the perturbation of identity. Since our aim is a gradient descent method, a space ofadmissible shapes and an according inner product has to be determined. The shape derivatives canthen be represented as gradients with respect to this inner product and used as shape updates inthe sense of a descent direction. t is a common approach (see, e.g., [8, 26]) to locally approximate the shape space by H (Ω , R d )-deformations to the initial shape Ω and equip this with the bilinear form originating from an ellipticPDE. Here, we focus on linear elasticity, i.e. we choose an inner product on H (Ω , R d ) by a ( v, w ) := (cid:90) Ω σ ( v ) : ∇ w dx. (15)Let dJ (Ω)[ w ] be the shape derivative evaluated in Ω in direction w . In order to obtain thedeformation field v ∈ H (Ω , R d ) we then solve a ( v, w ) = dJ (Ω)[ w ] ∀ w ∈ H (Ω , R d ) . (16)This step is the representation of the derivative dJ (Ω)[ · ] with respect to the inner product a ( · , · ) in H (Ω , R d ). Note that deformations in H (Ω , R d ) are not admissible in general, since a Lipschitzdomain deformed via F t in terms of v may loose this property.A consequence of a gradient descent method applied to a situation as the one depicted in Figure1 is that overlappings are likely to appear during the optimization. Moreover, since the domain Ωand the inclusions Ω cell , , . . . , Ω cell ,N are represented by edges and faces in a finite element mesh,the topology can not be changed. It is thus of interest to modify a ( · , · ) such that large gradients inthe resulting deformation field are avoided. This has two effects: first, the quality of the FEM meshis retained to some extend even under large deformations. This is crucial for the applied linearsolver and the overall approximation quality of the FEM method. Second, we can move modelingaspects into a ( · , · ). The ultimate goal of our method is to reflect the behavior of biological tissue inthe human skin. Selecting descent directions which avoid large gradients promotes shape updateswhere the thickness of channels between cells does not tend to zero too fast.Recall that (16) is the necessary and sufficient optimality condition of the minimization problemmin v ∈ H (Ω , R d ) j ( v ) := 12 a ( v, v ) − dJ (Ω)[ v ] . (17)Our aim is to modify the selected descent directions to prefer translations over compression of thematerial in Ω. We thus enrich (17) by a penalty term involving the gradients of the shape update.We thus propose the following modified formulationmin v ∈ H (Ω , R d ) j ( v ) := j ( v ) + j ( v ) := j ( v ) + ν penalty (cid:90) Ω (cid:16)(cid:0) ∇ v : ∇ v − b (cid:1) + (cid:17) dx (18)for a bound b > ν penalty > · ) + denotes the positive part, i.e. max { , (cid:107)∇ v (cid:107) F − b } with (cid:107)∇ v (cid:107) F = ∇ v : ∇ v . For largepenalty parameter ν penalty , gradients exceeding a certain threshold are thus avoided. A generalizedderivative evaluated in a direction w ∈ H (Ω , R d ) follows as dj ( v )[ w ] = ν penalty (cid:90) Ω (cid:0) ∇ v : ∇ v − b (cid:1) + ∇ v : ∇ w dx. (19)We can then again formulate the necessary optimality conditions for (18) and obtain the variationalformulation: Find v ∈ H (Ω , R d ) such that g ( v, w ) := a ( v, w ) + dj ( v )[ w ] = dJ (Ω)[ w ] for all w ∈ H (Ω , R d ) . (20)Note that g ( · , · ) is nonlinear and not even differentiable. The left hand side g ( · , · ) only constitutesa suitable inner product in H (Ω , R d ), if the term involving the positive part is inactive, i.e. when ∇ v : ∇ v − b ≤ v exceed the threshold b , the formulation (20)switches to a nonlinear behavior and can no longer be represented by a symmetric, coercive bilinearform such as a ( · , · ).Finally, we solve (20) for v via a semi-smooth Newton method and use the derivative in directions w, h ∈ H (Ω , R d ) given by j (cid:48)(cid:48) ( v )[ w, h ] = (cid:90) Ω σ ( h ) : ∇ w dx + ν penalty (cid:90) Ω γ (cid:2) ∇ v : ∇ h ) · ( ∇ v : ∇ w ) + (cid:0) ∇ v : ∇ v − b (cid:1) ∇ h : ∇ w (cid:3) dx, (21) lgorithm 1 Gradient-penalized optimization algorithm Choose initial domain Ω ; choose ν elast , ν vol , ν peri ; choose ν penalty , b, t for k = 0 , , , . . . do Compute displacement field u k , i.e. solve model equations (3), for the current Ω k Compute shape derivative as given by the equations (11), (13), (14), i.e. assemble dJ (Ω k ) = ν elast dJ elast ( u k , Ω k ) + ν vol dJ vol (Ω k ) + ν peri dJ peri (Ω k ) Reset shape derivative away from surface, i.e. set dJ (Ω k )[ φ ] = 0 if supp( φ ) ∩ Γ int = ∅ Compute descent direction v , i.e. solve g ( v, w ) = dJ (Ω k )[ w ] for all w ∈ H (Ω , R d ) Update shape: Ω k +1 = { x + tv ( x ) : x ∈ Ω k } end for where L ∞ (Ω) (cid:51) γ = (cid:40) , if (cid:107)∇ v (cid:107) F − b > , else . (22)For theoretical background on Newton’s methods for semi-smooth functions we refer the readerto [31]. Furthermore, the threshold on gradients of state variables is inspired by the theory ofoptimal control with gradient constraints (see for instance [11]). Algorithm 1 shows the resultingoverall algorithm.Since the outer shape of the domain Ω should not be modified, we restricted ourselves to a H (Ω , R d )-approximation of descent directions in the shape space. However, for the numericalresults in the next section, we slightly change the space for the gradient representation: we wantto allow tangential movement of nodes along the boundary to allow for better mesh quality. In thetwo dimensional setting of Figure 1 it is thus only required that the x -component at Γ side and the y -component at Γ bottom ∪ Γ top is zero. Note that this is sufficient to guarantee an unchanged outerboundary ∂ Ω. We test the proposed method using a setup that is inspired by the topology of the upper human skin(see e.g. [20] and the references therein). We emphasize that the presented numerical results are farfrom being biologically meaningful. In contrast, these simulations are solely used to demonstratethe numerical properties of the method, and we have chosen artificial material parameters anddomain sizes for this example. However, the setup exposes a reasonable similarity with biologicallyobserved patterns and we leave more realistic simulations for future work.For our test case, the topology of the domain Ω = (0 , is shown in the initial configurationin Figure 2. The outer domain (red) contains several cell inclusions which differ in their materialproperties. Cells of the same color have the same material properties and we number the eightcell types from top to bottom ( k = 0 , . . . , λ, µ ) = (1 . , . , on Ω out , (23)for the outer space material, and the different cell types Ω cell,k , k = 0 , . . . , λ, µ ) = ( λ k , µ k ) , on Ω cell,k , with (24) λ k = (1 − k · λ top + k · λ bottom , λ top = 1 . , λ bottom = 2 . µ k = (1 − k · µ top + k · µ bottom , µ top = 0 . , µ bottom = 0 . f = 0 .
1, at the top boundary Γ top , and usethe plain strain assumption. For the computation of the shape derivate, we use globally constantvalues ( λ, µ ) = (0 . , .
0) on the whole domain Ω. The coarse grid meshing has been performedwith ProMesh [2], visualization is by ParaView [1], and all simulations are carried out with thesimulation toolbox UG4 [32]. . . . . O b j ec t i v e v a l u e Optimization step J ( = ν elast J elast + ν vol J vol + ν peri J peri ) ν elast J elast (Energy) ν vol J vol (Volume Ω out ) ν peri J peri (Perimeter) Figure 3: Value of the objective function for the different components for b = 0 . The results of an optimization for ν elast = 100, ν vol = 1, ν peri = 0 . ν penalty = 5 · , b = 0 . t = 1, and 2 levels of refinement for the FEM mesh is shown in Figure 4. As solver, we used theNewton scheme without linesearch for the nonlinearity, and the linear system of equations are solvedusing the BiCGStab method with a geometric multigrid preconditioner employing 3 Block-Jacobipre- and postsmoothing steps and a LU base solver.The simulation results depict the change of the surface shape for selected optimization steps.As seen from the reference configuration, the objective first quickly increases all cells in order toreduce the outer volume, where the penalty term helps to still keep intercell channels open duringthe optimization. Once the domain is almost completely filled with cells, the optimization turns toan increase of the stiffer cells at the bottom and decreases the volume of the less stiff cells at thetop. From the displacement figures, we can observe that the overall displacement diminishes withthe optimization which reflects the fact that now the material has become stiffer due to the newsize and location of the cell domains.This behavior is also reflected in the behavior of the objective functional and we show the valuesof the individual components in Figure 3. As can be observed, the computed shape movement isa descent direction for the objective value. Once the outer volume contribution has been removed onsiderably, the energy contribution is balanced against the surface term by decreasing the size ofless stiff cells. In Figure 5, we show the influence of the penalization of the norm of the gradient with respectto the parameter b . With increasing b , i.e. allowing larger deformations within each optimizationstep, the cells are inflated more rapidly. The boundaries of the different shape surfaces thereforeapproach quicker with two effects: first, the channels between the cells become very close and themesh is thereby pretty distorted. As a consequence, the optimization stops already after a few (25for b = 0 .
1, 3 for b = 1 .
0) steps due to a failure of the iterative solver. Choosing a direct solver likeLU factorization in such a case can clearly enforce the computation of a solution, however, it simplyresults in a penetration of the FEM mesh elements and therefore to meaningless results. Second,due to the oversized (i.e. non-penalized) shape gradient, the optimization does not enter the regimewhere cell stiffnesses are balanced, but only remains in the inflation part of the optimization.
For the numerical solution in this example, we have chosen a Newton iteration to resolve thenonlinearity for shape gradient computation which is stopped at an absolute or relative reductionof 10 − (whatever applies first). Each linearization within a Newton step is solved by the geometricmultigrid preconditioned BiCGstab method down to a relative residuum reduction of 10 − whichempirically shows to be sufficient to maintain the quadratic reduction rate of the Newton method.The iteration count for the Newton method is restricted to 200 steps and otherwise the simulation isaborted. For the linear elasticity problem, we apply the geometric multigrid accelerated BiCGStabsolver directly and solve down to an absolute or relative reduction of the order 10 − while imposinga maximum iteration count of 2000 steps.The solver behavior for each optimization step is depicted in Figure 6. For the individual op-timization steps and different values of b , We show the iteration counts for the Newton solver, theaverage iteration count of the linear solver within a Newton step, and the linear solver iterationcount for the elasticity problem. A stopped curve corresponds to a simulation run where the conver-gence criteria could not be met within the imposed maximal iteration count and the computationwas aborted. Generally two observations can be stated: first, with larger values of b (allowinglarger deformations for one optimization step), the iteration count for the nonlinear as well as thelinear solver quickly increases dramatically and the simulation is stopped. Second, with sufficientpenalization, the simulation can run for a reasonable amount of steps. Again, the iteration countincreases towards the end of the simulation which we attribute to the distorted finite element mesh.However, with penalization and up to 100 optimization steps (see Fig. 4 for a visualization), theiteration count for the Newton method roughly remains constant (approx. 15 steps for each opti-mization step), and the average linear iteration count for the relative reduction within a newtonstep increases from 3 to a moderate number of 40 steps. For the linear elasticity solver, the iterationcount increases from 20 for the nicely shaped mesh at the beginning to roughly 75 steps for thedistorted mesh at optimization step 100. For the proposed method, the topology of the finite element mesh is fixed and not altered during thesimulation. This clearly has the drawback that element shapes become unpleasant, and with thedeteriorating mesh quality also a reduction in FE solution quality as well as in the solver efficiencyhas to be expected. This has been observed for the iteration counts in our examples, however, theintroduced penalization helps to mitigate this behavior and to allow for a reasonable optimization.On the other hand, the fixed FE mesh topology has an important advantage: with respect toparallelization, the mesh has not to be rebalanced or redistributed among involved processes. Aninitial, well-balanced partition of the finite element mesh onto the processes will remain for the wholeoptimization range and this reasoning applies also to the load balance of coarse mesh hierarchies.Therefore, using available and established linear solver implementations, a good scaling is to beexpected once a good initial load balancing can be achieved.For this study, we use the parallel geometric multigrid approach with hierarchically distributedmeshes [24]. This approach has shown close to optimal scalability on hundred thousands of cores. b : (a) 0 .
01, (b) 0 .
05, (c) 0 .
1; LU enforced solution (d) with mesh penetration for b = 0 . I t e r a t i o n c o un t li n e a r s o l v e r ( a v g ) Optimization step b = 0 . b = 0 . b = 0 . b = 1 . (a) Average number of linear iteration stepsto reduce the residuum by a factor of 10 − within a Newton step I t e r a t i o n c o un t n o n li n e a r s o l v e r Optimization step b = 0 . b = 0 . b = 0 . b = 1 . (b) Number of steps for the nonlinear itera-tion (Newton method) I t e r a t i o n c o un t( li n e a r e l a s t i c i t y ) Optimization step b = 0 . b = 0 . b = 0 . b = 1 . (c) Number of iterations for the linear solverof the elasticity problem M e s h q u a li t y Optimization step b = 0 . b = 0 . b = 0 . b = 1 . (d) Mesh quality (ratio between circum-sphere and inscribed sphere radius of an el-ement); shown are maximal (solid) and me-dian (dashed) values Figure 6: Iteration counts for linear and nonlinear solver (initial mesh 2 times refined)
Our mesh partitioning is based on the ParMETIS [15] software. Here, we perform a strong andweak scaling study for the optimization components and measure accumulated times for the first10 optimization steps. We use the BiCGStab method preconditioned with the geometric multigridusing a Block-Jacobi smoother with 3 pre- and postsmoothing steps, V-cycle, assembled coarsegrid matrices, and a gathered LU base solver. All measurements have been performed on the
Statik&Dynamik -Cluster at the Department of Civil and Environmental Engineering of the RuhrUniversity Bochum. The Cluster consists of 92 nodes, each with two Intel Xeon Skylake Gold 6148(40 cores per node) and 192GB memory, i.e. a maximum of 3680 cores. Our study tries to fullyutilize the allocated nodes and thus starts at 40 cores and increases the core count by a factor of4 (refining the 2d mesh with each increase for the weak scaling study) up to 2560 cores. In order o account for jitter, we have performed five runs for each measurement and present the arithmeticmean.The results are presented in Figure 7 for the strong scaling, and in Figure 8 for the weak scalingstudy. We present wallclock times (accumulated over 10 optimization steps) and the gained speeduprelative to 40 cores. Both studies show a good scalability and speedup, and the number of Newtonsteps remains almost constant even for the weak scaling study. A surprising fact is the increase inthe iteration count of the linear solver. While this is to be expected for the weak scaling (a slightlydifferent numerical problem is solved with each mesh refinement), this should not be the case forthe strong scaling study. Note that we employ a Jacobi smoother and the mathematical algorithmis thereby exactly the same even for parallel computations. A detailed investigation revealed thatsmall differences in round-off errors for different parallel partitions accumulate for these simulationsand explain the slightly different iteration counts. T i m e [ s ] Number of processes T ass T init T solve Sp ee dup ( s t r o n g ) Number of processes S ass S init S solve S ideal (a) Linear elasticity problem T i m e [ s ] Number of processes T lin,init T lin,solve T newton Sp ee dup ( s t r o n g ) Number of processes S lin,init S lin,solve S newton S ideal (b) Computation of shape derivative Procs Refs DoFs Linear solver Newton solver Linear solver(elasticity) (shape derivative) (shape derivative)40 7 37,153,922 305 140 567160 7 37,153,922 297 140 567640 7 37,153,922 308 140 6272560 7 37,153,922 306 140 630 (c) Accumulated iteration counts for the strong scaling study
Figure 7: Strong scaling: accumulated wallclock time for the first 10 optimization steps (left),speedup relative to 40 cores (right), and accumulated iteration counts (table). Shown arethe fine level matrix assembling ( ass ), multigrid initialization ( init ) and solution ( solve )phase, and the overall newton algorithm ( newton )12 T i m e [ s ] Number of processes T ass T init T solve Sp ee dup ( w e a k ) Number of processes S ass S init S solve S ideal (a) Linear elasticity problem T i m e [ s ] Number of processes T lin,init T lin,solve T newton Sp ee dup ( w e a k ) Number of processes S lin,init S lin,solve S newton S ideal (b) Computation of shape derivative Procs Refs DoFs Linear solver Newton solver Linear solver(elasticity) (shape derivative) (shape derivative)40 6 9,291,330 303 132 544160 7 37,153,922 297 140 567640 8 148,592,898 301 140 6252560 9 594,326,018 293 140 625 (c) Accumulated iteration counts for the weak scaling study
Figure 8: Weak scaling: accumulated wallclock time for the first 10 optimization steps (left),speedup relative to 40 cores (right), and accumulated iteration counts (table). Shown arethe fine level matrix assembling ( ass ), multigrid initialization ( init ) and solution ( solve )phase, and the overall newton algorithm ( newton )13
Conclusion and outlook
The results in this article propose a shape optimization algorithm that serves two purposes: First,towards a modeling of cellular composites, it is possible to maintain the thickness of inter-cellularchannel width by promoting shape descent directions with penalized gradients. We incorporate thisas a modification into well-established inner products in a locally linearized description of shapevariations. Second, the gradient penalization prevents early mesh degeneration. Thereby multigridscalability can be retained, which is a mandatory ingredient in our simulations towards realistic,large-scale simulations of cell structures in the human skin. We want to emphasize that our proposedmethod is not limited to linear elastic models with a limited degree of realism. The optimizationalgorithm immediately applies to more complex, potentially nonlinear elasticity models without anychanges in methodology. Under the assumption of hyperelastic material laws, FE models for thehuman skin are a vivid field of research (e.g., [23, 16]). It is thus our next step to apply the proposedoptimization algorithm to refined models with experimentally validated material parameters.
Acknowledgements
This work has been supported by the DFG Priority Program 1648
Software for Exascale Comput-ing (SPPEXA) in the project
Exasolvers . We gratefully acknowledge the computing time on the
Statik&Dynamik -Cluster at the Department of Civil and Environmental Engineering of the RuhrUniversity Bochum, Germany.
References [1] .[2] .[3] G. Allaire.
Shape optimization by the homogenization method , volume 146. Springer Science &Business Media, 2012.[4] G. Allaire, F. Jouve, and G. Michailidis. Thickness control in structural optimization via alevel set method.
Structural and Multidisciplinary Optimization , 53(6):1349–1382, 2016.[5] M.P. Bendsoe and O. Sigmund.
Topology Optimization: Theory, Methods, and Applications .Springer Berlin Heidelberg, 2011.[6] R. Correa and A. Seeger. Directional derivative of a minmax function.
Nonlinear Analysis:Theory, Methods & Applications , 9(1):13–22, 1985.[7] M.C. Delfour and J.-P. Zol´esio.
Shapes and Geometries: Metrics, Analysis, Differential Cal-culus, and Optimization , volume 22 of
Advances in Design and Control . SIAM, 2nd edition,2001.[8] P. Gangl, A. Laurain, H. Meftahi, and K. Sturm. Shape optimization of an electric motorsubject to nonlinear magnetostatics.
SIAM Journal on Scientific Computing , 37(6):B1002–B1025, 2015.[9] M.B. Giles and N.A. Pierce. An introduction to the adjoint approach to design.
Flow, turbulenceand combustion , 65(3-4):393–415, 2000.[10] J. Haslinger and R.A.E. M¨akinen.
Introduction to Shape Optimization: Theory, Approximation,and Computation , volume 7 of
Advances in Design and Control . SIAM, 2003.[11] M. Hinterm¨uller and K. Kunisch. PDE-constrained optimization subject to pointwise con-straints on the control, the state, and its derivative.
SIAM Journal on Optimization ,20(3):1133–1156, jan 2010.[12] J. A. Iglesias, K. Sturm, and F. Wechsung. Two-dimensional shape optimization with nearlyconformal transformations.
SIAM Journal on Scientific Computing , 40(6):A3807–A3830, jan2018.[13] A. Jameson. Aerodynamic shape optimization using the adjoint method.
Lectures at the VonKarman Institute, Brussels , 2003.[14] O. G. Jepps, Y. Dancik, Y. G. Anissimov, and M. S. Roberts. Modeling the human skinbarrier - towards a better understanding of dermal absorption.
Advanced drug delivery reviews ,65(2):152–168, 2013.
15] George Karypis, Kirk Schloegel, and Vipin Kumar. Parmetis: Parallel graph partitioning andsparse matrix ordering library. 01 1997.[16] Georges Limbert. Mathematical and computational modelling of skin biophysics: a review.
Proceedings. Mathematical, physical, and engineering sciences , 473 2203:20170257, 2017.[17] P.M. Michor and D. Mumford. Vanishing geodesic distance on spaces of submanifolds anddiffeomorphisms.
Documenta Mathematica , 10:217–245, 2005.[18] P.M. Michor and D. Mumford. Riemannian geometries on spaces of plane curves.
Journal ofthe European Mathematical Society , 8(1):1–48, 2006.[19] B. Mohammadi and O. Pironneau.
Applied Shape Optimization for Fluids . Oxford UniversityPress, 2001.[20] A. N¨agel, M. Heisig, and G. Wittum. A comparison of two-and three-dimensional models for thesimulation of the permeability of human stratum corneum.
European Journal of Pharmaceuticsand Biopharmaceutics , 72(2):332–338, 2009.[21] A.A. Novotny, R.A. Feij´oo, E. Taroco, and C. Padra. Topological sensitivity analysis forthree-dimensional linear elasticity problem.
Computer Methods in Applied Mechanics andEngineering , 196(41-44):4354–4364, sep 2007.[22] P. Penzler, M. Rumpf, and B. Wirth. A phase-field model for compliance shape optimization innonlinear elasticity.
ESAIM: Control, Optimisation and Calculus of Variations , 18(1):229–258,2012.[23] B. Querleux.
Computational Biophysics of the Skin . Pan Stanford Publishing, 2014.[24] Sebastian Reiter, Andreas Vogel, Ingo Heppner, Martin Rupp, and Gabriel Wittum. A mas-sively parallel geometric multigrid solver on hierarchically distributed grids.
Comp. Vis. Sci. ,16(4):151–164, 2013.[25] S. Schmidt, E. Wadbro, and M. Berggren. Large-scale three-dimensional acoustic horn opti-mization.
SIAM Journal on Scientific Computing , 38(6):B917–B940, 2016.[26] V. Schulz, M. Siebenborn, and K. Welker. Efficient PDE constrained shape optimization basedon Steklov–Poincar´e-type metrics.
SIAM Journal on Optimization , 26(4):2800–2819, 2016.[27] V.H. Schulz and M. Siebenborn. Computational comparison of surface metrics for PDE con-strained shape optimization.
Computational Methods in Applied Mathematics , 16(3):485–496,2016.[28] James Albert Sethian.
Level set methods and fast marching methods: evolving interfaces incomputational geometry, fluid mechanics, computer vision, and materials science , volume 3.Cambridge university press, 1999.[29] M. Siebenborn and K. Welker. Algorithmic aspects of multigrid methods for optimization inshape spaces.
SIAM Journal on Scientific Computing , 39(6):B1156–B1177, jan 2017.[30] J. Sokolowski and J.-P. Zol´esio.
Introduction to Shape Optimization , volume 16 of
Computa-tional Mathematics . Springer, 1992.[31] M. Ulbrich. Semismooth newton methods for operator equations in function spaces.
SIAMJournal on Optimization , 13(3):805–841, jan 2002.[32] A. Vogel, S. Reiter, M. Rupp, A. N¨agel, and G. Wittum. UG 4: A novel flexible softwaresystem for simulating PDE based models on high performance computers.
Comp. Vis. Sci. ,16(4):165–179, 2013.,16(4):165–179, 2013.