# PDE-Driven Shape Optimization: Numerical Investigation of Different Descent Directions and Projections Using Penalization and Regularization

Peter Philip

► To cite this version:

Peter Philip. PDE-Driven Shape Optimization: Numerical Investigation of Different Descent Directions and Projections Using Penalization and Regularization. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.237-246, 10.1007/978-3-662-45504-3_23. hal-01286429

HAL Id: hal-01286429

https://inria.hal.science/hal-01286429

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PDE-driven shape optimization: Numerical investigation of different descent directions and projections using penalization and regularization

Peter Philip

Department of Mathematics, Ludwig-Maximilians University (LMU) Munich
Theresienstrasse 39, 80333 Munich, Germany

philip@math.lmu.de

http://www.math.lmu.de/~philip

**Abstract.** We consider shape optimization problems with elliptic partial differential state equations. Using regularization and penalization, unknown shapes are encoded via shape functions, turning the shape optimization into optimal control problems for the unknown functions. The method is designed to allow topological changes in a natural way. Based on convergence and differentiability results, numerical algorithms are formulated, using different descent directions and projections. The algorithms are assessed in a series of numerical experiments, applied to an elliptic PDE arising from an oil industry application with two unknown shapes, one giving the region where the PDE is solved, and the other determining the PDE's coefficients.

**Keywords:** shape optimization, optimal control, fixed domain method, elliptic partial differential equation, numerical simulation

## 1 Introduction

We study an elliptic shape optimization problem motivated by the oil industry application studied in [12], where one aims at monitoring the interior of a pipeline. The cross section through the pipeline is modeled by a set $D \subseteq \mathbb{R}^2$, consisting of a liquid region $\Omega$ (such that $D \setminus \Omega$ represents air) that is part oil (region $O \subseteq \Omega$) and part water ($\Omega \setminus O$). The shape optimization needs to reconstruct the parts $O$ and $\Omega$ from given measurements $y_d$ of a quantity (e.g. voltage) taken in a region $E$ adjacent to the boundary of $D$. Thus, we are led to the following problem, previously formulated in [5]:

$$ \min_{\Omega, O} \frac{1}{2} \int_E |y - y_d|^2 dx + \frac{1}{2} \int_E |\nabla y - \nabla y_d|^2 dx, \quad (1a) $$

subject to

$$ \begin{aligned} & \int_{\Omega} [a_1 \chi_O + a_2(1 - \chi_O)] \nabla y \cdot \nabla v \, dx + \int_{\Omega} [b_1 \chi_O + b_2(1 - \chi_O)] y v \, dx \\ & = \int_{\Omega} f v \, dx, \quad \forall v \in H_0^1(\Omega), \end{aligned} \quad (1b) $$
$$y - \xi \in H_0^1(\Omega), \qquad (1c)$$

where the sets $E \subseteq O \subseteq \Omega \subseteq D \subseteq \mathbb{R}^2$ all are bounded and open, $D$ also connected; $\chi_O$ denotes the characteristic function of $O$; $a_1, a_2, b_1, b_2 > 0$, $f \in L^2(D)$, $y_d, \xi \in H^1(D)$ all given. In particular, depending on $\xi$, (1c) can mean homogeneous or nonhomogeneous Dirichlet conditions.

The employed method (previously used, e.g., in [4, 5, 8]) is based on the introduction of shape functions $g$ and $p$, defined on $\bar{D}$ and encoding the unknown sets $\Omega$ and $O$, respectively, and on a technique for the approximation and regularization of characteristic functions. Assuming $g, p: \bar{D} \to \mathbb{R}$ to be continuous, the corresponding sets are obtained via

$$\Omega_g = \text{int}\{x \in D : g(x) \ge 0\}, \quad O_p = \text{int}\{x \in D : p(x) \ge 0\} \qquad (2)$$

($\Omega_g$ and $O_p$ then are open Caratheodory sets, not necessarily connected). Enforcing the constraints $E \subseteq O \subseteq \Omega$ translates into the set of admissible pairs

