---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.19.1
  main_language: python
---

```{eval-rst}
.. include:: ../../../global.rst
```

(demo_deflation)=
# Computing Multiple Local Minimizers of Topology Optimization Problems - Five-holes Double-pipe

## Problem Formulation

Usually, topology optimization problems attain multiple local minimizers. In this demo,
we compute multiple local minimizers of the five-holes double-pipe optimization problem.
This problem has been investigated previously, e.g., in
[Papadopoulos, Farrell, Surowiec - Computing multiple solutions of topology optimization
problems](https://doi.org/10.1137/20M1326209) and [Baeck, Blauth, Leithaeuser, Pinnau,
Sturm - A Novel Deflation Approach for Topology Optimization and Application for Optimization
of Bipolar Plates of Electrolysis Cells](https://doi.org/10.1137/24M1670913). The problem can be written
as follows

$$
\begin{align}
&\min_{\Omega,u} J(\Omega, u) = \int_\mathrm{D} \nabla u : \nabla u +
\alpha_\Omega u \cdot u \text{d}x \\
&\text{subject to} \qquad
\begin{alignedat}{2}
    -\Delta u + \nabla p + \alpha_\Omega u &= 0 \quad &&\text{ in } \mathrm{D},\\
    \nabla \cdot u &= 0 \quad &&\text{ in } \mathrm{D},\\
    u &= u_D \quad &&\text{ on } \partial \mathrm{D},\\
    V_L &\leq |\Omega| \leq V_U.
\end{alignedat}
\end{align}
$$

+++

Here, {math}`u` and {math}`p` denote the fluids velocity and pressure, respectively.
Moreover,
{math}`\alpha` denotes the viscous resistance or inverse permeability of the material,
which is used to distinguish between fluid (where {math}`\alpha` is low) and solid
(where {math}`\alpha` is high). As in the previous demos, e.g. {ref}`demo_pipe_bend`,
{math}`\alpha` represents a jumping coefficient between the considered materials,
i.e., it is given by {math}`\alpha_\Omega(x) =
\chi_\Omega(x)\alpha_\mathrm{in} + \chi_{\Omega^c}(x) \alpha_\mathrm{out}`, so that
{math}`\Omega` models our fluid and {math}`\Omega^c` models the solid part.

On the outer boundary of the hold-all domain,
Dirichlet boundary conditions are specified,
indicating, where the fluid enters and exits. Moreover, the goal of the optimization
problem is to minimize the energy dissipation of the fluid while achieving a certain
volume of the fluid region. For more details on this problem, we refer the reader
to [Baeck, Blauth, Leithaeuser, Pinnau, Sturm - A Novel Deflation Approach for
Topology Optimization and Application for Optimization of Bipolar Plates of
Electrolysis Cells](https://doi.org/10.1137/24M1670913).

The generalized topological derivative for this problem is given by

$$
\mathcal{D}J(\Omega)(x) = \left(\alpha_\mathrm{in} - \alpha_\mathrm{out}\right)
u(x)\cdot (u(x) + v(x)).
$$

:::{note}
We do not specify the adjoint equation here, as this is derived automatically
by cashocs as in previous demos.
:::

## Implementation

The complete python code can be found in the file {download}`demo_deflation.py
</../../demos/documented/topology_optimization/deflation/demo_deflation.py>`,
and the corresponding config can be found in {download}`config.ini
</../../demos/documented/topology_optimization/deflation/config.ini>`.

The first part of the code is completely analogous to {ref}`demo_pipe_bend`
and therefore we omit the details here and just state the corresponding code

```python
from fenics import *

import cashocs

cfg = cashocs.load_config("config.ini")
mesh, subdomains, boundaries, dx, ds, dS = cashocs.import_mesh("mesh_deflation.xdmf")

v_elem = VectorElement("CG", mesh.ufl_cell(), 2)
p_elem = FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_elem * p_elem)
CG1 = FunctionSpace(mesh, "CG", 1)
DG0 = FunctionSpace(mesh, "DG", 0)

alpha_in = 2.5 / (100 ** 2)
alpha_out = 2.5 / (0.01 ** 2)
alpha = Function(DG0)

psi = Function(CG1)
psi.vector().vec().set(-1.0)
psi.vector().apply("")

up = Function(V)
u, p = split(up)
vq = Function(V)
v, q = split(vq)

F = (
    inner(grad(u), grad(v)) * dx
    - p * div(v) * dx
    - q * div(u) * dx
    + alpha * dot(u, v) * dx
)
```

For the Dirichlet boundary conditions, we specify that we have two inflows at
left part of the boundary as well as two outflows on the right part,
with the following expressions. Altogether, the boundary conditions are specified
using the code

```python
v_1 = Expression(('-144*(x[1]-1.0/6)*(x[1]-2.0/6)', '0.0'), degree=2)
v_2 = Expression(('-144*(x[1]-4.0/6)*(x[1]-5.0/6)', '0.0'), degree=2)
bcs = cashocs.create_dirichlet_bcs(V.sub(0), Constant((0, 0)), boundaries, [1])
bcs += cashocs.create_dirichlet_bcs(V.sub(0), v_1, boundaries, [3, 5])
bcs += cashocs.create_dirichlet_bcs(V.sub(0), v_2, boundaries, [2, 4])
```

Additionally, we specify the volume constraint by

```python
vol_des = 0.5
```

We define the cost functional of our problem as well as its corresponding
generalized topological derivative with the lines, where the penalization
for the volume constraint is not considered here in comparison to
{ref}`demo_pipe_bend`

```python
J = cashocs.IntegralFunctional(
    inner(grad(u), grad(u)) * dx
    + alpha * dot(u, u) * dx
)
dJ_in = Constant(alpha_in - alpha_out) * (dot(u, v) + dot(u, u))
dJ_out = Constant(alpha_in - alpha_out) * (dot(u, v) + dot(u, u))
```

::::{note}
As in {ref}`demo_poisson_clover`, the generalized topological derivative for this
problem is identical in {math}`\Omega` and {math}`\Omega^c`, which is usually not the
case. For this reason, the special structure of the problem can be exploited with the
following lines in the configuration file
```{code-block} ini
:caption: config.ini
[TopologyOptimization]
topological_derivative_is_identical = True
```
::::

+++

As in the previous demos, we have to specify the update routine of the level-set
function `update_level_set`, which we do as follows:

```python
def update_level_set():
    cashocs.interpolate_levelset_function_to_cells(psi, alpha_in, alpha_out, alpha)
```

That is, in the `update_level_set` function, the jumping coefficient is updated
with the {py:func}`interpolate_levelset_function_to_cells
<cashocs.interpolate_levelset_function_to_cells>` function.

### The Deflation Approach

To compute multiple local minizers of this topology optimization problem, we use
the approach presented in [Baeck, Blauth, Leithaeuser, Pinnau,
Sturm - A Novel Deflation Approach for Topology Optimization and Application for Optimization
of Bipolar Plates of Electrolysis Cells](https://doi.org/10.1137/24M1670913). Here,
the distance to previously found local minimizers is penalized in the objective
function. The deflated topology optimization problem
{py:class}`DeflatedTopologyOptimizationProblem <cashocs.DeflatedTopologyOptimizationProblem>`
can then be defined via the lines

```python
dtop = cashocs.DeflatedTopologyOptimizationProblem(
    F, bcs, J, up, vq, psi, dJ_in, dJ_out, update_level_set, config=cfg,
    volume_restriction=vol_des
)
```

::::{note}
The definition of the deflated topology optimization problem is similar to the
definition of a usual topology optimization problem, see, e.g.,
{ref}`demo_shape_stokes`.
::::

+++

::::{note}
The volume constraint is handled as in {ref}`demo_projection` by a projection of the
level-set function.
::::

Now, we can solve the topology optimization problem with the following lines

```python
delta = 1000000.
gamma = 0.7
dtop.solve(tol=1e-6, it_deflation=4, gamma=gamma, delta=delta, inner_rtol=0., inner_atol=0., angle_tol=1.0)
```

Here, `it_deflation` denotes the number of iterations of the deflation procedure, and
`gamma` and `delta` are the penalty parameters of the penalty function. The Parameters
`gamma` denotes the local support of the penalty function, and the parameter `delta`
is the penalty parameter. Overall, the code results in a maximum of two times
`it_deflation` solves of a topology optimization problem.

:::{note}
The choice of the penalty parameters `gamma` and `delta` is not trivial and was chosen by
a numerical parameter studies. For a more detailed discussion about the choices of the
parameters, we refer to [Baeck, Blauth, Leithaeuser, Pinnau,
Sturm - A Novel Deflation Approach for Topology Optimization and Application for Optimization
of Bipolar Plates of Electrolysis Cells](https://doi.org/10.1137/24M1670913).
:::

:::{note}
In density based topology optimization, the optimization problems reduces to an optimal
control problem. A deflated optimal control problem can be defined in a similar fashion
with
{py:class}`DeflatedOptimalControlProblem <cashocs.DeflatedOptimalControlProblem>`.
The solution of the optimization problem is then analog to the topology optimization
case.
:::

### Visualization

We visualize the local minimizers found by the deflation procedure with

```python
from matplotlib import colors
```

```python
import matplotlib.pyplot as plt
import numpy as np

rgbvals = np.array([[0, 107, 164], [255, 128, 14]]) / 255.0
cmap = colors.LinearSegmentedColormap.from_list(
    "tab10_colorblind", rgbvals, N=256
)

for i in range(0, len(dtop.control_list_final)):
    plot(dtop.control_list_final[i], cmap=cmap, vmin=-1e-10, vmax=1e-10)
    plt.tight_layout()
    plt.savefig("./deflation_shape_{i}.png".format(i=i), bbox_inches='tight')
    # plt.show()
```

and the results look as follows
![](/../../demos/documented/topology_optimization/deflation/deflation_shape_0.png)
![](/../../demos/documented/topology_optimization/deflation/deflation_shape_1.png)
![](/../../demos/documented/topology_optimization/deflation/deflation_shape_2.png)
![](/../../demos/documented/topology_optimization/deflation/deflation_shape_3.png)
![](/../../demos/documented/topology_optimization/deflation/deflation_shape_4.png)