Simultaneous Shape and Mesh Quality Optimization using Pre-Shape Calculus

Throughout section 2, we take a look at a general prototype shape optimization problem

$$\underset{\Gamma\in B_{e}^{n}}{\min}\mathcal{J}(\Gamma).$$

(1)

Here we only assume that the shape functional $\mathcal{J}:B_{e}^{n}\rightarrow{\mathbb{R}}$ is first order shape differentiable.

Now we reformulate eq. 1 in the context of pre-shape optimization by use of the canonical projection $\pi:\operatorname{Emb}(M,{\mathbb{R}}^{n+1})\rightarrow B_{e}^{n}$. The canonical projection $\pi$ maps each pre-shape $\varphi\in\operatorname{Emb}(M,{\mathbb{R}}^{n+1})$ to an equivalence class $\pi(\varphi)=\Gamma\in B_{e}^{n}$, which consists of all pre-shapes mapping onto the same image shape in ${\mathbb{R}}^{n+1}$. We remind the reader that we can associate $\pi(\varphi)$ with the set interpretation $\varphi(M)$ of the shape $\Gamma$. So $\pi(\varphi)$ can be interpreted as the set of all parameterizations of a given shape $\Gamma\in B_{e}^{n}$ in ${\mathbb{R}}^{n+1}$. The pre-shape formulation of eq. 1 takes the form

$$\underset{\varphi\in\operatorname{Emb}(M,{\mathbb{R}}^{n+1})}{\min}(\mathcal{J}\circ\pi)(\varphi).$$

(2)

It is important to notice that the extended target functional of eq. 2 is pre-shape differentiable in the sense of [12, Definition 3], since we assumed $\mathcal{J}$ to be shape differentiable (cf. [12, Prop. 1]).

Next, we need the pre-shape parameterization tracking functional introduced in [12, Prop. 2]. In the discrete setting, this functional allows us to track for given desired sizes of surface and volume mesh cells. For this, let us assume functions $g^{M}:M\rightarrow(0,\infty)$ and $f^{M}_{\varphi}:\varphi(M)\rightarrow(0,\infty)$ to be smooth and fulfill the normalization condition

$$\int_{\varphi(M)}f^{M}_{\varphi}(s)\;\mathrm{d}s=\int_{M}g^{M}(s)\;\mathrm{d}s\quad\forall\varphi\in\operatorname{Emb}(M,{\mathbb{R}}^{n+1}).$$

(3)

Further, we assume $f$ to have shape functionality (cf. [12, Definition 2]), i.e.

$$f^{M}_{\varphi\circ\rho}=f^{M}_{\varphi}\quad\forall\rho\in\operatorname{Diff}(M).$$

(4)

Shape functionality for $f^{M}$ simply means, that $f^{M}_{\varphi_{1}}=f^{M}_{\varphi_{2}}$ for all pre-shapes $\varphi_{1},\varphi_{2}\in\pi(\varphi)$ corresponding to the same shape $\Gamma=\pi(\varphi)$.

Finally, the pre-shape parameterization tracking problem takes the form

$$\underset{\varphi\in\operatorname{Emb}(M,{\mathbb{R}}^{n+1})}{\min}\frac{1}{2}\int_{\varphi(M)}\Big{(}g^{M}\circ\varphi^{-1}(s)\cdot\det D^{\tau}\varphi^{-1}(s)-f^{M}_{\varphi}(s)\Big{)}^{2}\;\mathrm{d}s=:\mathfrak{J}^{\tau}(\varphi)$$

(5)

where $D^{\tau}$ is the covariant derivative. For each shape $\Gamma\in B_{e}^{n}$ there exists a pre-shape $\varphi\in\pi(\varphi)=\Gamma$ minimizing eq. 5, which is guaranteed by [12, Prop. 2].

As we need the pre-shape derivative of the parameterization tracking functional $\mathfrak{J}^{\tau}$ in order to devise our algorithms, we state it in the style of the pre-shape derivative structure theorem [12, Thrm. 1]