$$U_{\text{ad}} := \{(g, p) \in C(\bar{D}) \times C(\bar{D}) : g \ge p \text{ on } D \text{ and } p \ge 0 \text{ on } E\}. \qquad (3)$$

Let $H: \mathbb{R} \to \mathbb{R}$ denote the Heaviside function. Then $H(g), H(p): \bar{D} \to \mathbb{R}$ are the characteristic functions of $\Omega_g$ and $\bar{O}_p$, respectively. We use the differentiable regularization of the Heaviside function given by

$$H_{\varepsilon}(r) := \begin{cases} 1 & \text{for } r \ge 0, \\ \frac{\varepsilon(r+\varepsilon)^2 - 2r(r+\varepsilon)^2}{\varepsilon^3} & \text{for } -\varepsilon < r < 0, \\ 0 & \text{for } r \le -\varepsilon, \end{cases} \qquad (4)$$

and obtain the following regularized fixed domain approximation of (1) (with the abovementioned constraints $E \subseteq O \subseteq \Omega$):

$$\min_{(g,p) \in U_{\text{ad}}} \frac{1}{2} \int_E |y_\varepsilon - y_d|^2 \, dx + \frac{1}{2} \int_E |\nabla y_\varepsilon - \nabla y_d|^2 \, dx, \qquad (5a)$$

$$\int_D \left[ [a_1 H_{\varepsilon}(p) + a_2(1 - H_{\varepsilon}(p))] \nabla y_{\varepsilon} \cdot \nabla v + [b_1 H_{\varepsilon}(p) + b_2(1 - H_{\varepsilon}(p))] y_{\varepsilon} v \right] dx \\ + \frac{1}{\varepsilon} \int_D (1 - H_{\varepsilon}(g)) y_{\varepsilon} v \, dx = \int_D f v \, dx, \quad \forall v \in H_0^1(D), \qquad (5b)$$

$$y_{\varepsilon} - \xi \in H_0^1(D). \qquad (5c)$$

For $\varepsilon > 0$ small, the penalty term with the $1/\varepsilon$ in (5b) forces the state $y_\varepsilon$ to be close to 0 outside $\Omega_g$ (for precise, rigorous versions of this statement see [5, Th. 2], [8, Th. 2.2, Th. 3.1]). It is noted that, even though the above-described method encodes $\Omega_g$ and $O_p$, in fact, as level sets of the functions $g, p$, respectively, our method is essentially different from the well-known level set method of [6], since no time dependence of the functions $g, p$ and no time evolution of the
corresponding open sets $\Omega_g$, $O_p$ is assumed. In particular, we do not need to solve any Hamilton-Jacobi equations in the process.

We now consider a triangular finite element partition of $\overline{D}$, $\overline{D} = \bigcup_{T_h \in \mathcal{T}_h} T_h$, $h > 0$, assuming the grid in $D$, restricted to $E$, provides a finite element mesh in $E$ as well. Let $V_h$, $\tilde{V}_h$ denote the corresponding finite element spaces in $D$ constructed with piecewise affine continuous functions (with 0 trace on $\partial D$ for elements of $V_h$). Defining

$$ U_{\text{ad}}^{h} := \{(g, p) \in \tilde{V}_{h} \times \tilde{V}_{h} : g \ge p \text{ on } D \text{ and } p \ge 0 \text{ on } E\}, \quad (6) $$

the discretized form of (5) reads

$$ \min_{(g_h, p_h) \in U_{\text{ad}}^h} j(g_h, p_h) := \frac{1}{2} \int_E |y_{\varepsilon,h} - y_{d,h}|^2 dx + \frac{1}{2} \int_E |\nabla y_{\varepsilon,h} - \nabla y_{d,h}|^2 dx, \quad (7a) $$

