{'fraction_non_alphanumeric': 0.04650752963793656, 'fraction_numerical': 0.012512015379685998, 'mean_word_length': 4.386228319958581, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pre-Shape Calculus: Foundations and Application to Mesh Quality Optimization

Deformations of the computational mesh arising from optimization routines usually lead to decrease of mesh quality or even destruction of the mesh. We propose a theoretical framework using pre-shapes to generalize classical shape optimization and calculus. We define pre-shape derivatives and derive according structure and calculus theorems. In particular, tangential directions are featured in pre-shape derivatives, in contrast to classical shape derivatives featuring only normal directions. Techniques from classical shape optimization and -calculus are shown to carry over to this framework. An optimization problem class for mesh quality is introduced, which is solvable by use of pre-shape derivatives. This class allows for simultaneous optimization of classical shape objectives and mesh quality without deteriorating the classical shape optimization solution. The new techniques are implemented and numerically tested for 2D and 3D.


Introduction
Solutions of PDE constrained optimization problems, in particular problems where the desired control variable is a geometric shape, are only meaningful, if the state variables of the constraint can be calculated with sufficient reliability. A key component of reliable solutions is given by quality of the computational mesh. It is well-known that poor quality of elements affect the stability, convergence, and accuracy of finite element and other solvers.
We propose a unified framework using so-called pre-shapes. In this setting, both shape optimization and mesh quality optimization problems can be formulated at the same time. We give a class of problems called pre-shape parameterization tracking problems, which can act as regularizations for shape optimization problems. These problems can be solved for volume and surface meshes with arbitrary dimension, and yield numerical algorithms similar to so called mesh deformation methods, which reallocate nodes of numerical meshes according to targeted element volumes to improve mesh quality. At the same time, the proposed framework is suitable to derive a calculus mimicking classical shape calculus. This enables formulation of efficient routines solving shape optimization problems, which at the same optimize quality of the surface mesh representing the shape, without noticeable additional computational cost or interference with the original shape optimization problem.
With this, desired surface and surrounding volume mesh quality are ensured during shape optimization.
In this paper we will establish the theoretical foundations of pre-shape optimization and calculus, and show its connection to classical shape calculus. We implement these methods in form of a pre-shape gradient descent to achieve a targeted quality of volume and surface meshes without a shape optimization target. In [34], we build on achievements of this paper, and introduce theoretical and numerical results to solve shape optimization problems while simultaneously improving volume and shape mesh quality according to targeted node distributions. The techniques elaborated in [34] leave optimal shapes invariant and offer minimal additional numerical costs. They also permit use of different metrics to represent gradients. We compare pre-shape mesh regularizations for various metrics in [34] for a hard to solve shape optimization problem.
Literature Review We give a brief overview of techniques related to the ones treated in this article. Our methods are not related to mesh untangling andrelaxation, edge swapping or remeshing strategies such as [18,19,1,25]. Of course, as there is a vast amount of literature concerning mesh generation and improvement strategies thereof, we can only give a selective overview. Two very prominent classes of techniques for mesh quality improvement are the so called mesh deformation method and methods based on Laplacian smoothing. Mesh deformation methods go back to a theoretical result initially proposed by Moser in [36], extended by Banyaga [3] and by Dacorogna and Moser [11]. This gave rise to mesh deformation methods pioneered in [31] by Liao and Anderson, which redistribute mesh vertices such that uniform cell volumes are achieved. Mesh deformation methods are powerful, because they prevent mesh tangling while offering precise control over the element volumes. The original method was further developed in various directions by Liao et. al [7,32,8,51]. These advances allow to target non-uniform cell volume distributions, make deformation methods applicable to time-dependent problems, and shows its use in higher order mesh generation methods. Also, the combination of multigrid-and mesh deformation methods was analyzed and implemented by Turek et. al. in [48,21,22]. A different family of mesh quality improvement techniques are those based on Laplacian smoothing [17,18,44,50]. They do not necessarily track for cell volume distributions, but instead improve quality by averaging or smoothing vertex coordinates more or less specifically. Several strategies for increasing mesh quality not based on mesh deformation methods mentioned or Laplacian smoothing exist. In the context of shape optimization and -morphing, these include correcting for errors in Hadamard's theorem due to discretization [16], adding non-linear advection terms in shape gradient representations [37], approximating shape morphing by volume-preserving mean-curvature flows [28], and use of techniques related to centroidal Voronoi reparameterization in combination with eikonal equation based non-linear advection terms for representation of shape gradients [41].

General Theory for Pre-Shape Calculus
Pre-Shape Spaces In order to provide theoretical grounds for the numerical procedures we are about to elaborate in the subsequent paper, we need to specify a framework for the objects, 'shapes', for which we seek to optimize. Several possible theories and techniques exist in order to precisely formulate the notion of shapes. For example, shapes can be viewed as sets together with corresponding characteristic functions in an ambient space, leading to an approach which emphasizes geometric measure theory as in [12] by Delfour and Zolésio. We choose a differ-ent setting, namely a shape space approach using infinite dimensional differential geometry, since it naturally permits to view the shape and its parameterization at the same time. This is key to extend numerical routines in a way that optimizes parameterizations, i.e. meshes, in a desired way without interfering with the shape optimization taking place. For an excellent overview of shape spaces we refer to [4], from which we borrow several definitions for the following introduction to shape spaces considered in this article.
For the rest of this article, let M be an n-dimensional, orientable, path-connected and compact C k,α -or C ∞ -manifold. Further, we will use R n+1 as the ambient space for building our theory. In particular only shapes of codimension 1 are considered.
Denote by Diff(M ) the regular Lie-group of C ∞ -diffeomorphisms of M onto itself. Then the space B e of unparameterized C ∞ -shapes in R n+1 is defined by (cf. [39]) where Emb(M, R n+1 ) is the space of all C ∞ -embeddings of M into R n+1 and Diff(M ) is acting on the right. For avoidance of confusion, we remember that in this context a C ∞ -embedding is an injective, smooth map ϕ : M → R n+1 , which has injective first derivative everywhere, i.e. ϕ is an injective immersion. The resulting space B e (M, R n+1 ), also called non-linear Grassmannian (cf. [27, 44.21ff]) or differentiable Chow-variety. It forms a smooth Hausdorff manifold, whose elements can be regarded as unparameterized hypersurfaces of R n+1 (cf. [35,Corollary 3.3]). In the following we will abbreviate B e (M, R n+1 ) by B n e , still having the implicit relation to the manifold M and its dimension in mind. This space can be equipped with various metrics, which means shape optimization can be regarded as optimization on an infinite dimensional Riemannian manifold (cf. [42]). Notice that B n e is not a manifold if C ∞ -regularity is replaced by Hölder-or Sobolev-regularity. In this more general setting, resulting spaces have a diffeological structure (cf. [49]). To acquire intuition, a graphical visualization of B n e is given in fig. 1. For our purposes it is not enough to view shape optimization as optimization in B n e . Instead, we will exploit additional structures on the space Emb(M, R n+1 ) induced by the action of Diff(M ) and base our framework as optimization in Emb(M, R n+1 ). Elements ϕ ∈ Emb(M, R n+1 ) can be interpreted as parameterized shapes in the ambient space R n+1 , whereas elements of Diff(M ) acting on the right can be seen as reparameterizations. The authors of [4, Ch. 1.1] call Emb(M, R n+1 ) pre-shape space, an expression we will borrow for the techniques we will build in this paper. Notice that the term pre-shape space is used differently depending on the literature, e.g. in [26], where the authors use this term for the space of labeled landmarks which are equivalent under translation and scaling.
The additional structure of parameterized shapes ϕ ∈ Emb(M, R n+1 ) compared to unparameterized shapes Γ ∈ B n e enables to not just view the shape itself, but also to distinguish various types of discretizations in the ambient space and corresponding numerical meshes. Even further, this concept enables to control the parameterization of the hold-all domain itself, allowing for control of the way volume meshes are discretized. The structure for this is given by the fact (cf. [27,Thrm. 44.1], [6]) that the quotient map makes Emb(M, R n+1 ) the total space of the smooth principal fibration with Diff(M ) acting as the structure group or standard fiber, and B n e being the base space, which goes back to [6]. As a reminder, a fiber bundle is a manifold, which locally looks like a product space B × F , where B corresponds to the base space, and F corresponds to the standard fiber. In our context, this means the pre-shape space Emb(M, R n+1 ) is the collection of parameterized shapes, which locally looks like 'Shape'×'Parameterization'. However, this relationship holds only locally, and the global structure of the pre-shape space Emb(M, R n+1 ) is much more complex. The situation is graphically sketched in fig. 2. An application of the bundle projection π to a parameterized shape ϕ ∈ Emb(M, R n+1 ) results in its unparameterized shape π(ϕ) ∈ B n e in the base space. Hence, we can view the fiber π(ϕ) as the collection of all parameterizations of the shape ϕ(M ). It is important to avoid confusion of ϕ(M ) and π(ϕ), which are both called shapes. The first interprets shapes as subsets ϕ(M ) ⊂ R n+1 , the latter as equivalence-classes, i.e.
The equivalence class interpretation is the collection of parameterizations corresponding to a certain shape in R n+1 . In order to formulate an analogue of shape calculus in Emb(M, R n+1 ), we need to characterize the tangential bundles T Emb(M, R n+1 ) and T B n e , as well as their relations. For this, we make use of results by Michor and Kriegl [27].
Since we assume M to be compact, the respective tangent bundle of the preshape space is isomorphically given by The fiber-bundle structure leads to a decomposition of the tangent bundle of the total space T Emb(M, R n+1 ) into the so called vertical bundle, defined as ker T π ⊂ T Emb(M, R n+1 ), and the horizontal bundle. Since we only deal with compact and orientable n-dimensional manifolds M , the existence of outer normal vector-fields n on ϕ(M ) ⊂ R n+1 is guaranteed. In the following, let ., . 2 denote the L 2 -scalar product. Thus we obtain π(ϕ 0 ) Figure 2: Illustration of the pre-shape space Emb(M, R n+1 ) for M = S 1 . To illustrate the parameterization interpretation of fibers π(ϕ), the same four points are mapped from M to ϕ i (M ). where is the space of tangential vector fields on ϕ(M ) in the ambient space R n+1 and being the space of normal vector fields on ϕ(M ). The tangential fields are parts of the vertical bundle, whereas the normal fields constitute the horizontal bundle part. This also gives the well-known characterization of the tangential bundle of the classical shape space B n e via normal vector fields, i.e.
As previously, we also visualize the situation for tangential bundles in pictures fig. 3 and fig. 4.
Pre-Shape Calculus Next, we introduce a suitable notion of objective functionals. We are inspired by [12,Ch. 4.3.1], where shape functionals J are functions on a set of admissible shapes A, which are considered to be a subset of the power set A ⊆ P(R n+1 ). This is a set-theoretic approach, because the power set P(R n+1 ) is Figure 3: Illustration of a shape tangential vector from T B n e for M = S 1 .
the set of all subsets of R n+1 . Since we can canonically associate every equivalence class π(ϕ) ∈ B n e with its set ϕ(M ) ⊂ R n+1 , we get the following definition for the special set of admissible shapes B n e .
Definition 1 (Shape and Pre-Shape Functionals). Let M be an n-dimensional, orientable, path-connected and compact C ∞ -manifold. Consider the shape space B n e as defined in eq. (1) and the space of embeddings Emb(M, R n+1 ). Then a function J : B n e → R is called shape functional, and a function is called pre-shape functional.
The nomenclature pre-shape functional for functions as in eq. (10) is motivated by regarding Emb(M, R n+1 ) as a pre-shape space as done by Michor et. al. in [4,Ch. 1.1]. Since optimization is classically taking place in shape spaces as opposed to pre-shape spaces, we will highlight some of their correspondences and relations. The following definition is motivated by the construction of the shape space B n e in eq. (1).
If J has shape functionality for all ϕ ∈ Emb(M, R n+1 ), we say J has shape functionality.
In order to give optimality criteria for the pre-shape optimization problems and to formulate derivative based optimization algorithms, we need to introduce a shape derivative analogue for Emb(M, R n+1 ). Also, it is desirable that the analogue is compatible with the classical Eulerian derivative, as for example found in [12,Ch. 4.3.2] or [43,Ch. 2.1]. This motivates us to proceed in the fashion of classical shape Figure 4: Illustration of vectors from T Emb(M, R n+1 ) with pure tangential/vertical components for M = S 1 . Note that four points are added to illustrate the parameterization interpretation of fibers π(ϕ).
optimization by defining a pre-shape derivatives based on families of deformations perturbing the image space. We show their relation to classical shape derivatives, and then give a structure theorem similar to the Hadamard-Zolésio structure theorem (cf. [12, Ch. 9, Thrm. 3.6]). Shape calculus or sensitivity analysis of classical shape optimization (cf. [23,Ch. 3], [5]) will carry over to pre-shape spaces naturally.

Remark 1 (Validity of Pre-Shape Theory for different Regularities of Shapes).
We want to remind the reader, that the choice of C ∞ -regularity for pre-shapes in Emb(M, R n+1 ) is not necessary to introduce the concepts of this section, but merely serves as an exemplary situation. It is clear, that the same definitions apply for embeddings ϕ of Sobolev-or Hölder-regularity. In these cases test functions and directions V of course need to have according regularity.
Definition 3 (Perturbation of Identity and Pre-shape Derivatives). Let J be a preshape functional (not necessarily having shape functionality), ϕ ∈ Emb(M, R n+1 ) and V ∈ C ∞ (R n+1 , R n+1 ). Then the family of functions is called perturbation of identity of ϕ in direction V for t ∈ [0, τ ) and some τ > 0.

The limit
is called pre-shape derivative for J at ϕ ∈ Emb(M, R n+1 ) in direction V , if it exists and is linear and bounded in V .
The perturbation of identity for shapes at Notice that this is a set, in contrast to the perturbation of identity for pre-shapes eq. (12) which is a function in Emb(M, R n+1 ). Shape derivatives of a shape functional J are given by The difference quotients defining pre-shape and shape derivatives use completely different objects, therefore their difference is significant. Their relationship is explored in proposition 1 and theorem 1. The next proposition shows a result relating shape differentiability of classical shape optimization and pre-shape derivatives.
Proposition 1 (Shape Differentiability implies Pre-Shape Differentiability). Consider a shape functional J : B n e → R. Then it has a canonical extension to a pre-shape functional where π is the bundle projection as in eq. (2). Further, there is a one-to-one correspondence of shape functionals J and pre-shape functionals J with the property of shape functionality. Additionally, if J is shape differentiable in the classical sense, then its extension J is pre-shape differentiable.
Proof. One-to-one correspondence of pre-shape functionals with the property of shape functionality as in definition 2 and classical shape functionals as in eq. (9) is clear. On the one hand, every canonical extension eq. (16) of a classical shape functional J has shape functionality. On the other hand, every pre-shape functional J having shape functionality gives rise to a shape functional J fulfilling eq. (16), as every fiber π(ϕ) is the orbit of a ϕ by Diff(M ) acting from the right. The pre-shape differentiability assertion in proposition 1 holds, since J • π extends the values of J constantly onto the fibers of Emb(M, R n+1 ). For a fixed ϕ ∈ Emb(M, R n+1 ), a case analysis for directions V ∈ C ∞ (R n+1 , R n+1 ) being either tangential or normal at ϕ(M ) can be made. In case of horizontal, i.e. normal, directions V we recover the classical shape derivative DJ . On the other hand, the pre-shape derivative in vertical directions can be represented as a differential via curves on ϕ(M ), which combined with J • π extending J constantly on fibers gives a vanishing pre-shape derivative. Linearity and boundedness of D(J • π)(ϕ)[V ] in V are easy to see due to its vanishing for tangential components of V together with linearity and boundedness of the classical shape derivative DJ by assumption.
We can now situate classical shape optimization problems in the context of optimization in pre-shape spaces Emb(M, R n+1 ) for suitable manifolds M . But first, we observe that a unique solution ϕ(M ) of a shape optimization problem has multiple parameterizations in general. For shape optimization problems posed in the pre-shape space Emb(M, R n+1 ), this leads to non-uniqueness of solutions, which might at first seem like a disadvantage. However, due to non-uniqueness up to elements in the solution fiber π(ϕ), it is possible to demand additional properties for the pre-shape solution. This gives several opportunities to enhance numerical shape optimization routines, while at the same time narrowing down the amount of non-uniqueness of pre-shape solutions to a reasonable level. We will exploit this in upcoming works, such as [34]. For example, increasing mesh quality while not changing the shape at hand can be viewed as a condition posed on a shape optimization problem selecting a pre-shape in a given fiber. Proposition 1 also offers a possibility to transfer results concerning shape differentiability of classical shape functionals to the pre-shape setting without the need for new proofs. In particular, existence of stationary points in B n e is carried over to Emb(M, R n+1 ) as existence of stationary fibers. Hence proposition 1 shows that pre-shape optimization is in some sense a canonical generalization of classical shape optimization.
The definition of material-and shape derivatives found in classical shape optimization and structural sensitivity analysis literature (cf. [5, Definition 1, Definition 2], [23, Ch. 3.3.1]) possesses useful properties for practical applications. In particular, through the use of material derivatives, it is often straightforward to derive expressions for shape derivatives of integral quantities. We proceed by extending the notion of material derivatives from the classical context to the pre-shape calculus framework to harness these practical benefits.

Definition 4 (Pre-Shape Material Derivative). Consider a family of functions {f
if the limit exists. Here, ϕ t is the perturbation of identity for pre-shapes (cf. eq. (12)) and The careful reader might notice the similarity of classical shape and pre-shape material derivatives. However, the main difference is a possible dependence of functions f on parameterizations of shapes/domains they are defined for. Still, both notions coincide if the pre-shape functional has shape functionality, as we will see in corollary 1 coming from the main structure theorem 1.
The definition of the material derivative can be generalized to functions and domains of weaker regularity, such as Sobolev functions and open subset Ω ⊂ R n+1 with Lipschitz boundaries. A necessity for this comes from the fact, that state solutions stemming from PDE constrained shape optimization problems need a well-defined material derivative for sensitivity analysis to be convenient. This is done in the same way as with the classical shape material derivative (cf. [23, p. 111]).
It is important to notice that the family f ϕ can be seen as a function f : In the first component, the perturbation of identity for pre-shapes comes into play, which differs from the classical shape-material derivative. This leads to the following decomposition of the material derivative, which is similar to classical shape calculus, e.g. given by Haslinger and Mäkinen in [23, p. 111 In the following we give a characterization of the pre-shape derivative in the style of the Hadamard-Zolésio structure theorem as found in [

Theorem 1 (Structure Theorem for Pre-Shape Derivatives).
Let J : Emb(M, R n+1 ) → R be a pre-shape differentiable pre-shape functional (not necessarily having shape functionality) and let ϕ ∈ Emb(M, R n+1 ). Denote by n ϕ(M ) the outer normal vector field of a shape ϕ(M ) for a ϕ ∈ Emb(M, R n+1 ). Then the following holds: (i) The support of DJ(ϕ) is given by (ii) There exist continuous linear functionals g T : and where is the trace operator and n the outer unit normal vector field on ϕ(M ).
(iii) If J has shape functionality, then for all ϕ ∈ Emb(M, R n+1 ), we have g T = 0 and where J : B n e → R is the natural shape functional corresponding to J by J • π = J. In particular g N corresponds to the distribution of the classical Hadamard-Zolésio structure theorem.
Proof. Since we are in a different situation than the classical Hadamard-Zolésio structure theorem [12, Ch. 9, Thm. 3.6] for (i) and (ii), we give proofs for these ourselves. For Consider the perturbation of identity ϕ t for V of ϕ as in eq. (12). (13), which immediately gives us (i).
For (ii), the proof, in some extent, follows analogous reasoning as in [12, Ch. [40,Def. 6.22] for definition of supports of distributions). In addition, we can use that it is contained in the Schwartz space C ∞ c (R n+1 , R n+1 ) due to compactness of ϕ(M ). This gives us the tempered distribution property. Then g T and g N are defined by restriction to the vertical and horizontal part of V , recurring on decomposition of the tangent bundle in horizontal and vertical components eq. (5), giving us (ii).
For (iii), let J have shape functionality and let V ∈ C ∞ (R n+1 , R n+1 ). On ϕ(M ), we can decompose V into normal and tangential components. For the tangential part, we can follow analogous arguments as in the proof of proposition 1, giving us a curve ϕ t in the fiber π(ϕ) generating the pre-shape derivative as a differential at a given ϕ. As ϕ t is running on the fiber of ϕ and J has shape functionality (cf. definition 2) in ϕ by assumption, J(ϕ t ) is constant for all t, rendering g T = 0 by eq. (20). Further, by eq. (20) the pre-shape derivative DJ(ϕ)[V ] reduces to g N acting on normal directions. With proposition 1, shape functionality of J leads to a well defined shape functional J : B n e → R with J • π = J. As the tangential part of V has no impact on DJ(ϕ)[V ], we can find a horizontal curve ϕ t generating DJ(ϕ) [V ]. The representative ϕ t either creates a trivial curve π(ϕ t ) in B n e , which leads to eq. (23) being 0 on both sides, or a non-trivial curve π(ϕ t ) in B n e . If π(ϕ t ) is non trivial, we have for the shape derivative DJ (π(ϕ)) and pre-shape derivative DJ(ϕ), resulting in eq. (23). By association of ϕ(M ) with π(ϕ) ∈ B n e this also shows that g N corresponds to the distribution in the classical Hadamard-Zolésio structure theorem (cf. Structure theorem 1 gives an intuitive way to understand the pre-shape derivative eq. (13) and the relation between shape functionals and pre-shape functionals. Part (i) of theorem 1 has the same meaning as in the classical Hadamard-Zolésio structure theorem for shape derivatives, namely that deformations of the hold-all domain only influence the pre-shape functional if they deform the (pre-) shape ϕ(M ).
The difference to classical shape derivatives is illustrated in equation eq. (20), where the effect of deformations on the objective is split into normal and tangential components. The normal part g N can be understood as the shape optimization part of DJ, i.e. J depending on the change of interface ϕ(M ). This is also reflected by the structure of its support given in eq. (21), which states that only normal directions V ∈ N ϕ(M ) deforming ϕ(M ) have an effect on g N . On the other hand, g T is interpretable as the part of DJ being sensitive to reparameterization of the shape ϕ(M ), which is shown by the structure of its support in eq. (22), where only tangential vector field V ∈ T ϕ(M ) play a role. In classical shape optimization, tangential vectors are always in the nullspace of the shape derivative. But in the more general pre-shape case both components discussed can have non-trivial effects. This is also reflected by theorem 1 (iii), stating that pre-shape functionals having shape functionality have vanishing tangential part of the pre-shape derivative DJ, meaning that they are only supported by normal components of the deformation field V ∈ C ∞ (R n+1 , R n+1 ) just as in classical shape optimization theory. On the other hand, if 'shape derivatives' are not found to vanish in tangential directions, the 'shape functional' at hand is actually a true pre-shape functional. This is the case for mesh optimization techniques we will introduce.
Also, for shapes ϕ(M ) ⊂ R n+1 being bounded and topologically closed C k+1submanifolds of R n+1 with non-empty interior, the classical Hadamard-Zolésio structure theorem was generalized in [46,Thrm. 5.5]. In the special case of the objective J having shape functionality, i.e. vanishing tangential component of the pre-shape derivative as by theorem 1, the results of [46,Corollary 4.2] coincide with results in theorem 1.
Before we come to some exemplary pre-shape derivatives and their decompositions, we formulate a simple corollary, which connects the classical material derivatives to their pre-shape versions.
Then the material derivative decomposes to In particular, if f has shape functionality, then for the corresponding shape dependent functionf : holds for all ϕ ∈ Emb(M, R n+1 ), which means the pre-shape and classical material derivative coincide.
Proof. To get decomposition eq. (25), we simply use formula eq. (18) and apply decomposition eq. (20) of the structure theorem to the occurring pre-shape derivatives for fixed x 0 ∈ R n+1 . If f has shape functionality, we can apply part (iii) of theorem 1 to decomposition eq. (25) and see that the pre-shape material derivative equals the classical material derivative D mf (π(ϕ))[V ], since g N corresponds to the distribution from the classical Hadamard theorem.
With the last part of this corollary, we can apply pre-shape material derivatives to shape differentiable functions, yielding the same results as the classical material derivative by association with pre-shape extensions via proposition 1.
We have now seen that structure theorem 1, proposition 1 and corollary 1 guarantee validity of classical shape calculus formulae and results in the context of pre-shapes. Pre-shape calculus can be applied to objects from shape optimization if they are associated with their corresponding pre-shape counterparts, leading to the same derivatives and thus optimization methods. Even further, it is possible to apply pre-shape calculus to mixed shape and pre-shape problems, where the shape part is treated just as if shape calculus was applied, with the key difference that a pre-shape component would otherwise be non-accessible. In the following, we show some simple examples which are not accessible by classical shape calculus.
The pre-shape functional J measures the difference of a targetφ to another parameterized shape ϕ. Its pre-shape derivative can be calculated for directions V ∈ C ∞ (R 2 , R 2 ) by using elementary techniques We can choose S 1 with canonical parameterization as a starting pre-shape or point of reference by considering In order to formulate the Hadamard-Zolésio-type decomposition eq. (20), we need the outer normal vector field and an oriented tangential vector field, which for S 1 are given by Now we can examine the problem for several different parameterized target shapes ϕ ∈ Emb(S 1 , R 2 ). First, we can consider rescaling by a factor α ∈ (0, ∞), which lets S 1 contract or expand. The target for this is given bỹ Using eq. (28), the pre-shape derivative becomes This shows that rescaling of S 1 has vanishing parameterization part g T ≡ 0, whereas the remaining shape component g N is in the style of the classical Hadamard-Zolésio representation given above. In particular, only vector fields V acting in normal direction are supported. Next, let us consider a rotation of the circle. For this, we let α ∈ [0, 2π) and consider target rotations Plugging this into eq. (28) and doing some reformulations, the according decomposition becomes Here we see both components of the decomposition, the first corresponding to the normal g N , and the second corresponding to the tangential part g T . Notice that the normal component vanishes exactly for trivial rotations, whereas the tangential part also vanishes for the reflection at origin case α = π. Finally, we can also translate S 1 by some fixed z ∈ R 2 , which gives a target The decomposition of the pre-shape derivative becomes where decomposition into g N and g T depend on normal and tangential components of z on S 1 . Next, we give a summary of several useful pre-shape calculus formulae. Due to proposition 1 and corollary 1, they are also true for shape derivatives and functionals.
Corollary 2 (Pre-Shape Calculus Rules). Let f, g : Emb(M, R n+1 ) × R n+1 → R be pre-shape differentiable and differentiable in the second component, and let h : R → R be differentiable. Let Ω ⊆ R n+1 be an open, bounded domain with Lipschitz boundary, Γ an n-dimensional C ∞ -submanifold of R n+1 . Consider ϕ ∈ Emb(M, R n+1 ) and V ∈ C ∞ (R n+1 , R n+1 ). Then the following set of rules applies for pre-shape and -material derivatives, including the special case of shape derivatives.
with div Γ (V ) being the tangential divergence of V on ϕ, and Dh being the total derivative of h.
Proof. Let the assumptions stated above hold. Identity (i) was already discussed in corollary 1.
The product-and chain-rule (ii) and (iii) are simple consequences of the definition of the material derivative eq. (17).
Example 2: Pre-Shape Parameterization Tracking Problem With this problem class we introduce a non-trivial example for pre-shape optimization problems. Its pre-shape derivative in fact is no classical shape derivative, and thus not tractable by shape calculus. Later on, this problem class will enable us to optimize the overall mesh quality of discretizations and representations of shapes similar to deformations methods going back to [31]. To be specific, the fiber bundle structure of Emb(M, R n+1 ) offers the possibility to modify a given shape optimization problem in such a way, that the original solution is maintained, while at the same time optimization for the parameterization can take place. This gives rise to several different mesh regularization algorithms, and justifies the numerical optimization procedures we will establish in future works.
Before we further elaborate on this, we introduce necessary vocabulary and notation to formulate our problem in Emb(M, R n+1 ) for submanifolds M ⊂ R n+1 . First, we need the concept of local frames, which are local orthonormal bases of tangential vectors on M (cf. [30,Ch. 8]). For an open subset U ⊆ M , a smooth local frame is a tuple of dim(M ) tangential vector fields τ := (τ 1 , . . . , τ n ), such that for each p ∈ U the tangential vectors τ i (p) ∈ T p M are linearly independent. If we have a Riemannian metric on M , then we can additionally demand τ (p) = (τ 1 (p), . . . , τ n (p)) to be orthonormal with respect to the Riemannian metric for all p ∈ U , thus calling the frame (τ 1 , . . . , τ n ) local orthonormal frame. Note that local orthonormal frames always exist, due to simple use of the Gram-Schmidt algorithm in tangential spaces (cf. [30,Lemma 8.13]).
To achieve a natural and numerically tractable formulation of a pre-shape parameterization tracking problem, we also need to introduce the covariant derivative of an embedding ϕ ∈ Emb(M, R n+1 ). For this, we use similar ideas as in [12, Ch. 9.5.6] and modify them to our situation using local orthonormal frames. Given a ϕ ∈ Emb(M, R n+1 ), let τ : U → (T M ) n be a smooth local orthonormal frame on U ⊆ M and let τ ϕ : V → (T ϕ(M )) n be a local orthonormal frame on V ⊆ ϕ(M ). Without loss of generality, we can assume V = ϕ(U ), since we can choose V = V ∩ ϕ(U ). Then we define the local covariant derivative representation for ϕ under choice of frames τ and τ ϕ by . . . Dϕτ n,p , τ ϕ n,ϕ(p) where Dϕ is the Jacobian matrix of ϕ and ., . the Euclidean scalar product of R n+1 . We want to make clear that the covariant derivative D τ ϕ should not be mistaken for the tangential derivative D Γ ϕ, which is given by (cf. [12, Ch. 9.5.2]) Having numerical implementations in mind, we are also interested in the case of shapes with non-trivial boundaries. However, when the boundary is non-trivial, we impose restrictions on the pre-shapes permitted. Specifically, embeddings leaving the boundary invariant are sufficient, i.e.
For numerical routines this means that a specified boundary ∂M of the starting shape is left fixed, whereas the interior of the shape is able to deform and change its shape and parameterization. In the case of empty boundary, the pre-shape space becomes Emb(M, R n+1 ), meaning shapes are allowed to move freely. With the introduction of covariant derivatives and appropriate pre-shapes for the boundary case, we can formulate a pre-shape parameterization tracking problem inspired by a least-squares formulation of the deformation method for mesh element volume optimization as found in [8,10,21]. We remind the reader, that mesh deformation methods track for a specified target cell volume f by changing coordinates of nodes. For our formulation, we take a slight twist by using inverse Jacobians, which changes interpretation of optimal ϕ and targets f ϕ . In our case, f ϕ describes the desired local density of mesh nodes, whereas the authors of [8,10,21] use targets f describing the local cell volume. For this reason the mentioned authors need to use reciprocals of f , instead of reciprocals of Jacobians. Still, both formulations are equivalent by inverting the solutions, Jacobians and targets. In addition to targets f ϕ , we also incorporate a positive function g M : M → (0, ∞), which will act as the distribution of nodes for the initial mesh.
The following proposition gives the definition, well-definedness and existence of solutions of the pre-shape parameterization tracking problem.
Then the following problem is called pre-shape parameterization tracking problem. It is well-defined and independent of choice of local orthornormal frames τ , τ ϕ on M and ϕ(M ). Further, in each fiber π(ϕ) there exists a global C ∞ -solution to problem eq. (41), i.e. an embeddingφ satisfying Proof. The main ingredient of this proof is a theorem by Moser [36]  Fix an orientation for M , ϕ ∈ Emb ∂M (M, R n+1 ) and let τ , τ ϕ be local orthonormal frames for each ϕ ∈ Emb ∂M (M, R n+1 ). Well-definedness of the integrand in problem eq. (41) is clear, since g M is positive and ϕ is an immersion, making D τ ϕ ∈ GL(n, R) and thus det D τ ϕ non-vanishing. In the case of non-empty boundary, the integral is well-defined by using the interior of ϕ(M ), as the boundary ∂M is a set of measure zero. Independence of choice of orientation preserving local orthonormal bases inducing the covariant derivative (see eq. (37)) is also clear, since an orientation preserving change of the orthonormal base can be realized by multiplications with orthogonal matricesB, B ∈ SO(n), and hence by the determinant product rule det D τ ϕ remains invariant. Further, if no global orthonormal frame exists, well-definedness and independence of choice of local orthonormal frames is guaranteed as well. This can be ensured by using a partition of unity, which covers ϕ(M ) with open domains of local orthonormal frames, and linearity of the integral in eq. (41) together with the previous argument about the change of orthonormal bases.
Remark 3 (Hölder Regularity Case). Existence and well-definedness results from proposition 2 also hold in the more general context of C k,α -Hölder regularity. For given k ∈ N and α ∈ (0, 1), if M , f ϕ and g M have C k,α -regularity, and ∂M has C k+3,α regularity, then solutions ϕ with C k+1,α -regularity exist in each fiber. This stems from the regularity results in [11].
Having guaranteed existence of solutions, we want to turn our attention to the pre-shape derivative of the general parameterization tracking problem eq. (41). This is serves several different purposes in our studies. For one, it is of numerical interest, since we will construct several algorithms for improvement of mesh quality in shape optimization routines based on derivatives. At the same time, eq. (41) serves as a non-trivial example to illustrate the application of pre-shape calculus techniques developed in section 2. In particular, we will see that the derivative to the general parameterization tracking problem eq. (41) is not accessible via classical shape calculus techniques. In the following, we will leave the target density functions f ϕ general and only assume enough regularity for existence of the pre-shape derivative. Later in this article, we will propose an explicit way to construct f ϕ (cf. monitor function study [9]), while also ensuring existence of its material derivatives with a closed form (cf. eq. (62)).
As we permit non-empty boundaries left to be invariant in the parameterization tracking problem, the space of possible test functions V is altered in an according way. In particular, due to invariance of ∂M , vector fields for Emb ∂M (M, R n+1 ) are given by vector fields vanishing on the boundary (cf. [15,Thrm. 8.2], [45,Thrm. 3.19]), i.e.
Of course, the same is true for Hölder and Sobolev regularities. With this in mind, we can derive the pre-shape derivative of the parameterization tracking problem.
Theorem 2 (Pre-Shape Derivative of the Pre-Shape Parameterization Tracking Problem). Let the assumptions of proposition 2 hold and denote by J τ the objective functional of the general parameterization tracking problem eq. (41). Also, assume enough regularity for f ϕ , such that material derivatives exist. Then, for fixed ϕ ∈ Emb ∂M (M, R n+1 ) and V ∈ C ∞ ∂ϕ(M ) (R n+1 , R n+1 ), the preshape derivative of eq. (41) is given by with D m (f ϕ ) being the pre-shape material derivative of f ϕ and div Γ the tangential divergence on ϕ(M ). The pre-shape derivative does not depend on the choice of oriented local orthonormal frames τ, τ ϕ for representing the covariant derivative D τ .
Proof. For the proof we rely on pre-shape calculus rules we have established in section 2. In particular, we will make use of formulae found in corollary 2. So let M fulfill the assumptions made in theorem 2. Fix a ϕ ∈ Emb ∂M (M, R n+1 ) and let V ∈ C ∞ ∂ϕ(M ) (R n+1 , R n+1 ). The following arguments are all valid for C k,α -regularity. The use of pre-shape material derivative makes sense for families of differentiable functions on varying domains {f ϕ : ϕ(M ) → R} ϕ∈Emb ∂M (M,R n+1 ) depending smoothly on ϕ, since the limit D m f (ϕ) [V ] involves the term f ϕt (x t ). An easy check reveals that the term is well defined due to definition 4 of moving points x t and the perturbation of identity for pre-shapes eq. (12) coinciding (45) However, note that in this case there is no decomposition of type eq. (18) for D m f (ϕ) , since x t / ∈ ϕ(M ) and x 0 / ∈ ϕ t (M ) in general. With this in mind, we can apply corollary 2 (v) to eq. (41) in order to get For simplification of the material derivative of the integrand, we employ our assumption on existence of material derivatives for f ϕ , together with the chain-and product rule for material derivatives (cf. corollary 2), to see In the following we examine the remaining material derivatives in eq. (47), except for , as we let f ϕ remain general. To avoid confusion, we remind the reader that we are confronted with mappings h taking two arguments, one explicitly being a pre-shape, making them operators of the form We will use the following relationship of embeddings and the domain perturbation of identity where ϕ t is the perturbation of identity for pre-shapes (cf. eq. (12)) for sufficiently small t > 0. If material derivatives of h are assumed to exist, eq. (49) leads to the following elementary but interesting identity Applying this to the first remaining material derivative in eq. (47), we get since g M is does not depend on choice of ϕ ∈ Emb(M, R n+1 ). Next, we apply analogous techniques to the second material derivative. Hence, for calculation of the material derivative of det D τ ϕ • ϕ −1 it is sufficient to calculate its pre-shape derivative. Also, since the flow (ϕ t ) t∈[0,ε) given by the perturbation of identity in direction V (cf. eq. (12)) is differentiable in t, we can employ Jacobi's formula for the derivative of the determinant at t 0 = 0 to arrive at where Adju(·) is the adjugate matrix and tr(·) is the trace operator for matrices.
Knowing D τ ϕ is invertible for all p ∈ M due to ϕ ∈ Emb(M, R n+1 ), we can use Cramer's rule to express the adjugate in terms of inverses. Also, we can use invariance of the trace operator under permutations of multiplicative order of matrices, giving us Using eq. (51) and eq. (53) in eq. (47), and plugging the material derivative into eq. (46), we arrive at which is the desired pre-shape derivative eq. (44). The covariant derivative D τ ϕ, and hence also the pre-shape derivative eq. (44), is independent of choice of orthonormal frames by similar reasoning as in the first part of the proof to proposition 2.
In general cases, DJ τ (ϕ)[V ] is non-vanishing for vector fields V tangential to ϕ(M ). By structure theorem for pre-shape derivatives theorem 1, globally vanishing tangential pre-shape derivatives indicate a functional which is almost of classical shape functional type. If we take the form of the pre-shape derivative eq. (44), this clearly means eq. (41) cannot be formulated as a shape optimization problem, nor is it tractable by classical shape calculus.
In light of the main structure theorem 1 for pre-shape derivatives, we can further refine the representation of the pre-shape derivative eq. (44) by decomposing it into normal and tangential parts. Interestingly, if the user wants to optimize for mesh quality by using pre-shape derivative based parameterization tracking, it is not recommendable to use the full pre-shape derivative found in equation eq. (44). Instead, by decomposing the pre-shape derivative, we will see that the tangential component is sufficient for this task. If a special case of the pre-shape derivatives normal component is used, we recover numerical methods solving Plateau's problem by constructing minimal surfaces (cf. [13,38,14]). In some sense orthogonal to this, the use of tangential components gives algorithms resembling the deformation method for opimization of mesh quality (cf. [31,21,10]).
Decomposing the Pre-Shape Derivative To derive this decomposition, notice the following informal relationship between tangential divergence and the mean curvature κ for hypersurfaces for V ∈ C ∞ (R n+1 , R n+1 ) (cf. [29,Defi. 4 due to orthogonality of tangential gradients ∇ Γ ( V, n ) and the outer normal vector field n on the interior of ϕ(M ). Also, let us briefly assume f ϕ : R n+1 → R mapping from the whole ambient space, which simplifies using normal derivatives of f ϕ for the decomposition. With this, and the assumption of constant target parameterizations for each fiber, i.e. f ϕ = f ϕ•ρ for all ρ ∈ Diff(M ), we can refine eq. (44) to with shape (i.e. normal) component and pre-shape (i.e. tangential) component Here, D(f ϕ ) is the classical shape derivative of f ϕ . The first integral corresponds to the classical shape derivative component g N of decomposition eq. (20) acting on normal directions. Next, the second integral acts on tangential directions and therefore corresponds to the parameterization part g T in eq. (20).
Normal Component: Minimal Surfaces For illustration, let us deviate from normalization of the target (cf. eq. (40)) by chosing f ϕ = 0 and g M = 1 for all ϕ ∈ Emb ∂M (M, R n+1 ). In this situation the classical shape derivative component of DJ τ (ϕ) is given by Since ϕ are embeddings and M compact, according Jacobians are bounded and non-vanishing. In our special situation this means the horizontal component of the pre-shape derivative eq. (59) is vanishing exactly for shapes with vanishing mean curvature κ. Put differently, minimal surfaces and their higher dimension analogues are exactly the stationary points for this horizontal component. Hence to this observation, a gradient ascend using eq. (59) resembles an algorithm for evolutionary surfaces proposed by Dziuk in [13], which solves Plateau's problem by approximating a mean curvature flow. Note that an ascend is necessary, since our formulation of the pre-shape parameterization tracking problem involves inverse Jacobians, which is connected to Plateau's problem by (60) This also illustrates qualitative properties of a steepest descent using the complete pre-shape derivative. Briefly stated, if less vertices are desired at a location, the gradient descent in normal direction tends to blow up the shape, increasing distances of neighboring vertices. If more vertices are desired it tends to locally flatten the shape, driving the nodes together. This shows that application of the full pre-shape derivative distorts the shape in normal direction, hence interfering with the actual shape optimization problem to be regularized.
Use of the tangential component in eq. (56) leads to algorithms similar to the mesh deformation methods from [31,21,10]. Such a routine is discussed and implemented in the numerical section 3. We will use it in upcoming works to construct various regularization methods for shape optimization problems.
An illustrative numerical example of a pre-shape derivative for the parameterization tracking problem and its decomposition are shown in fig. 5 for target f ϕ (x, y, z) ≡ 1 ϕ(M ) 1 ds · x and a sphere centered at (0.5, 0.5, 0.5).
A class of externally defined targets f ϕ As we have left the target f ϕ unspecified during derivation of the pre-shape derivative eq. (44), we want to give an constructive example, which can be implemented in numerical routines. To do so, we have to keep in mind that the normalization requirement eq. (40) on f ϕ has to be fulfilled. One way to accomplish this, is by defining f ϕ using a given globally defined function q : R n+1 → (0, ∞). By assuming H 2 -regularity for q, existence of pre-shape derivatives and their closed form as in eq. (44) is guaranteed. The according target vertex point density on a shape ϕ(M ) is then given by which is well-defined due to the trace theorem for Sobolev functions.
If the target f ϕ is chosen such that normalization eq. (40) is not fulfilled, solutions to the parameterization tracking problem eq. (41) might still exist. Depending on whether M g M ds is greater or smaller M f ϕ ds, the gradient flow generated by the pre-shape derivative eq. (44) locally shrinks or blows up the shape ϕ(M ) in normal direction to compensate for the difference.
Next we calculate D m (f ϕ )[V ], which exists since we have q |ϕ(M ) ∈ H 1 ϕ(M ) . For this, we apply pre-shape calculus rules established in corollary 2. Also, since the external force q : R n+1 → (0, ∞) is defined on the entire ambient space, we can make direct use of decomposition for pre-shape material derivatives corollary 1. (c) Complete pre-shape gradient in volume mesh representation (d) Slice through xy-plane at center of complete volume pre-shape gradient Figure 5: Negative pre-shape gradient of J τ on a sphere scaled by 0.02 using target f ϕ (x, y, z) ≡ 1 ϕ(M ) 1 ds · x, which is depicted by color, and gradient representation by a linear elasticity metric. Color shifting towards red means higher desire fore more volume/vertex allocation.
Together with the fact that q and g M do not depend on ϕ, this gives where the last equality stems from Stokes theorem and our assumption on V to vanish on the boundary ∂ϕ(M ).
Next, we illustrate how the closed pre-shape derivative formula can be used to derive additional important properties of pre-shape optimization problems. In particular, we will see that local and global solutions parameterization tracking in each fiber coincide.
Proposition 3 (Characterization of Global Solutions by Fiber Stationarity). Let assumptions of theorem 2 be satisfied. Then the following statements are equivalent: where int ϕ(M ) is the interior of ϕ(M ) and n is the outer normal field on int ϕ(M ).
(ii) ϕ is a global solution to eq. (41), and in particular it satisfies (iii) the complete pre-shape derivative of J τ vanishes in ϕ, i.e.
Additionally, its normal component g N ϕ in ϕ vanishes. In particular, the necessary first order condition regarding only directions V tangential to ϕ(M ) is already a sufficient condition for being a global minimizer to J τ .
Proof. Let us assume the setting of theorem 2. We show equivalence of all assertions by a circular argument '(i) =⇒ (ii) =⇒ (iii) =⇒ (i) .
The non-trivial part is to prove '(i) =⇒ (ii)'. Let us assume (i) by fixing a ϕ ∈ Emb(M, R n+1 ) satisfying fiber stationarity eq. (63). With the pre-shape derivative formula from theorem 2 at hand, we can apply an integration by parts on manifolds (cf. [47, Ch. 2.2, Proposition 2.3]), either using that M is closed or V is vanishing on the boundary, to get where we also used the assumption in proposition 2 that f ϕ is constant in each fiber and Tr |ϕ(M ) [V ], n 2 = 0 to reformulate D m (f ϕ )[V ]. Due to assumption eq. (63) of fiber stationarity, we know eq. (66) equals zero for all V tangential on ϕ(M ) up to the boundary. So on the interior of ϕ(M ) we get By assumption we have g M > 0, so together with non-vanishing determinant by ϕ ∈ Emb(M, R n+1 ), this implies Since involved functions are Lipschitz continuous, and as proposition 3 assumes M to be smooth and path connected, we can derive constancy of the involved term. However, after using a pull-back and normalization assumption eq. (40), we see that the discussed constant is 0, giving us Since ϕ is chosen from Emb ∂M (M, R n+1 ), it leaves the boundary fixed, finally giving '(i) =⇒ (ii)'. Lastly, we easily see that '(iii) =⇒ (i)'. Since the complete vanishing of the pre-shape derivative eq. (65) in particular implies its vanishing for directions V tangential to ϕ(M ), which is eq. (63), and concludes the proof.
proposition 3 tells us that there are no stationary points other than global solutions to the pre-shape parameterization tracking problem eq. (41). This strongly resembles the situation for convex optimization problems, where the only candidates for local optimality are indeed global solutions.
Notice that proposition 3 gives existence of stationary points ϕ for each shape via existence result proposition 2. Since stationary points are global solutions, we can simply use the existence result proposition 2 for this.
Additionally, proposition 3 guarantees that optimization with the tangential component of pre-shape derivative eq. (44) is sufficient to reach a globally optimal solution for (41). This permits design of regularizations for shape optimization algorithms using pre-shape parameterization tracking with the property to leave the shape at hand invariant, while at the same time finding an optimal parameterization of the respective shape.

Numerical Tests of Parameterization Tracking involving Pre-Shape Derivatives
We have now finished our introduction of pre-shape calculus and its application to parameterization tracking problems. In order to test our theoretical results, we present three implementations of pre-shape gradient descent methods for the parameterization tracking problem. For this we use the open-source finite-element software FEniCS (cf. [33,2]). Construction of meshes is done via the free meshing software Gmsh (cf. [20]). We use a single core of an Intel(R) Core(Tm) i3-8100 CPU at 3.60 GHz featuring 16 GB RAM. The single core runs at 800 MHz while the code is executed on a virtual machine.
In the following, we show three implementations solving the parameterization tracking problem eq. (41) by using the tangential component of the pre-shape derivative seen in decomposition eq. (56). The solution process also features a simple backtracking line search, which scales the initial gradient of the current iteration U i according to a given factor c and rescales it by 0.5 if no sufficient decrease in J is apparent. In order to apply a descent algorithm, we are in need of pre-shape gradients. Because gradients are defined with respect to a sufficient bilinear form, we have to choose a form which fits our application. Since we are in infinite dimensions, there is a multitude of non-equivalent choices to represent derivatives as gradients. These can differ in resulting regularity of the gradients, and also in computational expense. As a bilinear form, we choose the weak formulation of the linear elasticity as proposed in [43], which gives us H 1 -regularity of pre-shape gradients. By only using the shear component of the linear elasticity featuring the second Lamé parameter µ elas , and adding a zero order term, the gradient U is calculated by solving its representing system Here, we choose the weights α LE , α L 2 > 0. For µ max , µ min > 0, the second Lamé parameter µ elas is chosen as the solution of the Poisson problem Solving eq. (70) on the entire hold-all domain D gives us a volume representation U of the pre-shape derivative DJ τ . The pre-shape gradient system eq. (70) is assembled in FEniCS and solved with a sparse LU method from PETSc used as a linear algebra backend. The first example shows an application of the parameterization tracking problem to improve the quality of a given hold-all domain D = [0, 1] 2 ⊂ R 2 . This is realized by using an unstructured 2-dimensional volume mesh created via Gmsh featuring 4262 triangular cells and 2212 nodes. Then we distort the mesh quality of this unstructured mesh by applying as a deformation to the interior of D. The deformed initial mesh ϕ 0 (D) is depicted in fig. 6. Notice that in this scenario the initial model M is given by the hold-all domain D = [0, 1] 2 with non-trivial boundary ∂D. Therefore we are in the situation where the boundary ∂D is left invariant (cf. eq. (39)). Also, there is no normal component of the pre-shape derivative in this case, as the codimension of D ⊂ R 2 is zero.
To formulate the parameterization tracking problem eq. (41), we need to specify an initial point distribution g M and target f ϕ . Here, the target chosen is given by the constant This ensures that a uniform cell volume distribution of the hold-all domain is targeted. The initial point distribution g M is represented by using a continuous Galerkin Ansatz with linear elements. Degrees of freedom are situated at the mesh vertices and set to the average of inverses of surrounding cell volumes, i.e.
Here p i is a mesh vertex and C i is the set of its neighboring cells C. Finally, the resulting function is normed to satisfy the demanded normalization condition eq. (40) of the parameterization tracking problem. The initial point distribution estimated by this procedure is shown in fig. 6 (a).
With both g M and f ϕ specified, the target J τ and its pre-shape derivative DJ τ can be assembled. For the gradient representation we use weights α LE = 0.02, α L 2 = 1 and Lamé parameters µ max = µ min = 1, resulting in constant µ elas = 1. An initial scaling factor of c = 0.01 for the negative gradient during line search is applied. The method successfully exits after 37.07s and 45 iterations. Results of the pre-shape gradient descent using the tangential component of the pre-shape derivative and the described methodology are shown in fig. 6 and fig. 8.
As our second example, we use the exact same parameters as in the first example. Note that in particular, the starting mesh and therefore its initial volume distribution g M are the same as in the first example. We can see an illustration in fig. 6 (a). To show the general applicability of parameterization tracking, we replace the uniform target f ϕ from eq. (73) by a more complicated non-uniform target (75) The pre-shape gradient descent for this non-uniform target achieves convergence after 38.12s and 46 iterations. We visualize an intermediate mesh, and the final mesh in fig. 6 (c) and (d). The target function values J τ (ϕ i ) and pre-shape gradient norms are shown in fig. 8. Interestingly, notice that the intermediate mesh (c) looks like a superposition of the final and initial meshes (d) and (a). Essentially, this is an illustration of snapshots from a discretized flow in the fiber of Emb(D, D) corresponding to the shape D, which is abstractly visualized in fig. 2. We see in fig. 6 (d), that the prescribed non-uniform cell volume distribution is achieved, even though the initial mesh in fig. 6 (d) has degenerate cells distributed on vertical lines. Our third example applies the parameterization tracking problem to a sphere centered in the hold-all domain D = [0, 1] 3 ⊂ R 3 . It acts as the modeling manifold M and its initial parameterization ϕ 0 is given by the identity embedding it into the hold-all domain. The initial shape is a structured triangular surface mesh approximating a sphere centered in (0.5, 0.5, 0.5) with radius 0.3 using Gmsh. It consists of 6240 triangular cells and 3122 vertices on the surface. The sphere is embedded in a hold-all domain consisting of 21838 tetraedic cells and 4059 nodes.
For the third example we target a non-uniform surface cell volume distribution given by · sin(10 · 2π · x).
The target function is of the form eq. (61), which permits use of the material derivative formula eq. (62) for assembling the pre-shape derivative DJ τ . At the same time, it satisfies normalization condition eq. (40). Also, we set the initial vertex distribution to a constant In order to calculate covariant derivatives and associated Jacobian determinants, we apply a Gram-Schmidt algorithm to construct local tangential orthonormal bases.
Here, we choose weights α LE = 0.02, α L 2 = 1 and Lamé parameters µ max = 30, µ min = 5 for gradient representation. The line search employs an initial scaling factor c = 0.001 for the negative gradient. For this scenario, the gradient representation of the pre-shape derivative, and the resulting surface mesh with its associated vertex distribution are depicted in fig. 7. The method successfully exits after 1256.78s and 48 iterations. Target function values J τ (ϕ i ) and pre-shape gradient norms are shown in fig. 8. In light of proposition 3, we see that the gradient norm and target values converge simultaneously by using tangential components of DJ τ only. Also, the shape of the sphere is left invariant, which would not be the case if normal components or the full pre-shape derivative (cf. fig. 5) were used.

Conclusion and Outlook
In this work we introduced a unified framework to formulate shape optimization and mesh quality optimization problems. A calculus for pre-shape derivatives, which act in normal and tangential directions, and according structure theorems were derived. In particular, rules and problem formulations from classical shape optimization carry over to the pre-shape setting. These techniques were tested on a class of parameterization tracking problems. Resulting numerical implementations of a gradient descent method based on decomposed pre-shape derivatives show promising results for optimization of volume-and surface mesh quality. In forthcoming works we will derive efficient algorithms harnessing the opportunity to simultaneously solve shape optimization problems and improve mesh quality of shapes and ambient spaces. For this, can will design various targets for parameterization tracking, giving a desired type of mesh for the user. Also, we will define pre-shape Hessians to harness second order information. (Quasi-)Newton methods in the context of pre-shape optimization, as well as an optimal choice of pre-shape gradient representations, will enhance the performance of these algorithms.  (a) J τ (ϕi) (b) Ui L 2 (D,R n+1 ) Figure 8: (a) Values for the pre-shape parameterization tracking target J τ (ϕ i ) for iterates ϕ i of the tangential pre-shape derivative component based steepest descent method. Target for the 3D sphere case is scaled by 3.
(b) L 2 -norms U i L 2 (D,R n+1 ) of the gradient representations U i of pre-shape derivatives for each iterate ϕ i . Gradient norms for the 3D sphere case are scaled by 25.