$$\mathfrak{D}\mathfrak{J}^{\tau}(\varphi)[V]=\langle g_{\varphi}^{\mathcal{N}},V\rangle+\langle g_{\varphi}^{\mathcal{T}},V\rangle\quad\forall V\in C^{\infty}({\mathbb{R}}^{n+1},{\mathbb{R}}^{n+1}),$$

(6)

with normal/ shape component

$\displaystyle\begin{split}\langle g_{\varphi}^{\mathcal{N}},V\rangle=&-\int_{\varphi(M)}\frac{\operatorname{dim}(M)}{2}\cdot\Big{(}\big{(}g^{M}\circ\varphi^{-1}\cdot\det D^{\tau}\varphi^{-1}\big{)}^{2}-f_{\varphi}^{2}\Big{)}\cdot\kappa\cdot\langle V,n\rangle\\ &\qquad\qquad\quad+\Big{(}g^{M}\circ\varphi^{-1}\cdot\det D^{\tau}\varphi^{-1}-f_{\varphi}\Big{)}\cdot\Big{(}\frac{\partial f_{\varphi}}{\partial n}\cdot\langle V,n\rangle+\mathcal{D}(f_{\varphi})[V]\Big{)}\;\mathrm{d}s\end{split}$

(7)

and tangential/ pre-shape component

$\displaystyle\begin{split}\langle g_{\varphi}^{\mathcal{T}},V\rangle=&-\int_{\varphi(M)}\frac{1}{2}\cdot\Big{(}\big{(}g^{M}\circ\varphi^{-1}\cdot\det D^{\tau}\varphi^{-1}\big{)}^{2}-f_{\varphi}^{2}\Big{)}\cdot\operatorname{div}_{\Gamma}(V-\langle V,n\rangle\cdot n)\\ &\qquad\qquad\quad+\Big{(}g^{M}\circ\varphi^{-1}\cdot\det D^{\tau}\varphi^{-1}-f_{\varphi}\Big{)}\cdot\nabla_{\Gamma}f_{\varphi}^{T}V\;\mathrm{d}s.\end{split}$

(8)

Notice that we signify a duality pairing by $\langle\cdot,\cdot\rangle$, and with $\mathcal{D}(f_{\varphi})[V]$ the classical shape derivative of $f_{\varphi}$ in direction $V$. For the derivation of the pre-shape derivative $\mathfrak{D}\mathfrak{J}^{\tau}$ by use of pre-shape calculus, we refer the reader to [12, Thrm. 2, Cor. 2].

All ingredients necessary are now available to formulate a regularized version of a shape optimization problem. For $\alpha^{\tau}>0$, we add the parameterization tracking functional in style of a regularizing term to pre-shape reformulation eq. 2

$$\underset{\varphi\in\operatorname{Emb}(M,{\mathbb{R}}^{n+1})}{\min}(\mathcal{J}\circ\pi)(\varphi)+\alpha^{\tau}\cdot\mathfrak{J}^{\tau}(\varphi).$$

(9)

We already found that the pre-shape extension eq. 2 of our initial shape optimization problem is pre-shape differentiable. Hence eq. 9 is pre-shape differentiable as well. This is foundational, since a pre-shape gradient system for the regularized problem eq. 9 can always be formulated.

The pre-shape derivative $\mathfrak{D}\mathfrak{J}(\varphi)[V]$ of the parameterization tracking problem (cf. eq. 6) is well-defined for weakly differentiable directions $V\in H^{1}({\mathbb{R}}^{n+1},{\mathbb{R}}^{n+1})$. By assuming the same for the shape derivative $\mathcal{D}\mathcal{J}$ of the original problem eq. 1, we can create a pre-shape gradient system for eq. 9 using a weak formulation with $H^{1}$-functions. Given a symmetric, positive-definite bilinear form $a(.,.)$ (e.g. [17], [13]), such a system takes the form

$$a(U^{\mathcal{J}+\mathfrak{J}^{\tau}},V)=\mathcal{D}\mathcal{J}(\Gamma)[V]\;+\;\alpha^{\tau}\cdot\langle g_{\varphi}^{\mathcal{N}},V\rangle\;+\;\alpha^{\tau}\cdot\langle g_{\varphi}^{\mathcal{T}},V\rangle\quad\forall V\in H^{1}({\mathbb{R}}^{n+1},{\mathbb{R}}^{n+1}).$$