$$ \int_D [a_1 H_\varepsilon(p_h) + a_2(1 - H_\varepsilon(p_h))] \nabla y_{\varepsilon,h} \cdot \nabla v_h \, dx \\ + \int_D [b_1 H_\varepsilon(p_h) + b_2(1 - H_\varepsilon(p_h))] y_{\varepsilon,h} v_h \, dx \\ + \frac{1}{\varepsilon} \int_D (1 - H_\varepsilon(g_h)) y_{\varepsilon,h} v_h \, dx = \int_D f_h v_h \, dx, \quad \forall v_h \in V_h \subseteq H_0^1(D), \quad (7b) $$

$$ y_{\varepsilon,h} - \xi_h \in V_h \subseteq H_0^1(D), \quad (7c) $$

where (7b) constitutes the equation for the discretized state $y_{\varepsilon,h} \in \tilde{V}_h \subseteq H^1(D)$, $f_h \in \tilde{V}_h$ and $\xi_h \in \tilde{V}_h$ are given suitable approximations of $f$ and $\xi$, respectively, $g_h, p_h \in \tilde{V}_h$ are discretized shape functions corresponding to discretizations of $\Omega_g$ and $O_p$, respectively, and $y_{d,h}$ is a suitable given continuous and piecewise affine approximation of $y_d$.

Similar to [5, Prop. 2] and [8, Cor. 5.3], one obtains the directional derivative of the cost functional $(g, p) \mapsto j(g, p)$ with $j$ as in (7a) at $(g_h, p_h) \in \tilde{V}_h \times \tilde{V}_h$ in the direction $(w_h, u_h) \in \tilde{V}_h \times \tilde{V}_h$ as

$$ \frac{1}{\varepsilon} \int_D H'_\varepsilon(g_h) w_h y_{\varepsilon,h} q_{\varepsilon,h} \, dx - \int_D (b_1 - b_2) H'_\varepsilon(p_h) u_h y_{\varepsilon,h} q_{\varepsilon,h} \, dx \\ - \int_D (a_1 - a_2) H'_\varepsilon(p_h) u_h \nabla y_{\varepsilon,h} \nabla q_{\varepsilon,h} \, dx, \quad (8) $$

where $q_{\varepsilon,h} \in V_h \subseteq H_0^1(D)$ is the solution to the adjoint equation

$$ \int_D [a_1 H_\varepsilon(p_h) + a_2(1 - H_\varepsilon(p_h))] \nabla q_{\varepsilon,h} \cdot \nabla v_h \, dx \\ + \int_D [b_1 H_\varepsilon(p_h) + b_2(1 - H_\varepsilon(p_h))] q_{\varepsilon,h} v_h \, dx + \frac{1}{\varepsilon} \int_D (1 - H_\varepsilon(g_h)) q_{\varepsilon,h} v_h \, dx \\ = \int_E (y_{\varepsilon,h} - y_{d,h}) v_h \, dx + \sigma \int_E (\nabla y_{\varepsilon,h} - \nabla y_{d,h}) \cdot \nabla v_h \, dx, \quad \forall v_h \in V_h, \quad (9) $$
and the direction of steepest descent ($w_{d,0}$, $u_{d,0}$) is given by

$$ 
\begin{aligned}
w_{d,0} &:= -(1/\epsilon) H'_{\epsilon}(g_h) y_{\epsilon,h} q_{\epsilon,h}, \\
u_{d,0} &:= H'_{\epsilon}(p_h) (a_1 - a_2) \nabla y_{\epsilon,h} \cdot \nabla q_{\epsilon,h} + H'_{\epsilon}(p_h) (b_1 - b_2) y_{\epsilon,h} q_{\epsilon,h}.
\end{aligned} 
\quad (10) $$

While (10) is difficult to use in practise as $H'_{\epsilon}(g_h)$ and $H'_{\epsilon}(p_h)$ are typically nonzero only in a small neighborhood of $\partial\Omega_{g_h}$ and $\partial O_{p_h}$, respectively, multiplication by nonnegative coefficients yields the following alternative descent directions ($w_{d,1}$, $u_{d,1}$) and ($w_{d,2}$, $u_{d,2}$), without such support restrictions:

