Mesh quality preserving shape optimization using nonlinear extension operators

In this article, we propose a shape optimization algorithm which is able to handle large deformations while maintaining a high level of mesh quality. Based on the method of mappings we introduce a nonlinear extension operator, which links a boundary control to domain deformations, ensuring admissibility of resulting shapes. The major focus is on comparisons between well-established approaches involving linear-elliptic operators for the extension and the effect of additional nonlinear advection on the set of reachable shapes. It is moreover discussed how the computational complexity of the proposed algorithm can be reduced. The benefit of the nonlinearity in the extension operator is substantiated by several numerical test cases of stationary, incompressible Navier-Stokes flows in 2d and 3d.


Introduction
Shape optimization is a discipline in the field of optimization constrained by partial differential equations (PDEs). Here the contour of the domain Ω, where typically a PDE models the effects of interest, plays the role of the optimization variable. Possible variants are that the outer shape of Ω is to be determined, e.g. when Ω represents a solid body, or interior interfaces, which separate spatially discontinuous coefficients such as material properties. Shape optimization in general is nowadays an active field of research with applications ranging from electrostatics [22], interface identification in transmission processes [27,10,21], fluid-dynamics [25,2,7], acoustics [31], image restoration and segmentation [12] and composite material identification [28,23] to nano-optics [13].
In this article, we focus on shape optimization in fluid dynamics, which is also one of the pioneering applications in this field [19,15,9]. In general, the optimization problem can be formulated as where j is a shape functional depending on a state variable y and the shape of the domain Ω. Moreover, y fulfills the PDE constraint E, which itself depends on Ω. A typical example is an obstacle specimen Ω obs in a flow tunnel Ω as depicted in figure 1. One of the main questions is an appropriate choice of the set of admissible shapes G adm , in which optimization takes place. For problems of this type, two prominent approaches can be identified in the literature. On the one hand, the Hadamard-Zolésio structure theorem is applied, which allows to trace back changes in the objective j solely to variations of the boundary Γ obs (see for instance [30,18]). It is thus possible to define directional shape derivatives via variations of Γ obs in a direction normal to the boundary. Together with the choice of an appropriate shape and tangent space, this allows to represent the sensitivity for j w.r.t. Γ obs as a gradient. This is then interpreted as a deformation to Γ obs and a new discretization mesh for the resulting domain can be computed. By this step the mesh quality of the deformed domain can be ensured as pursued in, e.g., [32,5]. Alternatively, the definition of shape and tangent space includes the surrounding domain Ω, which immediately results in deformation information for the entire mesh (e.g. [22,26,6]) and makes the additional call to a mesh generator superfluous. Typical approaches consider interpreting the shape sensitivity as a force term in linear elastic models described over Ω. The resulting displacement field is then applied as a mesh deformation. Especially in recent works (see for instance [11,3,4]), linear elastic extension equations are considered and, in particular, a very small or even zero first Lamé constant is favored. Moreover, a descent method allows to control the mesh quality from one iteration to the next, i.e. for one deformation. Yet, in the limit of the sequence of design updates, quality is typically lost. This effect is described in, e.g., [23], where variable interfaces must be prevented from overlapping.
In this article we follow a different approach, which gives a higher level of control on the quality of the mesh around the optimal shape. Based on the method of mappings (cf. to [20]), the question for admissible shapes G adm := {F (Ω) : F ∈ F adm } in (1) is translated to the choice of appropriate function spaces, in which a deformation from reference to the optimal configuration is to be found. Here F adm denotes a set of admissible mappings. Starting from a reference configuration Ω, it is then optimized over the transformations F (Ω) yet without explicitly performing mesh deformations. For this purpose the PDE constrained is transformed to the virtual domain as E(y, F (Ω)). The optimization problem then turns into a classical optimal control in the form of This approach is a recent field of studies and applied in, e.g., [17,29,2]. Also based on this approach is the investigation in [11], which is the starting point for the consideration in the present article. Here the problem in (2) is formulated as in terms of a regularization parameter α > 0 and a bound η det > 0 on the determinant of the derivative of the mapping function F . The focus of the investigations therein is on the extension operator S. It is suggested to choose S to be the composition of mappings c → b → w. Here c → b is realized via the solution operator of a Laplace-Beltrami equation on Γ obs . The mapping to the actual displacement, i.e. b → w, is then chosen to be the solution operator of a vector-valued elliptic equation, such as a linear elastic model. It is proven that -under certain circumstances -the domain mapping F is locally a The main focus of our present article is a numerical study of different choices of the extension operator S. It turns out that optimization settings, where larger deformations are to be expected, are a limiting factor for linear operators S. This limitation is due to the fact that the structure of a shape space, that is as large as possible, can hardly be linear since this would require to explain what scalar multiples or sums of shapes are. Yet, with the method of mappings and a linear extension operator S we approximate the set of admissible shapes locally by a linear function space of admissible deformations to a reference configuration.
We thus suggest a nonlinear extension mapping and present numerical studies on the applicability. It should be mentioned that the theory developed so far is not applicable in this case. It only applies to the linear choice of S, which is a special case of the more general consideration in this article.
The motivation for the choice of S in this present work is the observation that, on the one hand, via the condition det(DF ) ≥ η det the local injectivity of F can be ensured. But on the other hand, this limits significantly the subset of admissible transformations F adm and thus affects optimal shapes as outlined in the last section of this article. It is thus the task to find an operator S which prevents det(DF ) ≥ η det from becoming active even for large deformations. We also discuss cases where the reference domain is not of circular shape and demonstrate the performance of the extension and influence on the mesh quality in a deformed domain. The intention of this experiment is to demonstrate that the set of shapes G adm , which is constructed via the mappings from F adm , can be extended significantly and the dependence on the choice of a reference domain Ω can be hidden. In particular, the studies illuminate whether large deformations in the optimization are possible for general reference configurations, which do not fulfill certain properties like convexity or an injective normal vector field.
This article is structured as follows: In section 2 the shape optimization problem is set up and formulated in terms of the method of mappings. Section 3 is devoted to the nonlinear extension model and, furthermore, the derivation of necessary optimality conditions and the presentation of an optimization algorithm. In section 4 numerical studies are conducted and discussed. The article closes in section 5 with a conclusion of the results.