(10)

With $\Gamma=\pi(\varphi)$, the right hand side of eq. 10 is indeed the full pre-shape gradient of eq. 9. This stems from the fact that the pre-shape extension $\mathcal{J}\circ\pi$ has shape functionality by construction, which makes the pre-shape derivative $\mathfrak{D}(\mathcal{J}\circ\pi)$ equal to the shape derivative $\mathcal{D}\mathcal{J}$ by application of [12, Thrm. 1 $(iii)$].

The full pre-shape gradient system eq. 10 already achieves simultaneous solution of the shape optimization problem and regularization of shape mesh quality. Namely, it is not required to additionally solve linear systems to create a mesh quality regularized descent direction $U^{\mathcal{J}+\mathfrak{J}^{\tau}}$, since the original shape gradient system to problem eq. 1 is modified by adding the two force terms $g^{\mathcal{N}}$ and $g^{\mathcal{T}}$. A calculation of a descent direction $U^{\mathcal{J}+\mathfrak{J}^{\tau}}=U^{\mathcal{J}}+U^{\mathfrak{J}^{\tau}}$ by combining the shape gradient $U^{\mathcal{J}}$ and pre-shape gradient $U^{\mathfrak{J}^{\tau}}$ solving the decoupled systems

$$a(U^{\mathcal{J}},V)=\mathcal{D}\mathcal{J}(\Gamma)[V]\quad\forall V\in H^{1}({\mathbb{R}}^{n+1},{\mathbb{R}}^{n+1})$$

(11)

and

$$a(U^{\mathfrak{J}^{\tau}},V)=\alpha^{\tau}\cdot\langle g_{\varphi}^{\mathcal{N}},V\rangle+\alpha^{\tau}\cdot\langle g_{\varphi}^{\mathcal{T}},V\rangle\quad\forall V\in H^{1}({\mathbb{R}}^{n+1},{\mathbb{R}}^{n+1})$$

(12)

is in fact equivalent to solution of eq. 10.

However, by using eq. 10 for gradient based optimization, the first demanded property of leaving the optimal or intermediate shapes invariant is not true in general. This issue comes from involvement of the shape component $g^{\mathcal{N}}_{\varphi}$ (cf. eq. 7) of the pre-shape derivative to the parameterization tracking functional $\mathfrak{J}^{\tau}$. It is acting solely in normal directions, therefore altering shapes by interfering with the shape derivative $\mathcal{D}\mathcal{J}$ of the original problem.

For this reason, we deviate from the full gradient system eq. 10 consistent with the regularized problem eq. 9 by using a modified system

$$a(\tilde{U},V)=\mathcal{D}\mathcal{J}(\Gamma)[V]+\alpha^{\tau}\cdot\langle g_{\varphi}^{\mathcal{T}},V\rangle\quad\forall V\in H^{1}({\mathbb{R}}^{n+1},{\mathbb{R}}^{n+1}).$$

(13)

We project the pre-shape derivative $\mathfrak{D}\mathfrak{J}^{\tau}$ onto its tangential part, which is realized by simply removing the shape component $g^{\mathcal{N}}$ from the right hand side of the gradient system. By this we still have numerical feasibility, while also achieving invariance of optimal and intermediate shapes. Of course invariance is only true up to discretization errors. From a traditional shape optimization perspective, this stems from considering eq. 13 as a shape gradient system with additional force term acting on directions $V$ in the kernel of the classical shape derivative $\mathcal{D}\mathcal{J}$. In particular, by Hadamard’s theorem or the structure theorem for pre-shape derivatives [12, Thrm. 1], directions tangential to shapes $\Gamma$ are in the kernel of $\mathcal{D}\mathcal{J}$. Using the pre-shape setting, we interpret these directions as vector fields corresponding to the fiber components of $\operatorname{Emb}(M,{\mathbb{R}}^{n+1})$.

To sum up and justify the use of the pre-shape regularized problem eq. 9 and its modified gradient system eq. 13, we provide existence of solutions and consistency of the modified gradients with the regularized problem.