$$ w_{d,1} := -y_{\epsilon,h}q_{\epsilon,h}, \quad u_{d,1} := (a_1 - a_2) \nabla y_{\epsilon,h} \cdot \nabla q_{\epsilon,h} + (b_1 - b_2) y_{\epsilon,h} q_{\epsilon,h}, \quad (11) $$

$$ w_{d,2} := w_{d,1} \chi_S, \quad u_{d,2} := u_{d,1} \chi_S, \quad (12) $$

where $\chi_S$ denotes the characteristic function of

$$ S := \{x \in D : w_{d,1}(x) \ge u_{d,1}(x)\} \cup \{x \in E : w_{d,1}(x) \ge 0 \text{ and } u_{d,1}(x) \ge 0\}. \quad (13) $$

Using (12) has the advantage of maintaining the conditions $g \ge p$ on $D$ and $p \ge 0$ on $E$.

For the numerical results presented below, we employ four variants of an algorithm of gradient with projection type making use of (approximations of) the descent directions (11) and (12). The two variants based on (11) will be called A1a and A1b, whereas the variants based on (12) will be called A2a and A2b. Moreover, variants A1a and A2a will use the admissible set $U_{ad}^h$ of (6), whereas A1b and A2b will use the modification

$$ U_{ad,b}^h := \{(g, p) \in U_{ad}^h : |\nabla g|, |\nabla p| \le 1\}, \quad (14) $$

enforcing uniformly bounded gradients for the shape functions, a condition suggested by the results of [8, Sec. 5]. Variant A1a was previously considered in [8]; the remaining three variants are new.

The four algorithms are formulated below in Sec. 2.1, with a description of their implementation in Sec. 2.2. Numerical experiments comparing the performance of the four variants are then presented in Sec. 3.

# 2 Numerical algorithms

## 2.1 Formulation

In preparation for the numerical experiments of Sec. 3, we formulate the employed algorithms. As indicated at the end of the Introduction, we use four variants of the algorithm previously published in [8], built on the earlier version of [5]. As mentioned above, we denote the four variants by A1a, A1b, A2a, A2b, where A1a is precisely the algorithm used in [8]. A1a and A1b use the descent direction (11) for line searches followed by a projection step, whereas A2a and A2b use (12), which has the advantage of remaining within $U_{ad}^h$ during the line
search, avoiding the projection. Variants A1b and A2b project into the smaller space $U_{ad,b}^h$ of (14) after each line search. The algorithms consist of the following Steps (1) – (7):

(1): Set $n := 0$ and choose initial shape functions $(g_{h,0}, p_{h,0}) \in U_{ad}^h$.

(2): Compute the solution to the state equation $y_n := \theta_{\varepsilon,h}(g_{h,n}, p_{h,n})$, where $\theta_{\varepsilon,h} : \tilde{V}_h \times \tilde{V}_h \to \tilde{V}_h$ denotes the control-to-state operator corresponding to (7b), (7c); and compute the solution to the corresponding adjoint equation $q_n := \tilde{\theta}_{\varepsilon,h}(y_n)$, where $\tilde{\theta}_{\varepsilon,h} : \tilde{V}_h \to V_h$, $y_{\varepsilon,h} \mapsto q_{\varepsilon,h}$, denotes the solution operator corresponding to (9).

(3): Compute the descent direction $(w_d^n, u_d^n)$, where $w_d^n = w_{d,1}(y_n, q_n)$ and $u_d^n = w_{d,1}(y_n, q_n)$ according to (11) for A1a and A1b, whereas $w_d^n = w_{d,2}(y_n, q_n)$ and $u_d^n = w_{d,2}(y_n, q_n)$ according to (12) for A2a and A2b.

(4): Set $\tilde{g}_{h,n} := g_{h,n} + \lambda_n w_d^n$ and $\tilde{p}_{h,n} := p_{h,n} + \lambda_n u_d^n$, where $\lambda_n \ge 0$ is determined via line search, i.e. as a solution to the minimization problem