Optimization problem
We carry out our considerations based on a classical optimization problem in the field of fluid dynamics described in [19]. In a d-dimensional, bounded domain Ω with Lipschitz boundary, as sketched in figure 1, we consider the minimization of the following energy dissipation functional where the contour Γ obs of the obstacle Ω obs is assumed to be variable. The spatial dimension is chosen as d ∈ {2, 3}. In (4) the velocity field denoted by v is given in terms of the stationary, incompressible Navier-Stokes equations Together with Γ obs the fluid domain Ω is allowed to change, but the outer boundaries, i.e. Γ in , Γ out and Γ wall , of the experiment are fixed. In (5) p denotes the pressure, v ∞ describes the velocity profile at the inflow boundary, n is the outer normal vector and ν the viscosity. Furthermore, we assume that Γ obs ∩ (Γ in ∪ Γ wall ∪ Γ out ) = ∅ holds during the entire optimization. For the shape optimization of a specimen Ω obs with respect to functionals of type (4), it is essential to exclude trivial solutions. Here, shrinking Ω obs to a point or translations towards Γ wall represents undesired descent directions. Thus, the optimization problem has to be additionally constrained to geometrical conditions. Our benchmark problem is to find optimal shapes of a specimen with a given volume, which remains located in the center of the flow tunnel. This is achieved by fixing barycenter and volume of the obstacle Ω obs with the constraints vol(Ω obs ) = Ω obs 1 dx = const, (6) bc(Ω obs ) = 1 vol(Ω obs ) Ω obs x dx = const.
Since the computation for the barycenter involves the volume of Ω obs itself, these conditions are coupled in principle. Yet, if (6) is fulfilled, the term vol(Ω obs ) −1 in (7) is constant and can thus be factored out. By further assuming that the barycenter of the specimen Ω obs is 0 ∈ R d , it is thus sufficient to require Ω obs x dx = 0. In the following, for a vector-valued function f : R d → R d , we denote by Df the Jacobian matrix with the ordering Df = ∂fi (8) and consider the weak formulation of the PDE constraint (5): Note that within this article we are using the symbol δ · for test functions associated with a given variable. In order to reformulate the optimization problem (4) to (7) as an optimal control problem in appropriate function spaces, we consider from now on the domain Ω as a fixed reference configuration. Let F = id +w with w ∈ W 1,∞ (Ω, R d ) such that F results in an admissible deformation for Ω. For the method of mappings we then consider the state (9), objective (4) and the corresponding state variable v in terms of F (Ω).
By means of standard computations we obtain the weak formulation of the optimization problem pulled back to the reference domain Ω by for all test functions (δ v , δ p ) ∈ V ×Q. The optimization problem (10) to (14) still leaves open the question for the set of admissible mappings F adm . We thus follow the same approach as in [11] and translate it into the form of (3). By reformulating the constraint det(DF ) ≥ η det as a penalty term, we obtain the final optimal control problem min where (·) + denotes the positive-part function. The missing piece is now mapping from a scalar-valued boundary control c to admissible deformation fields w, which is the subject of the next section.