With theorem 1 we can rest assured that optimization algorithms using the tangentially regularized gradient $U$ from eq. 13 leave stationary points of the original problem invariant. Also, vanishing of the modified gradient $\tilde{U}$ indicates that we have a stationary shape $\Gamma=\pi(\varphi)$ with parameterization $\varphi$ having desired cell volume allocation $f_{\varphi}$. Of course, this is only true up to discretization error. We also see that the modified gradient system eq. 13 captures the same information as the shape gradient system eq. 11 and full pre-shape gradient system eq. 10 combined. This might seem counterintuitive at first glance, especially since necessary information is contained in one instead of two gradient systems. However, by application of pre-shape calculus to derive the fiber stationarity characterization of solutions [12, Prop. 3], we recognize this circumstance as a consequence of the special structure of regularizing functional $\mathfrak{J}^{\tau}$. The fact that pre-shape spaces are fiber spaces with locally orthogonal tangential bundles for parameterizations and shapes (cf. [12, Section 2]) is a fundamental prerequisite to this. We will discuss numerical results comparing optimization using standard shape gradients eq. 11 and gradients regularized by modified shape parameterization tracking eq. 13 in section 2.

Using pre-shape techniques to increase volume mesh quality is a different task than regularizing the shape mesh as we have seen in section 2.1. Using a pre-shape space $\operatorname{Emb}(M,{\mathbb{R}}^{n+1})$ for this endeavor does not make sense, since this space contains all shapes in ${\mathbb{R}}^{n+1}$ combined with their parameterization. For $\operatorname{Emb}(M,{\mathbb{R}}^{n+1})$, parameterizations are described via the $n$-dimensional modeling or parameter manifold $M$. We now need a second, different parameter space modeling the hold-all domain.

For this, let $\mathbb{D}\subset{\mathbb{R}}^{n+1}$ be a compact, orientable, path-connected, $n+1$-dimensional $C^{\infty}$-manifold. This hold-all domain will replace the unbounded hold-all domain ${\mathbb{R}}^{n+1}$ of our previous discussion. In most numerical cases, the hold-all domain has a non-empty boundary $\partial\mathbb{D}$. Hence we also permit non-empty $\partial\mathbb{D}$, but demand $C^{\infty}$-regularity of the boundary. We remind the reader that $\mathbb{D}$ is a compact manifold with boundary, so we have to distinguish from its interior, which we denote by $\operatorname{int}(\mathbb{D})$.

A first choice for a pre-shape space to the hold-all domain would be $\operatorname{Emb}_{\partial\mathbb{D}}(\mathbb{D},\mathbb{D})$, which is the set of all $C^{\infty}$-embeddings leaving $\partial\mathbb{D}$ pointwise invariant, i.e.

$$\operatorname{Emb}_{\partial\mathbb{D}}(\mathbb{D},\mathbb{D}):=\{\varphi\in\operatorname{Emb}(\mathbb{D},\mathbb{D}):\varphi(p)=p\quad\forall p\in\partial\mathbb{D}\}.$$

(21)

Leaving the outer boundary $\partial\mathbb{D}$ invariant can be beneficial from practical standpoint, since it might be part of an underlying problem formulation, such a boundary conditions of a PDE.

The pre-shape space eq. 21 in fact is the diffeomorphism group $\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D})$ leaving the boundary pointwise invariant. This is very important to notice, because it means $\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D})$ is a pre-shape space consisting exactly of one fiber. The shape of hold-all domain $\mathbb{D}$ is fixed, and as a consequence its tangent bundle $T\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D})$ consists of vector fields with vanishing normal components on $\partial\mathbb{D}$ (cf. [20, Thrm. 3.18]), i.e.

$$T_{\Phi}\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D})\cong C^{\infty}_{\partial\mathbb{D}}(\mathbb{D},{\mathbb{R}}^{n+1}):=\{V\in C^{\infty}(\mathbb{D},{\mathbb{R}}^{n+1}):\operatorname{Tr}_{|\partial\mathbb{D}}(V)=0\}\quad\forall\Phi\in\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D}).$$

(22)