$$ \min_{\lambda \ge 0} j(g_{h,n} + \lambda w_d^n, p_{h,n} + \lambda u_d^n). \quad (15) $$

(5): For A2a and A2b, set $(\tilde{g}_{h,n}, \tilde{p}_{h,n}) := (\tilde{g}_{h,n}, \tilde{p}_{h,n})$ (no projection is necessary to obtain $(\tilde{g}_{h,n}, \tilde{p}_{h,n}) \in U_{ad}^h$); for A1a and A1b, set $(\tilde{g}_{h,n}, \tilde{p}_{h,n}) := \pi_h(\tilde{g}_{h,n}, \tilde{p}_{h,n})$, where $\pi_h$ denotes the projection $\pi_h : \tilde{V}_h \times \tilde{V}_h \to U_{ad}^h$, obtained by first setting $\tilde{g}_{h,n}(x_i^h) := \max\{0, \tilde{g}_{h,n}(x_i^h)\}$ and $\tilde{p}_{h,n}(x_i^h) := \max\{0, \tilde{p}_{h,n}(x_i^h)\}$ for each node $x_i^h$ of the triangulation $\mathcal{T}_h$ such that $x_i^h \in \bar{E}$, and second setting $\tilde{p}_{h,n}(x_i^h) := \min\{\tilde{p}_{h,n}(x_i^h), \tilde{g}_{h,n}(x_i^h)\}$ for every node $x_i^h$ of the triangulation $\mathcal{T}_h$.

(6): For A1a and A2a, set $(g_{h,n+1}, p_{h,n+1}) := (\tilde{g}_{h,n}, \tilde{p}_{h,n})$ (no second projection necessary); for A1b and A2b, set $(g_{h,n+1}, p_{h,n+1}) := \pi_{h,b}(\tilde{g}_{h,n}, \tilde{p}_{h,n})$, where $\pi_{h,b}$ denotes the projection $\pi_{h,b} : U_{ad}^h \to U_{ad,b}^h$, obtained by dividing $\tilde{g}_{h,n}$ and $\tilde{p}_{h,n}$ by $\alpha$, defined as the max of the max-norms of $|\nabla g|$ and $|\nabla p|$, in case $\alpha > 1$.

(7): RETURN $(g_{h,fin}, p_{h,fin}) := (g_{h,n+1}, p_{h,n+1})$ if the change of $g, p$ and/or the change of $j(g, p)$ are below some prescribed tolerance parameter. Otherwise: Increment $n$, i.e. $n := n + 1$ and GO TO (2).

For all the numerical examples discussed below, we stopped the iteration and returned $(g_{h,fin}, p_{h,fin}) := (g_{h,n+1}, p_{h,n+1})$ if $|j(g_{h,n}, p_{h,n}) - j(g_{h,n+1}, p_{h,n+1})| < 10^{-5}$ AND $\|g_{h,n} - g_{h,n+1}\|_2 < 10^{-3}$ AND $\|p_{h,n} - p_{h,n+1}\|_2 < 10^{-3}$, where $|j(g_{h,n}, p_{h,n}) - j(g_{h,n+1}, p_{h,n+1})|/|j(g_{h,n+1}, p_{h,n+1})|$ is used for $|j(g_{h,n}, p_{h,n}) - j(g_{h,n+1}, p_{h,n+1})|/|j(g_{h,n+1}, p_{h,n+1})| > 1$ and analogous for $g_{h,n}$ and $p_{h,n}$.

## 2.2 Implementation

The state equations as well as the adjoint equations that need to be solved numerically during the above algorithms are discretized linear elliptic PDE with Dirichlet boundary conditions. The numerical solution is obtained via a finite volume scheme [7, Sec. 4]. More precisely, the software WIAS-HITNIHS<fnref n="1" />, originally designed for the solution of more general PDE occurring when modeling conductive-radiative heat transfer and electromagnetic heating [2], has been