Nonlinear extension operators
Consider the optimal control problem (15). The core of the reformulated shape optimization is the choice of the extension operator S, which links a scalar-valued boundary control c living on Γ obs to a vector-valued displacement field w in Ω. A domain transformation mapping F = id +w is then obtained by the so-called perturbation of identity. In particular, w has to fulfill certain regularity properties as investigated in [11]. It yet turns out in section 4 that for large deformations, i.e. when the reference domain and the optimal configuration differ significantly, linear operators S do not lead to satisfying results. Note that the choice of S significantly influences the set of reachable shapes G adm determined via F adm . It is thus our intention to find S which allows for large deformations without significantly restricting G adm . Simultaneously, the corresponding mesh deformations F (Ω) should exhibit high element qualities for further usage in numerical simulations. The focus of the present article is thus to propose and study nonlinear extensions S given in terms of the solution operator of the coupled PDEs In the equation above ∆ Γ obs denotes the vector-valued Laplace-Beltrami operator. Note that by solving (16) the scalar-valued control c is mapped to a vector valued quantity b. The benefit of this particular extension operator, and especially the nonlinearity η ext (w · ∇)w, which is in the focus of this article, becomes particularly visible for experiments with large deformations as pointed out in section 4.2. A popular choice, as discussed in the introduction, is to define the extension only via the linear term div(∇w + ∇w ). Yet, this restricts the set F adm significantly. This is visible especially for problems in fluid dynamics, as pointed out in section 4, where the reference shape is of spherical type but the optimum to be found is stretched and approximates non-smooth tip and back. Problems arise due to strong compressions of finite elements in the discretization orthogonal to the main deformation direction. This observation motivates to add the nonlinear advection term η ext (w·∇)w, which -geometrically speaking -promotes displacements w where nodes move along large gradients. This results in a homogeneous distribution of finite elements even around approximately non-smooth regions of Γ obs .
For (16), which specifies the mapping from boundary control to domain deformation, the weak formulation is given by: for all δ w ∈ H 1 0 (Ω, R d ) and in terms of η ext ≥ 0. Here D Γ obs denotes the derivative tangential to Γ obs . In (17) the scalar-valued boundary control c ∈ L 2 (Γ obs ) is multiplied with the outer normal vector field n to Ω at Γ obs . Then a vector-valued Laplace-Beltrami equation is solved over Γ obs . This is coupled with nonlinear (18) where the influence of the advection term is controlled via η ext . Note that the linear extension operators investigated in [11] arise as a special case of system (17) and (18). Now that the extension operator is chosen, we can combine the optimization problem (15) with the extension operator (17) and (18) to obtain the Lagrangian where ψ · denotes for each variable the associated multiplier. Note that there is no variable corresponding to the multipliers ψ vol ∈ R and ψ bc ∈ R d . These are the finite dimensional multipliers for the barycenter and volume condition (6) and (7).