Here, $\operatorname{Tr}_{|\partial\mathbb{D}}(V)$ is the trace of $V$ on $\partial\mathbb{D}$. Even further, the structure theorem for pre-shape derivatives [12, Thrm. 1] tells us that a pre-shape differentiable functional $\mathfrak{J}^{\mathbb{D}}:\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D})\rightarrow{\mathbb{R}}$ always has trivial shape component $g^{\mathcal{N}}=0$ of $\mathfrak{D}\mathfrak{J}$. In particular, ordinary shape calculus is not applicable for functionals of this type.

Our task to design mesh regularization strategies, which are numerically feasible and do not interfere with the original shape optimization problem. For a given shape $\varphi(M)\subset\mathbb{D}$, the latter is not guaranteed when a pre-shape space $\operatorname{Diff}_{\partial\mathbb{D}}$ is used to model parameterizations of the hold-all domain. A diffeomorphism $\Phi\in\operatorname{Diff}_{\partial\mathbb{D}}$ could change a given intermediate or optimal shape $\varphi(M)$, i.e. $\Phi\big{(}\varphi(M)\big{)}\neq\varphi(M)$ in general. Therefore we have to enforce invariance of shapes by further restricting the pre-shape space for $\mathbb{D}$. For this, we use the concept of diffeomorphism groups leaving submanifolds invariant (cf. [20, Ch. 3.6]). With our previous notation eq. 21, the final pre-shape space is $\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$, for a given fixed shape $\varphi(M)\subset\mathbb{D}$. This space is a subgroup of $\operatorname{Diff}_{\partial\mathbb{D}}(\mathbb{D})$ and $\operatorname{Diff}(\mathbb{D})$, with a tangential bundle consisting of vector fields vanishing both on $\partial\mathbb{D}$ and $\varphi(M)$ (cf. [20, Thrm. 3.18]).

Now we have a pre-shape space suitable to model reparameterizations of the hold-all domain $\mathbb{D}$, which can leave a given shape $\varphi(M)$ invariant. Next, we formulate an analogue of parameterization tracking problem eq. 9, which is tracking for the volume parameterization of the hold-all domain. Let us fix a $\varphi\in\operatorname{Emb}(M,\mathbb{D})$ and the corresponding shape $\varphi(M)\subset\mathbb{D}$. The volume parameterization tracking problem is then

$$\underset{\Phi\in\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})}{\min}\frac{1}{2}\int_{\mathbb{D}}\Big{(}g^{\mathbb{D}}\circ\Phi^{-1}(s)\cdot\det D\Phi^{-1}(s)-f^{\mathbb{D}}_{\varphi(M)}(s)\Big{)}^{2}\;\mathrm{d}x=:\mathfrak{J}^{\mathbb{D}}(\Phi).$$

(23)

Notice the similarities and differences to shape parameterization tracking problem eq. 9. Both pre-shape functionals track for a parameterization dictated by target $f$, and both feature a function $g$ describing the initial parameterization. The volume tracking functional $\mathfrak{J}^{\mathbb{D}}$ differs from $\mathfrak{J}^{\tau}$ by featuring a volume integral, instead of a surface one. Also, the covariant derivative of the Jacobian determinant in $\mathfrak{J}^{\mathbb{D}}$ is just the regular Jacobian matrix of $\Phi^{-1}$. The two most important differences concern their sets of feasible solutions and targets $f$. The target $f^{\mathbb{D}}_{\varphi(M)}$ does not depend on the pre-shape $\Phi$ for the hold-all domain, but instead depends on the shape $\varphi(M)$ which is left to be invariant. Having $f^{\mathbb{D}}$ depend on the shape of $\mathbb{D}$ does not make sense, since $\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$ consists only of one fiber, i.e. the shape of $\mathbb{D}$ remains invariant. Hence there is a dependence of both the feasible set of pre-shapes $\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$ and the target $f^{\mathbb{D}}_{\varphi(M)}$ on the shape $\varphi(M)$, because we desired $\varphi(M)$ to stay unaltered. For this reason, the volume parameterization tracking problem eq. 23 differs from eq. 9 to such an extent, that [12, Prop. 2] does not cover existence of solutions $\Phi\in\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$ for eq. 23. This makes it necessary to formulate a result guaranteeing existence of solutions for eq. 23 under appropriate conditions.