<footnote n="1">¹ High Temperature Numerical Induction Heating Simulator.</footnote>
adapted for use in the present context. WIAS-HITNIHS is based on the program package *pdelib* [1], it employs the grid generator *Triangle* [11] to produce constrained Delaunay triangulations of the domains, and it uses the sparse matrix solver *GSPAR* [3] to solve the linear system arising from the finite volume scheme.

The numerical scheme yields discrete $y_n$ and $q_n$ (cf. Step (2) of the above algorithms), defined at each vertex of the triangular discrete grid, interpolated piecewise affine, i.e. affinely to each triangle of the discrete grid. In consequence, the shape functions $g_{h,n}$ and $p_{h,n}$ are piecewise affine as well. Where integrals of these piecewise affine functions need to be computed (e.g. in Step (7) of the algorithms), they are computed exactly. A golden section search [10, Sect. 10.2] is used to numerically carry out the minimization (15). Note that the minimization (15) is typically nonconvex and the golden section search will, in general, only provide a local min $\lambda_n$.

For some numerical examples, the stated initial shape functions ($g_{h,0}$, $p_{h,0}$) are merely piecewise continuous (cf. the Introduction and [9]) and, thus, not in $U_{\text{ad}}^h$. However, the stated ($g_{h,0}$, $p_{h,0}$) are only used to determine the values $g_h(x_i^h)$, $p_h(x_i^h)$, at the nodes $x_i^h$ of the triangulation $\mathcal{T}_h$, and the resulting affinely interpolated functions are in $U_{\text{ad}}^h$. Moreover, in Step (3) of the algorithms, approximations of the descent directions are used, as for the gradients nodewise averages are computed, that are then affinely interpolated, and the conditions of (13) are enforced nodewise and affinely interpolated. In principle, in might occur that the approximated direction is no longer a descent direction, but such a case was not observed during our numerical experiments.

# 3 Numerical experiments

## 3.1 Numerical experiments with precomputed optimum

The numerical computations of the present section employ the circular fixed domain $D := \{(x_1, x_2) : x_1^2 + x_2^2 < 1\} \subseteq \mathbb{R}^2$ with fixed subdomain $E := \{(x_1, x_2) \in D : |x_1| > \frac{3}{4}, |x_2| < \frac{1}{2}\} \subseteq \bar{D}$ (note $E$ has two connected components). We use a fixed triangular grid provided by *Triangle* [11], consisting of 24458 triangles. The used regularization parameter is $\varepsilon = 10^{-5}$ (cf. [8, 9]). The settings for the remaining given quantities are $a_1 := 1$, $a_2 := 10$, $b_1 := 1$, $b_2 := 10$, $f(x_1, x_2) := 5$, $\xi(x_1, x_2) := 2$. The cost functional $j$ as in (7a) depends on the given function $y_{d,h}$. For the first set of numerical results, we precompute $y_{d,h} := y_{\varepsilon,h}$ numerically as the solution to the state equation (7b), (7c), using

$$ g_h(x_1, x_2) := \begin{cases} -1 & \text{if } (x_1, x_2) \notin E \text{ and } \|(x_1, x_2) - (-1, 0)\|_2 < 0.4, \\ -1 & \text{if } (x_1, x_2) \notin E \text{ and } \|(x_1, x_2) - (1, 0)\|_2 < 0.4, \\ 1 & \text{otherwise,} \end{cases} \quad (16a) $$

$$ p_h(x_1, x_2) := \begin{cases} 1 & \text{in } E, \\ -1 & \text{in } D \setminus E. \end{cases} \quad (16b) $$
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Precomputed $y_{d,h}$ used in all experiments of Sec. 3.1 (left, isolevels spaced at 0.2), obtained as the solution to the state equation (7b), (7c); with the corresponding $\Omega_g$ (middle) and $O_p$ (right).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Initial state, shapes used in all experiments of Sec. 3.1. Left: State isolevels spaced at 0.2. Middle: Shape $\Omega_g$. Right: Shape $O_p$. Cost: $j(g_{h,0}, p_{h,0}) = 19.0$.</figcaption></figure>