Lemma 1. The first order optimality system associated to the Lagrangian
Proof. The derivatives of L are obtained by utilizing standard rules of differentiation. Note that we particularly use the following identities For the derivative of the penalty term we utilize that Recall that the condition det(DF ) ≥ η det in the problem formulated in (3) is realized via the penalty term β 2 Ω ((η det − det(DF )) + ) 2 dx. The corresponding term in (20) of the optimality system in lemma 1 is non-differentiable due to the positive-part function (·) + . Following the discussions in [11, sec. 3.5] the mapping w → −β Ω (η det −det(DF )) + Tr((DF ) −1 Dδ w ) det(DF ) dx is semismooth and one can compute an element from the generalized derivative in direction δ as In the following we briefly present a solution algorithm for the optimality system (20) to (30). For this purpose we pursue a similar approach as in [11]. The core of this method is to solve the nonlinear shape optimization problem (15) for a decreasing sequence of regularization parameters α k , starting from α 0 = α init until the desired level α target is reached. Because each subsequent optimization problem k is nonlinear, this approach benefits from utilizing the known values y k := (w, v, p, b, c, ψ w , ψ v , ψ p , ψ b , ψ vol , ψ bc ) k as initial guess in the (k + 1)-th iteration. Algorithm 1 summarizes this procedure. Since parts of the optimality system are non-differentiable, we apply a semismooth Newton's method. In section 4 we

Numerical results
This section is devoted to different numerical case studies of stationary, incompressible Navier-Stokes shape optimization problems. The purpose is to illuminate features of the nonlinear extension operator S proposed in section 3. In particular, the benefit for optimization benchmark problems, which involve large deformations from the reference configuration to optimal shapes, is numerically investigated. It is moreover discussed how the local injectivity can be extended to globally injective transformation mappings by adding an artificial volume to the aerodynamic specimen. Furthermore, algorithmic solvability of the optimality system (20) to (30) is addressed in the end of this section.
The experimental settings for the tests are chosen to be comparable in 2d and 3d, respectively. The holdall domain G ∈ {G 2d , G 3d }, which reflects the flow tunnel in the experiment, is chosen as Let δ denote the diameter of the flow tunnel G. We then fix the velocity at the inflow boundary Γ in by In all experiments where the specimen Ω obs is a circle or sphere, the radius is given by r = 0.5 and bc(Ω obs ) = 0 ∈ R d . The discretization of all appearing PDEs is carried out with standard, piecewise linear P1 finite elements. In order to guarantee stability, we follow the pressure stabilized Petrov Galerkin approach (see for instance [14]), which utilizes an additional term for the pressure p and its adjoint variable ψ p . The system under consideration is thus enriched by the two equations All computations related to the finite element method are carried out using the GETFEM++ library [24]. We utilize a parallel version of the library, which relies on PARMETIS [16] for mesh partitioning and load balancing. All linear systems are handled via the parallel factorization solver MUMPS [1]. Both, 2d and 3d discretization meshes are produced with the GMSH toolbox [8] and the Delaunay algorithms therein. If not stated otherwise, all 2d experiments follow the strategy of algorithm 1 with the choice α init = 1e−4, α dec = 1e−1 and α target = 1e−10.   Our first numerical study demonstrates the nonlinear extension equation for large deformations of a non-convex shape. The reference domain Ω is chosen such that the specimen is described by a B-spline curve Γ obs given in terms of 6 control nodes. The situation is depicted in figures 2 and 3.

Non-convex shapes with large deformations
The relevance of this test case is to investigate the performance of the proposed approach for reference domains where the normal vector field n does not homogeneously point in all directions as for a circular shape. The influence of the normal vector is significant since it initially links the scalar-valued control c to a vector-valued quantity as can be seen in (16). Section 4.3 is devoted to an experiment where several directions are underrepresented in the discretization of the normal vector n due to the shape of Ω obs .
In figure 2 the magnitude of velocity v 2 , computed in the undeformed state Ω, i.e. when c = 0, is depicted. The right-hand side shows the velocity according to deformation F = id +w in terms of the optimal control c after solving the optimality system given by (20) to (30). Furthermore, the optimal mapping F can be seen in figure 3 in the displacement and deformation of discretization elements. The relocation of triangles shows the effect of the nonlinear advection in the extension operator. A deeper look on this effect and the resulting mesh quality is provided in section 4.5.
In this experiment the viscosity of the fluid is chosen as ν = 0.01. The holdall domain G is as described above. It consists of 382 surface segments on Γ obs and 10 106 triangles in Ω. Barycenter and volume of Ω obs in the reference configuration are given by bc(Ω obs ) = (0.030 778 4, −0.035 759) and vol(Ω obs ) = 0.809 041. Thus, an optimal shape also undergoes a small translation since bc(F (Ω obs )) = 0 is required.
The essential settings in terms of shape optimization are the parameter η det and η ext . Here we choose η ext = 3.0, which leads to the condition det(DF ) ≥ η det with η det = 5e−2 being inactive. We can explain the homogeneously and smoothly deformed mesh due to this fact. In contrast, section 4.2 shows examples where det(DF ) ≥ η det is active close to the tip of the optimal shape and how the displacement w and thereby the mesh quality in F (Ω) is affected. In this section we visualize the influence of the choice of η det and η ext on the optimization. This demonstrates how the set of admissible shapes F adm is determined thereby. The underlying experiment is a flow in G over a circular specimen as described at the beginning of section 4. The viscosity is again chosen to be ν = 0.01. The domain is discretized with 244 segments on Γ obs and 12 640 triangles in Ω.