As we want to regularize gradient descent methods for shape optimization, we need to explicitly specify the pre-shape derivative to the volume parameterization tracking problem eq. 23. However, it is also of theoretical interest, since it serves as an example of a problem with a pre-shape derivative having vanishing shape component. This means that neither the formulation of eq. 23 in the context of classical shape optimization, nor the derivation of a derivative using classical shape calculus are possible. Since the form of $\mathfrak{J}^{\mathbb{D}}$ is similar to $\mathfrak{J}^{\tau}$, we can mimic the steps from [12, Thrm. 2], which we do not restate for brevity. Then, for given $\Phi\in\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$, the pre-shape derivative of $\mathcal{J}^{\mathbb{D}}$ in decomposed form (cf. [12, Thrm. 1]) is given by

$$\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[V]=\langle g_{\Phi}^{\mathcal{N}},V\rangle+\langle g_{\Phi}^{\mathcal{T}},V\rangle\quad\forall V\in C_{\partial M\cup\varphi(M)}^{\infty}(\mathbb{D},{\mathbb{R}}^{n+1}),$$

(28)

with normal/ shape component

$\displaystyle\begin{split}\langle g_{\Phi}^{\mathcal{N}},V\rangle=0\end{split}$

(29)

and tangential/ pre-shape component

$\displaystyle\begin{split}\langle g_{\Phi}^{\mathcal{T}},V\rangle=&-\int_{\mathbb{D}}\frac{1}{2}\cdot\Big{(}\big{(}g^{\mathbb{D}}\circ\Phi^{-1}\cdot\det D\Phi^{-1}\big{)}^{2}-{f^{\mathbb{D}}_{\varphi(M)}}^{2}\Big{)}\cdot\operatorname{div}(V)\\ &\qquad\qquad\quad+\Big{(}g^{\mathbb{D}}\circ\Phi^{-1}\cdot\det D\Phi^{-1}-f^{\mathbb{D}}_{\varphi(M)}\Big{)}\cdot\mathfrak{D}_{m}\big{(}f^{\mathbb{D}}_{\varphi(M)}\big{)}[V]\;\mathrm{d}x.\end{split}$

(30)

Here, $n$ is the outer unit normal vector field of the invariant submanifold $\varphi(M)$. As previously mentioned, we signify a duality pairing by $\langle\cdot,\cdot\rangle$. It is also important to see that the restriction of $\operatorname{Diff}(\mathbb{D})$ to $\operatorname{Diff}_{\partial M\cup\varphi(M)}(\mathbb{D})$ does not influence the form of pre-shape derivative $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}$. Instead, it influences the space of applicable directions $V$ by restricting $C^{\infty}(\mathbb{D},{\mathbb{R}}^{n+1})$ to $C^{\infty}_{\partial M\cup\varphi(M)}(\mathbb{D},{\mathbb{R}}^{n+1})$. This stems from the relationship of tangential bundles of diffeomorphism groups with invariant submanifolds eq. 22. There is also another subtle difference of $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}$ and $\mathfrak{D}\mathfrak{J}^{\tau}$. Namely, the featured pre-shape material derivative $\mathfrak{D}_{m}\big{(}f^{\mathbb{D}}_{\varphi(M)}\big{)}$ depends on the invariant shape $\varphi(M)\subset\mathbb{D}$ instead of the volume pre-shape $\Phi$.

It is straight forward to formulate a pre-shape gradient system for $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}$ in the style of section 2.1 using Sobolev functions. For a symmetric, positive-definite bilinear form $a(.,.)$ it takes the form

$$a(U^{\mathfrak{J}^{\mathbb{D}}},V)=\alpha^{\mathbb{D}}\cdot\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[V]\quad\forall V\in H^{1}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D},{\mathbb{R}}^{n+1}).$$

(31)

Notice that the space of test functions enforces vanishing $U^{\mathfrak{J}^{\mathbb{D}}}$ on $\varphi(M)$. With a pre-shape gradient $U^{\mathfrak{J}^{\mathbb{D}}}$, it is possible to optimize for the volume parameterization tracking problem eq. 23 without altering the shape of $\varphi(M)$.