The computed $y_{d,h}$ with the corresponding $\Omega_g$ and $O_p$ is depicted in Fig. 1. Using the precomputed $y_{d,h}$ has the advantage that we actually know $y_{d,h}$ together with $g_h, p_h$ as in (16) provides an absolute minimum in the following numerical examples, employing the cost functional $j$ of (7a) with the precomputed $y_{d,h}$ from above. A series of four numerical experiments was conducted, all using the initial shape functions $g_{h,0}(x_1, x_2) := 1$, $p_{h,0}(x_1, x_2) := 1$ (see Fig. 2).

We refer to the experiments as 1:A1a, 1:A1b, 1:A2a, and 1:A2b, depending on which variant of the algorithm of Sec. 2.1 was used. The results for 1:A2a and 1:A2b are shown in Fig. 3. The final state and shapes for 1:A1a and 1:A1b were very similar to those of 1:A2a, with slightly higher final costs (0.69 and 0.28, respectively). All variants reduce the cost significantly, all resulting local minima being different and different from the absolute min. Variant A2a gives the best result, whereas A2b results in the highest final cost, where one also observes a symmetry breaking due to the discrete grid. Actually, for A2b, after the first line search, the cost is 0.29 with shapes resembling the final shapes of the other variants, but the projection of Step (6) can subsequently result in a cost increase, which occurs in this example.

## 3.2 Numerical experiments without precomputed optimum

In contrast to the experiments of the previous section, we now consider a setting, where we are no longer in the situation of a known precomputed optimum. For the following numerical results, the fixed domain $D$ is still the unit disk as in Sec. 3.1. However, the fixed subdomain $E$ is now at the bottom of $D$, defined by $E := \{(x_1, x_2) \in D : x_2 < -0.7\} \subseteq D$. The numerical computations employ a
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Final state, shapes for shape optimizations 1:A2a (1st row) and 1:A2b (2nd row) of Sec. 3.1. Left: State isolevels spaced at 0.2. Middle: Shapes $\Omega_g$. Right: Shapes $O_p$. Final costs $j(g_{h,fin}, p_{h,fin})$ are 0.053 for 1:A2a, 1.30 for 1:A2b. Required number of line searches: 6 for 1:A2a, 29 for 1:A2b.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Initial state, shapes used in all experiments of Sec. 3.2. Left: State isolevels spaced at 0.2. Middle: Shape $\Omega_g$. Right: Shape $O_p$. Cost: $j(g_{h,0}, p_{h,0}) = 24.8$.</figcaption></figure>

fixed triangular grid provided by Triangle [11], consisting of 24623 triangles. The parameter settings are as in Sec. 3.1, except for $f(x_1, x_2) := 10(x_1^2+x_2^2)+5$. The cost functional is as in (7a) with $y_{d,h}(x_1, x_2) := x_1+x_2$. A series of four numerical experiments was conducted, all using the initial shape functions $g_{h,0}(x_1, x_2) =$

$$p_{h,0}(x_1, x_2) := \begin{cases} 1 & \text{if } (x_1, x_2) \in E, \\ -1 & \text{otherwise} \end{cases} \quad \text{(see Fig. 4).}$$