Influence of factors η det and η ext
First, we observe the influence of η det on the set of admissible shapes F adm . Figure 4 visualizes how the condition det(DF ) ≥ η det acts on the optimal shape F (Γ obs ). Here we choose η det ∈ {0.5, 0.25, 0.2, 0.1} beginning with the largest and then decreasing values. This condition can be interpreted such that the allowed, local change of volume in Ω is relaxed from top to bottom in figure 4. In this experiment it turns out that in the last computation with η det = 0.1 the condition is inactive. Here η ext = 3.0 is fixed in all computations. The next experiment follows the same setup with the only difference that η det = 5e−2 is now fixed and η ext ∈ {0.0, 0.25, 0.5, 1.0, 2.0, 3.0} takes increasing values. The mesh deformations, resulting from the optimal solution F = id +w, are visualized in figure 5. Each of the subfigures shows a clip of size 0.1 × 0.08 around the tip of the optimal shape. Besides the significantly improved mesh qualities and more adequate set F adm we also observe that the Newton solver benefits from the appropriate choice of η det and η ext . As soon as the condition det(DF ) ≥ η det becomes active, the optimality system (20) to (30) is not differentiable any further and the solver switches to semismooth Newton's method. This effect is already documented in [11] for the case of Stokes flows and linear extension operator S.

Extending the local-only injectivity
This section focuses on an effect that is likely to appear for non-spherical reference domains. In particular, for the aerodynamic experiments considered in this article it might happen that the upper surface of the obstacle overlaps the lower one. Especially for large deformations from reference to optimal shape and for shapes that are streched parallel to the flow axis, we encounter effects as depicted in figure 6a for the experiment shown in figure 7. In other words, Ω → F (Ω) is not globally injective in this situation. This   is due to the fact that the condition det(DF ) ensures injectivity of F only locally but not globally. In the following we propose a modification of the extension operator S in order to extend the injectivity. Recall that in the setting followed up to here the obstacle domain Ω obs is treated as void and there is no discritization within. We now consider the operator S on the entire holdall domain G = Ω ∪ Ω obs in contrast to the state equation that remains in Ω. Moreover, the condition det(DF ) ≥ η det > 0 is now required on G. Thus, the displacement field is defined by w ∈ H 1 0 (G, R d ). Moreover, we reformulate the weak formulation of the extension operator S given in (18) to for all δ w ∈ H 1 0 (G, R d ) and appropriate η ext ≥ 0. Simultaneously, the penalty term, which enforces the local injectivity, in (19) changes to For this experiment we choose ν = 0.01 as in the previous sections. The holdall G has the same outer dimension and the specimen Γ obs is an ellipse with semimajor-axis r 1 = 2.7, semiminor-axis r 2 = 0.2 and barycenter bc(Ω obs ) = (0, 0) . Its surface is subdivided in 884 segments. Further, G is discretized by 36 360 triangles, 29 930 in Ω and 6430 in Ω obs . Figure 6 visualizes the effect of the mapping F on the discretization grid. In figure 6a the optimal solution for α = 1e−2 is shown. Note that in this particular case we stop the optimization for a larger value, since this already leads to singularities. Figure 6b depicts det(DF ), which is again bound away from zero by η det = 5e−2. It can be seen that, although this condition is inactive, the non-injective mapping can not be prevented.
The same experiment is then conducted with the changes proposed in the beginning of this section, which leads to the values of det(DF ) shown in figure 6c. Now Ω obs is discretized and S also acts on the interior of the specimen. Here the optimization is performed with the setting α init = 1e−4, α dec = 5e−1 and α target = 1e−10. The resulting optimal solution is visualized in figure 7 where figure 7a shows the reference domain and the velocity field computed for this configuration. Figure 7b depicts the domain F (G) and the velocity field computed on F (Ω). Note that the relatively fine grid is chosen at the front and the back of the shape due to the large curvature of Γ obs in these regions. This experiment turns out to be more challenging than, e.g., a spherical reference shape since on coarse grids the normal vector field in these areas tends to be underresolved. From a computational point of view, it is attractive to have a coarse grid in Ω obs , as chosen in the center of the specimen, to reduce cost for the solution of the operator S.