At this point, we have provided a suitable space $\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$ for pre-shapes representing the parameterization of hold-all domains $\mathbb{D}$, which leave a given shape $\varphi(M)\subset\mathbb{D}$ invariant. Also, we are able to guarantee existence for global minimizers to the volume version of parameterization tracking problem eq. 23 for all shapes $\varphi(M)$. In this subsection we are going to incorporate the volume mesh quality regularization simultaneously with shape optimization and shape mesh quality regularization.

To formulate a regularized version of the original shape optimization problem eq. 1, we have to keep the different types of pre-shapes involved in mind. These pre-shapes $\varphi\in\operatorname{Emb}(M,\mathbb{D})$ and $\Phi\in\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$ correspond to completely different shapes $\varphi(M)$ and $\mathbb{D}$. This is also illustrated by looking at the pre-shapes as maps

$$\varphi:M\rightarrow\mathbb{D}\quad\text{and}\quad\Phi:\mathbb{D}\rightarrow\mathbb{D}.$$

(32)

For this reason, we cannot simply proceed by adding $\mathfrak{J}^{\mathbb{D}}$ in style of a regularizer to increase volume mesh quality. From the viewpoint of problem formulation, this signifies a main difference in application of shape mesh quality regularization via $\mathfrak{J}^{\tau}$ and volume mesh quality regularization $\mathfrak{J}^{\mathbb{D}}$. To avoid this issue, we formulate the volume and shape mesh regularized shape optimization problem using a bi-level approach. We have already seen in remark 2, that simultaneous shape parameterization tracking and shape optimization can be put into the bi-level framework. In fact, this was seen to be equivalent to the added regularizer approach eq. 9 with regard to gradient systems.

Let us consider weights $\alpha^{\mathbb{D}},\alpha^{\tau}>0$. Then we formulate the simultaneous volume and shape mesh regularization of shape optimization problem eq. 1 as

$\displaystyle\begin{split}\underset{\Phi\in\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})}{\min}&\;\alpha^{\mathbb{D}}\cdot\mathfrak{J}^{\mathbb{D}}(\Phi)\\ \text{ s.t. }\varphi=&\;\underset{\varphi\in\operatorname{Emb}(M,\mathbb{D})}{\operatorname{arg}\min}\;\big{(}\mathcal{J}\circ\pi\big{)}(\varphi)+\alpha^{\tau}\cdot\mathfrak{J}^{\tau}(\varphi).\end{split}$

(33)

Of course, the bi-level problem eq. 33 can be formulated for $\alpha^{\mathbb{D}}=1$ without loss of generality. To stay coherent with section 2.1 regarding pre-shape gradient systems, which feature weighted force terms, we prefer to formulate eq. 33 with a factor $\alpha^{\mathbb{D}}>0$. Notice, that the regularization of shape optimization problem eq. 1 needs to use its pre-shape extension eq. 2. This is necessary in order to rigorously apply the pre-shape regularization strategies. In contrast to eq. 9 and eq. 17, the simultaneous volume and shape mesh quality regularized problem eq. 33 is not minimizing for one pre-shape, but for two different pre-shapes $\varphi\in\operatorname{Emb}(M,\mathbb{D})$ and $\Phi\in\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$. The lower level problem solves for a pre-shape corresponding to the actual parameterized shape solving the shape mesh regularized optimization problem eq. 9. On the other hand, the upper level problem looks for a pre-shape corresponding to the parameterization of the hold-all domain $\mathbb{D}$ with specified volume mesh quality. The set of feasible solutions $\operatorname{Diff}_{\partial\mathbb{D}\cup\varphi(M)}(\mathbb{D})$ to the upper level problem depends on the lower level problem, because the optimal shape to the latter is demanded to stay invariant.