We refer to the experiments as 2:A1a, 2:A1b, 2:A2a, and 2:A2b, depending on which variant of the algorithm of Sec. 2.1 was used. Results are shown in Fig. 5, except for 2:A2a, which converged after 10 line searches to a local min almost identical to the initial condition. All other variants reduce the cost significantly, all resulting local minima being different. Here, the lowest final cost is achieved for variant A1b. One notices significant changes in shapes (including topology changes) during the optimizations, where very different shapes can result in nearly identical costs. As in Sec. 3.1, symmetry breaking can occur due to the discrete grid.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Intermediate and final state, shapes for shape optimizations of Sec. 3.2, i.e. for Experiments 2:A1a (1st, 2nd row), 2:A1b (3rd, 4th row), and 2:A2b (5th, 6th row). Left: State isolevels spaced at 0.2. Middle: Shapes $\Omega_g$. Right: Shapes $O_p$. Costs at shown intermediate states are 2.64 for 2:A1a, 1.74 for 2:A1b, 1.70 for 2:A2b. Final costs $j(g_{h,fin}, p_{h,fin})$ are 1.72 for 2:A1a, 1.68 for 2:A1b, 1.69 for 2:A2b. Number of line searches for intermediate and for final state: 8 and 31 for 2:A1a, 2 and 29 for 2:A1b, 3 and 38 for 2:A2b.</figcaption></figure>
Conclusions

In a series of numerical experiments, we have studied four variants of an algo-
rithm of gradient with projection type for shape optimization problems driven
by elliptic PDE. The variants used different descent directions and different sets
of admissible shape functions. Except in one situation, all variants were effective
in finding local minima of significantly reduced costs. However, it did depend on
both the equation and on the initial condition, which variant showed the best
performance. Thus, further research seems warranted to further evaluate and
improve the different variants.

References

1. Fuhrmann, J., Koprucki, T., Langmach, H.: pdelib: An open modular tool box for the numerical solution of partial differential equations. Design patterns. In: Proceedings of the 14th GAMM Seminar on Concepts of Numerical Software, Kiel, January 23 -25, 1998. University of Kiel, Kiel, Germany (2001)

2. Geiser, J., Klein, O., Philip, P.: Numerical simulation of temperature fields during the sublimation growth of SiC single crystals, using WIAS-HiTNIHS. J. Crystal Growth 303, 352–356 (2007)

3. Grund, F.: Direct linear solvers for vector and parallel computers. In: Vector and Parallel Processing VECPAR’98, Lecture Notes in Computer Science, vol. 1573, pp. 114–127. Springer, Berlin (1999)

4. Mäkinen, R., Neittaanmäki, P., Tiba, D.: On a fixed domain approach for a shape optimization problem. In: Ames, W., van Houwen, P. (eds.) Computational and Applied Mathematics II: Differential Equations, pp. 317–326. North Holland, Amsterdam (1992)

5. Neittaanmäki, P., Pennanen, A., Tiba, D.: Fixed domain approaches in shape op-
timization problems with Dirichlet boundary conditions. Inverse Problems 25(5),
1–18 (2009)

6. Osher, S., Sethian, J.: Fronts propagating with curvature dependent speed: Al-
gorithms based on Hamilton-Jacobi formulations. J. Comput. Phys. 79(1), 12–49
(1988)

7. Philip, P.: Analysis, optimal control, and simulation of conductive-radiative heat transfer. Mathematics and its Applications / Annals of AOSR 2, 171–204 (2010)

8. Philip, P., Tiba, D.: A penalization and regularization technique in shape optimization problems. SIAM J. Control Optim. 51(6), 4295–4317 (2013)

9. Philip, P., Tiba, D.: Shape optimization via control of a shape function on a fixed domain: Theory and numerical results. In: Numerical Methods for Differential Equations, Optimization, and Technological Problems, Computational Methods in Applied Sciences, vol. 27, pp. 305–320. Springer, New York (2013)

10. Press, W., Teukolsky, S., Vetterling, W., Flannery, B.: Numerical Recipes. The Art of Scientific Computing. Cambridge University Press, New York, USA, 3rd edn. (2007)

11. Shewchuk, J.: Delaunay refinement algorithms for triangular mesh generation.
Comput. Geom. 22(1-3), 21–74 (2002)

12. Woo, H., Kim, S., Seol, J., Lionheart, W., Woo, E.: A direct tracking method for a grounded conductor inside a pipeline from capacitance measurements. Inverse Problems 22, 481–494 (2006)