Three-dimensional results
In this section we perform a three-dimesional optimization experiment as a proof of concept. As already observed in [11] for the Stokes experiment, more care has to be taken for the decrease-strategy of α in algorithm 1. Especially the semismooth Newton solver shows to be challenging w.r.t. to convergence when the condition det(DF ) ≥ η det becomes active.
The experiment shown in figure 8 is within the framework described at the beginning of section 4. The flow tunnel Ω is discretized by 632 093 tetrahedrons and the surface of the spherical specimen in the reference configuration Γ obs consist of 8558 triangles. Further, the viscosity is chosen to be ν = 0.01 and in algorithm 1 we set α init = 1e−4, α dec = 5e−1 and α target = 1e−6. The results shown here are obtained with an extension factor of η ext = 15. We visualize the impact of the optimization on the fluid by stream lines of the velocity field in figure 8. This figure also shows the effect of the particular operator S on the quality of the surface mesh when it undergoes the optimal deformation F . The combination of Laplace-Beltrami (17) and the nonlinear extension equation (17) leads to a homogeneous distribution of triangles on the surface Γ obs . It can be observed that this is due to tangential components in w| Γ obs . This is a benefit of a vector-valued extension equation over approaches which utilize a static extension of the normal vector field in order to extend the boundary control to the surrounding volume. Furthermore, figure 9 shows a zoom-in to the tip of the deformed domain F (Ω). Here we can see a crinkled clip in the x 1 x 2 -plane with x 3 = 0, which shows the quality of the tetrahedrons.

Quantification of the influence of η ext on mesh quality
This section presents numerical experiments which investigate the influence of the nonlinear extension operator S on the mesh quality in 2d and 3d (cf. figure 10). Recall that the discretization mesh is not actually deformed within the optimization. We though deform the reference domain Ω according to the optimal control w = S(c) and the corresponding deformation F = id +w. The 2d experiment is conducted on the same computational domain as before with a circular specimen, 312 surface segments and 6168 triangles in Ω. The fluid viscosity is chosen to be ν = 0.01. Figure 10a visualizes the influence of η ext ∈ [0, 1] on the mesh quality of F (Ω). It is measured by the ratio of radii of largest inscribed and smallest circumscribed circle, where the plot shows the value of the worst triangle.
This experiment quantifies the effect which is already visualized in figure 5. For a shape optimization with large deformations from reference to optimal configuration, i.e. w L 2 (Γ obs ) is relatively large, a pure linear extension operator S does not reliably lead to satisfying mesh qualities. Moreover, it can be seen that in this particular experiment there is a saturation effect of the nonlinearity in S starting at approximately η ext ≈ 1.5. Figure 10b shows the results of a similar experiment in 3d. Here a mesh is chosen with 6040 surface triangles on Ω obs and 147 385 tetrahedrons in Ω. Note that we decrease the viscosity to ν = 0.1 in this experiment in order to be able to obtain results for η ext < 3.0. In 3d quality   is measured by the radius ratio of smallest circumscribed sphere to the largest inscribed one. Again the worst element is visualized. Also note that the y-axis is in log-scale. In this setting it turns out, that the effect of compressed cells near the tip and back of the shape, which is stretching due to a decrease in α, is stronger than in 2d. We explain the solver failure due to the semismoothness in the optimality system, which becomes active in a significant number of finite elements in this situation. However, it can be observed that, starting with approximately η ext ≈ 8, a saturation is possible, where the mesh quality of F (Ω) remains adequate for further numerical computations.