In the remainder of this section we propose a regularized gradient system for simultaneous volume- and shape mesh quality optimization during shape optimization, and prove a corresponding existence and consistency result for the fully regularized problem. As done in section 2.1, we change the space of directions from $C^{\infty}$ to $H^{1}$, as it is more suitable for numerical application. We remind the reader that the pre-shape derivative $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[V]$ is defined only for directions $V\in H^{1}_{\partial M\cup\varphi(M)}(\mathbb{D},{\mathbb{R}}^{n+1})$ which vanish on $\partial M\cup\varphi(M)$. This is inevitable, since a criterion for successful application of volume mesh regularization for shape optimization routines has to leave optimal or intermediate shapes $\varphi(M)$ invariant. If $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[V]$ was used for regularizing the gradient in style of an added source term (cf. eq. 13) for general directions $V\in H^{1}(\mathbb{D},{\mathbb{R}}^{n+1})$, it would most certainly alter shapes and interfere with shape optimization. Also, it is not possible to put Dirichlet conditions on $\varphi(M)$, or to use a restricted space of test function as in eq. 31. Doing so would prohibit shape optimization itself. Hence we have to modify $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}$, such that general directions $V\in H^{1}(\mathbb{D},{\mathbb{R}}^{n+1})$ are applicable as test functions, while the shape at hand is preserved. To resolve this problem, we introduce a projection

$$\operatorname{Pr}_{H^{1}_{\partial M\cup\varphi(M)}}:H^{1}\rightarrow H^{1}_{\partial M\cup\varphi(M)},$$

(35)

which is demanded to be the identity on $H^{1}_{\partial M\cup\varphi(M)}$. We leave the operator projecting a given direction $V\in H^{1}$ onto $H^{1}_{\partial M\cup\varphi(M)}$ general. Suitable options include the projection via solution of a least squares problem

$$\operatorname{Pr}_{H^{1}_{\partial M\cup\varphi(M)}}(V):=\underset{W\in H^{1}_{\partial M\cup\varphi(M)}}{\operatorname{arg}\min}\frac{1}{2}\|W-V\|^{2}_{H^{1}}.$$

(36)

In practice, it is feasible to construct a projection eq. 35 by using a finite element representation of $V$ and setting coefficients of basis functions on the discretization of $\varphi(M)$ to zero. With this projection operator, we can extend the volume tracking pre-shape derivative $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}$ to the space $H^{1}_{\partial M\cup\varphi(M)}$ by using

$$\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[\operatorname{Pr}_{H^{1}_{\partial M\cup\varphi(M)}}(\cdot)]:H^{1}(\mathbb{D},{\mathbb{R}}^{n+1})\rightarrow{\mathbb{R}},\;V\rightarrow\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[\operatorname{Pr}_{H^{1}_{\partial M\cup\varphi(M)}}(V)],$$

(37)

with $\mathfrak{D}\mathfrak{J}^{\mathbb{D}}$ as in eq. 28.

Now we can formulate the fully regularized pre-shape gradient system for simultaneous volume-, shape mesh quality and shape optimization. We motivate the combined gradient system with the same formal calculations as in the bi-level formulation for shape quality regularization remark 2, which are inspired by [15]. Given a symmetric, positive-definite bilinear form $a(.,.)$, the gradient system takes the form

$$a(U,V)=\mathcal{D}\mathcal{J}(\Gamma)[V]\;+\;\alpha^{\tau}\cdot\langle g_{\varphi}^{\mathcal{T}},V\rangle\;+\;\alpha^{\mathbb{D}}\cdot\mathfrak{D}\mathfrak{J}^{\mathbb{D}}(\Phi)[\operatorname{Pr}_{H^{1}_{\partial M\cup\varphi(M)}}(V)]\quad\forall V\in H^{1}(\mathbb{D},{\mathbb{R}}^{n+1}),$$

(38)

where $g_{\varphi}^{\mathcal{T}}$ is the tangential component of shape regularization eq. 8 and $\mathcal{D}\mathcal{J}(\Gamma)$ is the shape derivative of the original shape objective with $\Gamma=\pi(\varphi)$. Notice that the fully regularized pre-shape gradient system eq. 38 looks similar to the shape gradient system eq. 11 of the original problem, differing only by two added force terms on the right hand side. These force terms can be thought of as regularizing terms to the original shape gradient. In practice, this means simultaneous volume and shape mesh quality improvement for shape optimization amounts to adding two terms on the right hand side of the gradient system. Hence they can also be viewed as a (pre-)shape gradient regularization by added force terms.