An iterative optimization algorithm
In the previous sections we solve the nonlinear, non-smooth optimality system with the direct solution strategy given in algorithm 1. Moreover, a direct solver library is applied to the resulting linear systems within semismooth Newton's method. This approach is clearly limited due to the high memory requirement. Especially, when the state equation results from a time-dependent problem, this procedure becomes impracticable. Hence, in this section we focus on a numerical study of decoupling system (20) to (30). This approach is summarized in algorithm 2. We demonstrate that it is possible to decouple the solution process of state (23) and (25), adjoint

9:
Solve (20), (21) and (26) to (30) for (w, b, c, ψ w , ψ b , ψ vol , ψ bc ) +1 with semismooth Newton's method and regularization parameter α k 10: (22) and (24) and shape related equations, i.e. (20), (21) and (26) to (30), from each other. On the one hand, this allows to reuse existing solvers for the state equation and embed them into the shape optimization framework. On the other, the memory requirement for linear solvers significantly reduces. Moreover, the semismooth part (20) is split from the other equations and a solver can be particularly tailored for this purpose. Algorithm 2 operates on the nonlinear optimality system as a fixpoint strategy. In an outer loop it is again iterated over a decreasing regularization parameter α as in algorithm 1. Thus, approximate solutions for the optimization problem according to α k are utilized as initial guess for the nonlinear solver in iteration k + 1. Yet, unlike in the direct approach, the subproblems are only solved approximately by a fixpoint iteration, which solves the decoupled equations of the optimalitiy system in turns. The termination criterion for this inner loop is the relative change in the control variable c measured in the L 2 (Γ obs )-norm.
In figure 11 the results of one run of algorithm 2 are shown. The underlying optimization experiment is a 2d computation on the same grid as in section 4.5 with 312 surface segments and 6168 triangles in Ω, α dec = 0.5, α init = 1.0, α target = 2e−7, ν = 0.1 and η ext = 1.5. Note that the initial value of α is significantly larger then the choices made for algorithm 1. Figure 11 shows the required inner iterations until the condition c +1 − c L 2 (Γ obs ) c +1 L 2 (Γ obs ) < is fulfilled for = 1e−2. Futhermore, the value of the objective J (cf. (15)) is visualized. It is computed in algorithm 2 in line 10 at the end of one inner loop. Notice the jumps in the objective function between iteration 5 and 20. In our experiments it turns out that this is an effect that both influences the minimal possible α dec and α init . In this setting a total number of 53 inner iterations, i.e. solutions of the state equation, are required to reach the optimal shape. This numerical study can thus be seen as a proof of concept how to reduce the computational costs of the large, coupled, nonlinear system (20) to (30). Thus, the proposed method is applicable to more complex problems, such as non-stationary Navier-Stokes flows.

Conclusion
In this article we have proposed and numerically demonstrated choices of nonlinear extension operators within the method of mappings for aerodynamic shape optimization. These operators are based on the idea that an additional, nonlinear advection term leads to a rearrangement of discretization cells along the major direction of deformations.
The main goal we have achieved is to circumvent mesh degeneracy effects that appear under large deformations when the extension of the boundary control is chosen according to linear elastic models. Especially in the underlying aerodynamic drag minimization, where optimal shapes tend to become stretched in flow direction and compressed in the orthogonal directions, we have numerically investigated how mesh quality can be preserved.
We have also demonstrated one possibility to decouple the solution process of the optimality system in order to overcome issues of computational complexity. Moreover, we have studied how the set of admissible shapes depends on the nonlinearity of the operator and how the local injectivity of mappings can be extended to large deformations. Since the proposed methodology is formulated in function spaces without taking a specific discretization into account, another benefit of this approach is that it naturally allows to introduce concepts like adaptivity. An important field for future investigations is a detailed description of properties of the set F adm , which is constructed in terms of the nonlinear extension operator S.