# Neural Level Set Topology Optimization Using unfitted Finite Elements

Connor N. Mallon1∗ ,  Aaron W. Thornton2 ,  Matthew R. Hill1,2 and  Santiago
Badia3∗

###### Abstract.

To facilitate the widespread adoption of automated engineering design
techniques, existing methods must become more efficient and generalizable. In
the field of topology optimization, this requires the coupling of modern
optimization methods with solvers capable of handling arbitrary problems. In
this work, a topology optimization method for general multiphysics problems is
presented. We leverage a convolutional neural parameterization of a level set
for a description of the geometry and use this in an unfitted finite element
method that is differentiable with respect to the level set everywhere in the
domain. We construct the parameter to objective map in such a way that the
gradient can be computed entirely by automatic differentiation at roughly the
cost of an objective function evaluation. Without handcrafted initializations,
the method produces regular topologies close to the optimal solution for
standard benchmark problems whilst maintaining the ability to solve a more
general class of problems than standard methods, e.g., interface-coupled
multiphysics.

1 Department of Chemical and Biological Engineering, Monash University,
Wellington Rd Clayton, 3800, Victoria, Australia. 2 CSIRO, Research Way
Clayton, 3168, Victoria, Australia, 3 School of Mathematics, Monash
University, Wellington Rd Clayton, 3800, Victoria, Australia. ∗ Corresponding
authors E-mails<EMAIL_ADDRESS>(Connor Mallon, Department of
Chemical and Biological Engineering, Monash University, Wellington Rd Clayton,
3800, Victoria, Australia<EMAIL_ADDRESS>(Santiago Badia, School
of Mathematics, Monash University, Wellington Rd Clayton, 3800, Victoria,
Australia)

##  1\. Introduction

After the birth of topology optimizations in the field of structural design
[Bendse1989], efforts have been made to increase the effectiveness of such
automated design approaches and allow for their deployment on a more general
class of problems [Sigmund2013, Guo2010].

A plethora of TO strategies exist through the literature, the most common of
which being density-based methods using the so-called SIMP (Solid isotropic
microstructure with penalization for intermediate densities) method
[Sigmund2001]. These involve varying a material distribution continuously
between 0 and 1 to introduce an artificial representation of the boundary.
Although simple for basic structural problems, a way to represent intermediate
design variables arising at the boundary must be included, which becomes
increasingly complex in multiphysics applications and makes imposing arbitrary
boundary conditions non-trivial [Yoon2014].

An alternative technique that can overcome some of the problems presented by
density methods and tackle a more general class of problems (e.g., interface-
coupling multiphysics and problems that involve surface PDEs on boundaries) is
the level set (LS) TO method [Osher1988, SethianJamesAlbert1999Lsma]. Using
this approach, the boundary is described by the zero iso-surface of an LS
function. It is instead this LS function that is varied to obtain optimized
designs. A precise location of the boundary is then available.

A variety of alternative implementations of the LS TO method have been made
[vanDijk2013]. They can be distinguished, among other things, by how they
update the topology at each iteration and their means of geometry mapping. The
methods to update geometries involve either updating the solution of Hamilton-
Jacobi equations by a velocity field based on sensitivity information
[Osher2001, Burman2018] or using a parameterization of the topology that is an
explicit function of the design variables of a steepest descent optimization
scheme. The latter approach allows one to leverage well-established nonlinear
programming techniques and is the method selected for this work.

Types of geometry mappings include using the LS function to define a conformal
mesh to the boundary (see e.g. [Ha2008, Yamasaki2011]) which requires re-
meshing at each iteration, density-based mappings (see e.g. [Allaire2004,
Wang2003, Dugast2020]), which recover some of the issues related to density
methods, or unfitted/immersed boundary techniques (see, e.g., [Parvizian2011,
Burman2015, badia_stokes_2018]). Unfitted methods rely on a fixed background
mesh and capture the precise location of the boundary using triangulations of
the cells cut by the LS function. By doing so, re-meshing is avoided yet an
accurate description of the interface is maintained. In contrast to density-
based mappings, or ersatz material approaches, boundary conditions other than
Neumann conditions can easily be imposed on the precise location of the
interface. This is critical for the handling of multiphysics problems where
one must impose transmission conditions. Furthermore, the sharp treatment of
the boundary in the method also means that there is no integration error on
any interface or boundary.

A known issue with unfitted techniques is the ill-conditioning problem
associated with small cut elements. The common XFEM [Kreissl2012,
Villanueva2017] approach uses a finite element (FE) space restricted to the
interior domain and cut cells for the solution and requires stabilization in
the vicinity of the boundary by, for example, ghost penalty terms [Burman2010]
or cell aggregation [badia_aggregated_2017, Badia2022-linking]. These methods
are consistent and can provide high-order approximation [Badia2022-high]
however the support of the stabilization terms changes depending on the
location of the cut cells, leading to potential non-differentiability in the
optimization problem which can harm the convergence of gradient-based
optimization algorithms. The specific unfitted TO technique used in this work
is instead a version of the finite cell method (FCM) [Parvizian2011], in which
a non-consistent penalty term is added everywhere in the fictitious domain
(outside the physical domain) to provide robustness. This stabilization is
suitable for TO because it is differentiable with respect to the level set
parameterization (see Section 5). An implementation of the FCM for TO is made
in [Parvizian2011], which uses a refined grid for the material boundary
compared to the solution to capture fine-scale geometry. We instead use
subgrid triangulations using the LS function as in [Kreissl2012] to capture
fine-scale structure in the integration and thus avoid the need to increase
the number of design variables parameterizing the geometry. The loss of
consistency of the FCM is not an issue in TO, where high-order approximations
are not very relevant.

When utilizing nonlinear programming techniques, the user is free to select a
particular parameterization of the geometry. With a mesh already defined for
the FE problem, it is natural to also use a FE function for a discrete
representation of the LS. Doing so, a parameterization is obtained with a
user-controlled resolution. A common approach to the optimization problem is
then to take the degrees of freedom (DOF) values of this FE function as the
design variables [Kreissl2011, Kreissl2012, Dijk2012]. This choice, however,
means that each parameter is only capable of a local influence on the
geometry. This can result in the optimizer making improvements locally without
working to find the most performant overall structure. Another option for
parameterization is to use B-splines as the basis functions for the level set
[Nol2020, Wang2019, Nol2022]. These approaches increase the level of
smoothness of the level set without having to use high order FEs. However,
although the increased smoothness eliminates some of the need for filtering
techniques, the design parameters are still only capable of influencing a
local region of the domain. Methods that allow for a larger influence of
design parameters on the LS function have also been proposed [Pingen2009,
Liu2020, Wang2006]. These methods use large supports for the basis functions
by utilizing, for example, radial basis functions to allow for wider influence
of a decreased number of design parameters. Using only parameters which have a
widespread influence can however result in an inability to describe small
spatial variations of the interface.

To incorporate design parameters that simultaneously optimize the geometry at
multiple scales and retain the advantages of both local and widespread design
parameter influence, we propose to parameterize the LS using a modified U-net
convolutional neural network. Here, the design variables have an influence at
a scale that depends on the layer that they appear within the network. Using
this method, we can discover geometries that have greater regularity than
those using only local design parameters by making improvements to large-scale
features during evolution. We also, however, maintain the ability to describe
small spatial variations of the interface. The nonlinearity of the U-net and
the nature of the connectivity of the parameters allow for the discovery of
complex relationships between features at different scales which can
ultimately lead to the emergence of high-performing regular geometries and
avoid the sub-optimal solutions of common parameterizations [Yan2018]. It is
noted that some regularity can be added to the previously discussed methods,
but it often requires perturbation of the objective function, see e.g.
[Dijk2012].

The combination of machine learning and TO was explored as early as the 1990s
[Adeli1995] but has gained massive momentum in recent years [Zhang2021,
Woldseth2022]. Neural networks and other ML techniques can be incorporated
into the TO process in many ways. Common data-driven approaches attempt to
train networks to map problem descriptions directly to a geometry [Hoang2022,
Yu2018, Li2019, Zheng2021]. These however require pre-training on already
optimized samples and suffer from a lack of generalisability [Woldseth2022].
Others replace some or all of the optimization loop for accelerated
convergence by training an auxiliary network [Kallioras2020, Joo2021]. These
approaches are based on the premise, which in general is not necessarily true,
that early iterations of the optimization contain the information to produce
performant optimal geometries. An alternative method is the inclusion of a NN
for a parameterization of the geometry [Deng2020, Chandrasekhar2020,
Hoyer2019]. These approaches typically optimize the parameters of a NN
representing a continuous function that maps positions in space to a density.
These approaches tend to focus on reducing the dimensionality of the design
space assuming that NNs can efficiently achieve expressiveness with a small
number of parameters [Barron1994]. A reduction in parameters does not,
however, necessarily lead to faster convergence for NNs [Chandrasekhar2020]
compared to the standard SIMP approach.

Instead of focusing on a neural parameterization that reduces the
dimensionality of the problem, we select a network description of the geometry
which is specifically designed to learn effectively on problems involving the
segmentation of a domain. The network used in this case is a modification of
the U-Net convolutional network. The U-Net architecture was originally
developed for biomedical image segmentation tasks [ronneberger2015unet] but
has proven successful for a variety of applications in which multi-scale
features and spatial correlation is important [Ulyanov2020]. These networks
are typically composed of encoding and decoding halves. The encoding section
maps the context of input images into a low dimensional latent space which is
localized in the upsampling section to provide segmentation at the desired
resolution. The work in [Hoyer2019] exploited the properties of this network
showing improved performance with a U-Net density parameterization for a SIMP
structural optimization problem. Similar to [Hoyer2019], we use a trainable
input vector for the network and feed this into the up-sampling half of the
U-Net. In contrast to most applications of this network, we have no input
image and therefore do not need the encoding of half of the network. It is the
up-sampling (or decoding) part of the network that provides the
parameterization of the multi-scale features which are important in this
context.

When using LS TO methods, the optimized designs can be dependent on the
initialization of the design parameters. This is because sensitivity
information exists only on the design boundary making it difficult to discover
where the introduction of a new hole away from the boundary in the design
domain would be beneficial. Depending on the application, different approaches
to initializing the domain can be made. For small-scale feature design, one
may choose to seed the geometry with irregular structure as in [Dapogny2019]
but for most level set methods, evenly spaced holes are usually the starting
point [Jenkins2016, Villanueva2017, Kreissl2012]. The final design in these
applications is then dependent on the size, shape and location of the holes
[Barrera2020]. In our approach, we alleviate some of the dependency of
initialization by taking many random seeds for the starting geometry.

Using unfitted/immersed techniques is currently the leading way to retain
accurate descriptions of the interface in topology optimization
[Alexandersen2020]. These methods are typically computationally expensive,
however, because finite difference schemes are required for the computation of
terms in the backward pass [Sharma2016]. To the best of our knowledge, we are
the first to implement a fully automatically differentiable framework for an
embedded method in topology optimization so that the backward pass takes
roughly the same time as solving the forward problem.

The Julia [Julia-2017] programming language is used to implement all aspects
of this project with the FE toolbox Gridap [Badia2020, Verdugo2022] being the
main package utilized. We also use the Julia machine learning library Flux
[Flux.jl-2018] for the implementation of the U-Net. Using these foundations,
we implement a routine combining NNs with an unfitted FE based TO. The main
contributions of this work are the presentation of:

  * •

An unfitted LS TO method with a NN parameterization that avoids sub-optimal
solutions and achieves regular optimized geometries without handcrafted
initializations and

  * •

a fully automatically differentiable unfitted LS TO method for multiphysics
problems with complex boundary conditions.

We present the overall framework as follows. First, in Section 2, we present
the entire optimization loop at a high level. We then go into more detail
about various stages in the loop. Details of the architecture of the neural
network are found in Section 3, details of the geometry processing are
presented in Section 4, the numerical discretization of the problem is
presented in Section 5 and the gradient implementation in Section 6. We then
benchmark the method against baseline methods and show the generality of the
method with an application to a multiphysics problem with complex boundary
conditions in Section 7.

##  2\. optimization Problem

In this section, we provide a succinct overview of the overall TO algorithm
proposed in this work. We aim to solve the problem:

$\displaystyle\underset{\mathbf{p}}{\text{min}}$ $\displaystyle\
J(\bm{u}(\mathbf{p}),\mathbf{p})$ (1) s.t $\displaystyle\
\mathscr{R}(\bm{u}(\mathbf{p}),\mathbf{p})$ $\displaystyle=0,$ $\displaystyle\
\quad\quad\ \mathscr{V}(\mathbf{p})$ $\displaystyle=0,$

where $\mathbf{p}$ are the parameters that describe our geometry,
$\mathscr{V}$ is an equality constraint (e.g. for the volume), $\mathscr{R}$
is the PDE residual, $\bm{u}$ is the solution of the PDE and $J$ is the
objective. If desired, further equality and inequality constraints can then be
imposed by adding penalty terms to the objective function.

Update
---
Start
---
$J$
---
min?
---
End
---
$\mathbf{p}$
---
$\phi=\mathscr{H}(\varphi)$
---
$\mathbf{\varphi}:\mathbf{\varphi}=\bm{\mathcal{N}}(\mathbf{p})$
---
$\bm{u}:\mathscr{R}(\bm{u}(\phi),\phi)=0$
---
Y
---
N
---
Figure 1. Computational graph of the optimization loop. We start with an input
to the system $\mathbf{p}$, and perform the forward pass by descending through
the blue boxes on the right side to obtain a performance measure $J$ as
explained in Section 2.1.1. The convergence criteria are used to decide
whether this $J$ represents an acceptable minimum. If not, the backward pass
is performed to compute an update for the parameters as explained in Section
2.1.2 and the loop is continued.

###  2.1. Optimization Loop

To optimize the parameters $\mathbf{p}$, we make use of a gradient-based
optimization strategy. To do so, we establish a map between $\mathbf{p}$ and
$J$ and a means to compute the gradient $\frac{dJ}{d\mathbf{p}}$ for parameter
updates. For the neural LS TO method, $\mathbf{p}$ represents the parameters
of a particular NN that outputs a vector $\mathbf{\varphi}$. This vector is
processed using the operator $\mathscr{H}$ to obtain the LS $\phi$ used in the
PDE and objective function. It is also only through $\phi$ that the PDE and
objective depend on the parameter vector $\mathbf{p}$. With these definitions,
we present the optimization loop for solving (1) in Figure 1.

####  2.1.1. Forward Pass

To solve the forward problem and get a performance measure $J$ for a set of
parameters $\mathbf{p}$ we descend through the light blue boxes on the right-
hand side of Figure 1 by performing the following steps:

  1. (1)

In the first light blue box, we evaluate the network
$\bm{\mathcal{N}}:\mathbf{p}\in\mathbb{R}^{N_{p}}\mapsto\mathbf{\varphi}\in\mathbb{R}^{N}$,
where $N_{p}$ is the number of parameters and $\bm{\mathcal{N}}$ is as defined
in Section 3.

  2. (2)

In the second light blue box, we process the output of the network
$\mathbf{\varphi}$ using the operator
$\mathscr{H}:\mathbf{\varphi}\in\mathbb{R}^{N}\mapsto\phi\in V_{h}^{1}$ to
obtain a suitable LS description of the geometry $\phi$, where $V_{h}^{1}$ is
the FE space for the LS defined in Section 4. This involves an interpolation
on a FE space, smoothing and the inclusion of an equality constraint on the
geometry so that the final LS function satisfies $\mathscr{V}(\phi)=0$. This
step is broken down in Section 4. We enforce the equality constraint here to
allow for the use of an unconstrained optimization method suitable for NNs.

  3. (3)

In the third light blue box, we solve the FE problem associated with the weak
form of the residual $\mathscr{R}(\bm{u}_{h}(\phi),\phi)=0$ on the domain
segmentation defined by $\phi$ for an approximate solution $\bm{u}_{h}$
obtained using a FE discretization. Details of the FCM used to solve the
problem are given in Section 5. We can then evaluate the objective
$J(\bm{u}_{h}(\phi),\phi)\in\mathbb{R}$.

####  2.1.2. Backwards Pass

To perform the update for the parameters using a steepest descent optimization
strategy, we require an evaluation of the gradient $\frac{\partial
J}{\partial\mathbf{p}}$. To do this efficiently at a cost roughly matching
that of the forward pass, we use reverse mode differentiation and define rules
to propagate sensitivities through each of the steps in the forward pass. We
use an adjoint rule for the PDE and use automatic differentiation for all of
the partial derivatives, including the derivative of integrals with respect to
the LS. The derivative is then passed onto a chosen optimizer to update the
parameters $\mathbf{p}$. The implementation of the gradient computation is
discussed in Section 6. Using this method, arbitrary loss functions can be
used for the physical problem operating on the simulation output fields and
update directions for the neural network parameters can be computed to improve
the performance at each iteration with respect to the chosen loss function.

##  3\. Network Architecture

Dense

Convolution

$\bm{\Theta}$
---
$w_{n}$
---
$\mathbf{\varphi}$
---
$h_{n}$
---
$h_{1}$
---
$c_{1}$
---
$w_{1}$
---
$h_{2}$
---
$c_{2}$
---
$w_{2}$
---
Figure 2. Architecture of the NN. A trainable input vector $\bm{\Theta}$ is
fed into the network. The light blue arrow involves a set of operations that
include a fully connected layer and the dark blue arrow involves a set of
operations that include a convolutional layer. The intermediate data
structures are of size $(c_{l},w_{l},h_{l})$ and the final output, after $n$
layers, gives $\mathbf{\varphi}$.

In this section, we present the NN architecture used in the geometry
parameterization. This describes the mapping between the parameters
$\mathbf{p}$ and the LS vector $\mathbf{\varphi}$. The network used here is
based on the U-net architecture and is mainly built from convolutional layers.
In contrast to a fully connected network, convolutional networks connect
smaller sets of neurons in each layer assuming that neurons in close proximity
have a more important relationship. This is a natural relaxation for the
processing of spatial data. Furthermore, convolutional networks are made
efficient by the assumption that features that are found in one local block
are likely to be found in a different local block, i.e. somewhere else in the
domain. This is done by sharing parameters amongst local blocks in the form of
a convolutional filter. The specific architecture used in this work is the one
presented in [Hoyer2019]. Because we are simply reparameterizing the pixel
values, there is no input into the network in the traditional sense. Most
approaches combining NNs and TO use fully connected layers and take the
network to provide a map between a spatial input $\mathbf{x}$ and a scalar
output $f_{\theta}(\mathbf{x})$ [Deng2020, Chandrasekhar2020]. In our case, we
only need a single evaluation of the network to output a vector that
represents the entire discrete LS function. This is because, as in
[Hoyer2019], the dimension of the output of the network matches the number of
nodal values of the LS function so that we have a 1:1 mapping between the
network and the geometry discretization. This follows the typical use of the
U-net, where the output image gives a segmented domain [ronneberger2015unet].
In contrast, however, to original applications of the U-net, our input vector
is also taken to be a set of trainable parameters $\mathbf{\bm{\Theta}}$.

The architecture of the network is illustrated in Figure 2. The first arrow
indicated with the label Dense in Figure 2 contains a fully connected layer
and a reshape:

$\mathbf{x}^{(1)}(\bm{\Theta})=\mathrm{reshape}(\mathbf{W}\bm{\Theta}+\mathbf{b}),\quad\mathbf{x}^{1}\in\mathbb{R}^{c_{1},w_{1},h_{1}},$
(2)

where $\mathbf{W}\in\mathbb{R}^{c_{1}w_{1}h_{1},N_{\bm{\Theta}}}$ is the dense
weight matrix, $\mathbf{b}\in\mathbb{R}^{c_{1}w_{1}h_{1}}$ is the bias and
$\mathrm{reshape}:\mathbb{R}^{c_{1}w_{1}h_{1}}\rightarrow\mathbb{R}^{c_{1},w_{1},h_{1}}$
is a reshaping map. The lengths $N_{\bm{\Theta}}$, $c_{1}$, $w_{1}$ and
$h_{1}$ represent the length of $\bm{\Theta}$, the initial number of channels,
the latent space image width and the latent space image height, respectively.
Note here that the first nonlinearity is imposed at the beginning of the next
layer.

The upsampling convolutional layers, depicted by the dark blue arrow in Figure
2, are defined as:

$\mathbf{x}^{(l+1)}(\mathbf{x}^{(l)})=\bm{\mathcal{P}}^{(l)}(\bm{\varrho}(\mathbf{\Phi}^{(l)}(\text{tanh}(\mathbf{x}^{(l)})))),\quad\mathbf{x}^{(l)}\in\mathbb{R}^{c_{l},w_{l},h_{l}},\
\ \mathbf{x}^{(l+1)}\in\mathbb{R}^{c_{l+1},w_{l+1},h_{l+1}},$ (3)

where
$\mathbf{\Phi}^{(l)}:\mathbb{R}^{c_{l},w_{l},h_{l}}\rightarrow\mathbb{R}^{c_{l},w_{l+1},h_{l+1}}$
is a bilinear resize, $\bm{\varrho}$ is a normalization to a mean of $0$ and
variance of $1$ across the channel dimension and
$\bm{\mathcal{P}}^{(l)}:\mathbb{R}^{c_{l},w_{l+1},h_{l+1}}\rightarrow\mathbb{R}^{c_{l+1},w_{l+1},h_{l+1}}$
is a convolutional operator with kernel size $(5,5)$ which is found to allow
for the required level of expressivity.

In this approach, our input to the network $\bm{\Theta}$ is taken to be
trainable. So to define our parameter vector $\mathbf{p}$, we collect the
parameters of $\mathbf{W}$, $\mathbf{b}$, $\bm{\Theta}$ and $\bm{P}^{i}$ into
a vector $\mathbf{p}\in\mathbb{R}^{N_{p}}$. Then, by composing the layers, we
obtain the function
$\bm{\mathcal{N}}:\mathbf{p}\in\mathbb{R}^{N_{p}}\mapsto\mathbf{\varphi}\in\mathbb{R}^{N}$:

${\bm{\mathcal{N}}}=\mathbf{x}^{(n)}(...(\mathbf{x}^{(2)}(\mathbf{x}^{(1)}(\bm{\Theta})))).$
(4)

The output image of size $(w_{n},h_{n})$ obtained by evaluating the network at
a set of parameters $\mathbf{\varphi}=\bm{\mathcal{N}}(\mathbf{p})$ is then
used to define the geometry for the problem as described in Section 4.
Importantly, the locality is preserved when defining this function. As we
traverse the network, we follow the design principle of the upsampling section
of the U-Net and trade-off channel depth for spatial resolution. This means
that, in general, the widths $w$ and heights $h$ will increase as we move
through the network and the number of channels $c$ will decrease as we move
through the network. The exact trade-offs here can vary, but we must set the
resizes and channel refinement to ensure $h_{n}w_{n}=N$ and $c_{n}=1$ so that
our output $\mathbf{\varphi}$ makes sense as a vector representing the DOF
values for the LS. We can increase the network size in the dense layer by
increasing $N_{\bm{\Theta}}$, $w_{1}$ and $h_{1}$ and in the convolutional
layers by increasing the elements in $c$ and the number of layers.

One approach to initializing the parameters of the NN is to pre-train the
network to output the manually selected geometry with holes, as in [Deng2021].
Using this as an initial guess, however, causes the geometry to converge
quickly to poor local minima. A more common approach when working with NNs and
a given objective map is to start with small random weights [RUMELHART1988].
In this case, we have high asymmetry in the weights and little activation
function saturation. It turns out that a random initialization of the NN with
the volume constraint gives an initial guess of a domain with a few holes in
random locations. This is in contrast to initializing the LS FE function DOF
values with random values which gives a geometry with many small holes and
fine features which is not necessarily desirable [Barrera2020]. As is common
in NN approaches, we can then easily take multiple seeds for the geometry
using different random initializations to alleviate initialization dependency.
The random initializations simply correspond to different size holes in
different locations. Using this method we eliminate the need for manual
geometry initialization.

We use Xavier uniform initialization for the weights
[Glorot2010UnderstandingTD], in this case the design variables, of the neural
network to facilitate convergence. To minimize the effect of the
initialization dependency of the optimized geometry, that is, wherever the
holes appear in the initial design based on the particular random
initialization, we take 100 random seeds of the parameters, run the
optimization and take the final design with the lowest objective function
value.

##  4\. LS Function Processing

In this section, we detail the LS description of the geometry used in the
method. We explain how the output vector from the NN $\mathbf{\varphi}$ is
used to define the LS to be used in the numerical method. The computation of
the LS function is performed in four steps. We introduce an interpolation step
in Section 4.1 to obtain a first LS function $\phi_{n(1)}$. In Section 4.2, we
smooth out this LS to obtain a smooth LS function $\phi_{f(2)}$. Then, we
perform a reinitialization of that LS in Section 4.3 to obtain an LS function
$\phi_{s(3)}$. Finally, in Section 4.4, we propose a volume correction
strategy to end up with the final LS function $\phi_{b(4)}$. In the subsequent
sections, we will refer to this final LS function $\phi_{b(4)}$ as $\phi$.

###  4.1. Interpolation

In this method, we work with discrete LS functions $\phi\in V_{h}^{1}$. The
DOF of $V_{h}^{1}$ are the values of the function at the vertices of
$\mathcal{T}_{h}$, which we denote with $\\{\bm{x}_{i}\\}_{i=1}^{N}$, where
$N$ is the number of mesh nodes. Thus, there is an isomorphism between
$V_{h}^{1}$ and $\mathbb{R}^{N}$. The output image vector of the NN $\varphi$
are the DOF values that uniquely determine the LS FE function.

###  4.2. Smoothing

Next, we convolve the function with a linear filter for smoothing:

${\phi_{f(2)}}_{i}=(\sum^{N}_{j=1}w_{ij})^{-1}(\sum^{N}_{j=1}w_{ij}{\phi_{n(1)}}_{j}),$
(5)

where $w_{ij}=\max(0,r_{f}-|x_{i}-x_{j}|)$, ${\phi_{n(1)}}_{j}$ is the
$j^{th}$ degree of freedom of the function $\phi_{n(1)}$, ${\phi_{f(2)}}_{i}$
is the $i^{th}$ degree of freedom of the function $\phi_{f(2)}$ and $r_{f}$ is
the smoothing radius.

###  4.3. Reinitialization

We then reinitialize the LS as a signed distance function. This is often done
to gain control over the spatial gradient of the LS function to improve
convergence [vanDijk2013]. In our case, it is of even greater importance as it
also guarantees that when we apply the translation to the LS to satisfy the
volume constraint, we do not artificially introduce volumes into the domain
far from the boundary where the LS is close to zero. This would add
discontinuity to the problem harming convergence. To perform this step, we
solve the reinitialization equation in [Xing2009] to obtain a signed distance
function $\phi_{s(3)}\in V_{h}^{1}$:

$\frac{\partial\phi_{s(3)}}{\partial\tau}+\text{sign}(\phi_{s(3)})(|\bm{\nabla}\phi_{s(3)}|-1)=0.$
(6)

The problem is solved using Picard iterations at steady state using the
initial point $\phi_{f(2)}$ so that only one solve of the adjoint equation is
required in the backward pass. Artificial viscosity is added to the problem
for stabilization and a surface penalty term is integrated on the embedded
boundary to prevent movement of the zero iso-surface in the reinitialization.
The weak form of the problem is then to find the solution $\phi_{s(3)}\in
V_{h}^{1}$ of the equation:

$\int_{\Omega}(\bm{w}\cdot\bm{\nabla}(\phi_{s(3)})v+c_{a}h|\bm{w}|\bm{\nabla}(\phi_{s(3)})\bm{\nabla}v-\text{sign}(\phi_{s(3)})v)d\Omega+\int_{\Gamma}(\phi_{s(3)},v)d\Gamma=0$
(7)

for all $v\in V_{h}^{1}$, where
$\bm{w}=\text{sign}(\phi_{s(3)})(\bm{\nabla}(\phi_{s(3)})/|\bm{\nabla}(\phi_{s(3)})|)$,
$c_{a}$ is the stabilization coefficient, set to $3$, and $h$ is the element
size.

###  4.4. Translation

The next task is to impose the volume constraint. This is done by applying a
translation to the entire LS function. Here we solve the nonlinear equation
$\mathscr{V}$ for the scalar bias $b$ to obtain a translation to the LS
function such that the domain $\Omega_{\rm in}$ satisfies the volume
constraint. A bisection method is used to find the root $b$ of the equation:

$\mathscr{V}=\int_{\Omega_{\rm in}}d\Omega-\mathscr{V}_{0}=0,$ (8)

where $\Omega_{\rm in}$ is defined by the LS $\phi_{b(4)}=\phi_{s(3)}+b$ and
$\mathscr{V}_{0}$ is the volume fraction given by the constraint. This is the
LS which is then used to define the boundary of the FE problem. Details of how
the LS is used to define domains is given in the next Section.

##  5\. Numerical Discretization

$\Omega_{\rm in}$
---
$\partial\Omega$
---
$\mathcal{R}(\bm{u},\bm{v})$
---
$\mathcal{R}^{\mathrm{stb}}(\bm{u},\bm{v})$
---
$\Omega_{\rm out}$
---
$\Gamma$
---
Figure 3. Illustration of the problem domains. The background domain is
segmented into the domains $\Omega_{\rm in}$ and $\Omega_{\rm out}$. In
$\Omega_{\rm in}$, the physical terms $\mathcal{R}(\bm{u},\bm{v})$ are
integrated. In $\Omega_{\rm out}$, the stabilization terms
$\mathcal{R}^{\mathrm{stb}}(\bm{u},\bm{v})$ are integrated.

Here we formulate the numerical discretization of the PDEs being solved in the
simulation step of the optimization loop. We consider a background polyhedral
bounded domain $\Omega$ with boundary $\partial\Omega$. A LS function
$\phi(x)$ is used to split $\Omega$ into two subdomains $\Omega_{\rm
in}(\phi)$ and $\Omega_{\rm out}(\phi)$ as follows (see Figure 3):

$\displaystyle\Omega_{\rm
in}(\phi)=\\{x\in\Omega:\phi(x)>0\\},\quad\Omega_{\rm
out}(\phi)=\\{x\in\Omega:\phi(x)<0\\}.$ (9)

We denote the interface between these two subdomains as
$\Gamma(\phi)\doteq\partial\Omega_{\mathrm{in}}(\phi)\cap\partial\Omega_{\mathrm{out}}(\phi)$.
Let us assume that
$\partial\Omega_{\mathrm{in}}(\phi)\cap\partial\Omega\neq\emptyset$. The
domain $\Omega_{\mathrm{in}}(\phi)$ is the one in which we consider our PDE
problem. The weak form of the continuous problem can be stated as follows:
find $\bm{u}\in V$ such that

$\int_{\Omega_{\mathrm{in}}(\phi)}\mathcal{R}(\bm{u},\bm{v})\mathrm{d}\Omega=0,\quad\forall
v\in V.$ (10)

where $V$ is a Hilbert space in which the problem is well-posed. We consider
zero flux Neumann boundary conditions on
$\partial\Omega_{\mathrm{in}}(\phi)\setminus\partial\Omega$ and (for
simplicity) homogeneous boundary conditions on
$\partial\Omega_{\mathrm{in}}(\phi)\cap\partial\Omega$; the generalization to
non-homogeneous Dirichlet boundary conditions is straightforward.

The domain $\Omega_{\mathrm{in}}(\phi)$ will change along the optimization
process. As a result, it is not practical to compute body-fitted unstructured
meshes for the geometrical discretization of $\Omega_{\mathrm{in}}(\phi)$.
Instead, we consider a background mesh, which can simply be a Cartesian
background mesh of $\Omega$ and make use of an unfitted FE discretization.
Unfitted (or embedded) discretizations relax the geometrical constraints but
pose additional challenges to the numerical discretization [Burman2010,
Badia2018]. The first issue is the integration over cut cells (for details,
see [Badia2022-Geometrical]). The other issue is the so-called _small cut
cell_ problem. Cut cells with arbitrary small support lead to ill-conditioned
systems [dePrenter2017]. Various techniques can be used to stabilize the
problem including the cut finite element method (CutFEM) [Burman2010], the
aggregated finite element method (AgFEM) [Badia2018] and the FCM
[Parvizian2007]. We select the FCM method in this case because it is
differentiable with respect to the LS everywhere in the domain. The CutFEM and
AgFEM are conversely not differentiable.

In order to state the discrete form of the continuous problem, we introduce
the unfitted FE space. Let $\mathcal{T}_{h}$ represent a conforming, quasi-
uniform and shape regular partition (mesh) of $\Omega$, $h$ being a
characteristic mesh size. $\Omega$ can be a trivial geometry, e.g., a square
or cube, and $\mathcal{T}_{h}$ can be a Cartesian mesh. We define a nodal
Lagrangian FE space of order $q\geq 1$ on $\mathcal{T}_{h}$ as:

$V_{h}^{q}=\\{\bm{v}_{h}\in\mathcal{C}^{0}(\Omega):\bm{v}_{h}|_{K}\in\mathcal{X}_{q}(K)\
\forall K\in\mathcal{T}_{h}\\},$ (11)

where $\mathcal{X}_{q}(K)$ is the space $\mathcal{Q}_{q}(K)$ of polynomials
with maximum degree $q$ for each variable when $\mathcal{T}_{h}$ is a
quadrilateral or hexahedral mesh and the space $\mathcal{P}_{q}(K)$ of
polynomials of total degree $q$ when $\mathcal{T}_{h}$ is a simplicial mesh.
In this work, we consider low-order spaces, which is the most reasonable
choice for TO applications.

The weak formulation of the problem solved by the method is now described. Let
us represent with $V_{h,0}^{q}=V_{h}^{q}\cap\mathcal{C}_{0}^{0}(\bar{\Omega})$
the nodal FE space that vanishes on the boundary $\partial\Omega$.

Now, we can define a first-order FCM discretization of (10) as follows: find
$\bm{u}_{h}\in V_{h}^{1}$ such that

$\int_{\Omega_{\mathrm{in}}(\phi)}\mathcal{R}(\bm{u}_{h},\bm{v}_{h})\mathrm{d}\Omega+\int_{\Omega\setminus\Omega_{\mathrm{in}}(\phi)}\alpha_{\mathrm{out}}\mathcal{R}^{\mathrm{stb}}(\bm{u}_{h},\bm{v}_{h})\mathrm{d}\Omega=0,\quad\forall
v\in V_{h}^{1},$ (12)

where $\alpha_{\mathrm{out}}\ll 1$ is the penalty parameter and
$\mathcal{R}^{\mathrm{stb}}$ is a _stabilizing_ differential operator on the
artificial domain.

###  5.1. Poisson Equation

The poisson equation is used to model the temperature $\theta(\bm{x})$ that
satisfies

$\hypertarget{eq:T}{}{-\bm{\nabla}\cdot(\bm{\kappa}\bm{\nabla}\theta)=f\quad\text{in}\
\Omega_{\mathrm{in}}(\phi),}$ (13)

where $\bm{\kappa}$ is the thermal conductivity tensor of the material and $f$
is a thermal source. A zero Dirichlet condition ($\theta=0$) is prescribed on
$\Omega_{\mathrm{in}}(\phi)\cap\partial\Omega$ and a zero flux condition
($\bm{n}\cdot(\bm{\kappa}\bm{\nabla}\theta)=0$) is prescribed on
$\Gamma(\phi)$.

The FCM approximation of this problem without a source term reads as: find
$\theta_{h}\in V_{h,0}^{1}$ such that

$\int_{\Omega_{\mathrm{in}}(\phi)}\bm{\kappa}\bm{\nabla}\theta_{h}\cdot\bm{\nabla}v_{h}\mathrm{d}\Omega+\int_{\Omega\setminus\Omega_{\mathrm{in}}(\phi)}\alpha_{\mathrm{out}}\bm{\kappa}\bm{\nabla}\theta_{h}\cdot\bm{\nabla}v_{h}\mathrm{d}\Omega=0,\quad\forall
v_{h}\in V_{h,0}^{1}.$ (14)

One can readily check that this method is weakly enforcing the zero flux
condition on $\Gamma(\phi)$ as $\alpha_{\mathrm{out}}\rightarrow 0$. We
observe that we use the same differential operator in the artificial domain
for stabilization purposes (times the scaling coefficient
$\alpha_{\mathrm{out}}$).

We consider the TO problem in which we aim at finding a level-set $\phi$ that
minimizes the integral of the temperature:

$J(\phi,\theta_{h}(\phi))=\int_{\Omega(\phi)}\theta_{h}(\phi)\
\mathrm{d}\Omega,$ (15)

where $\theta_{h}(\phi)$ is the solution of (14) given $\phi$.

###  5.2. Linear elasticity

We want to obtain the displacement $\bm{d}(\bm{x})$ that satisfies the linear
elasticity equation

$\hypertarget{eq:d}{}{\left.\begin{aligned}
-\bm{\nabla}\cdot\bm{\sigma}(\bm{d})&=\bm{f}&&\text{in}\ \Omega_{\rm
in}(\phi),\end{aligned}\right.}$ (16)

where $\bm{\sigma}=\lambda\text{tr}(\bm{\varepsilon})I+2\mu\bm{\varepsilon}$
is the stress tensor,
$\bm{\varepsilon}=\frac{1}{2}(\bm{\nabla}\bm{d}+(\bm{\nabla}{\bm{d}})^{\top})$
is the symmetric gradient,  $\ignorespaces\bm{d}$ is the displacement,
$\lambda$ and $\mu$ are the Lamé parameters given by
$\lambda=(E\nu)/((1+\nu)(1-2\nu))$ and $\mu=E/(2(1+\nu))$ and $\bm{f}$ is the
forcing term. A zero Dirichlet condition ($\bm{d}=\bm{0}$) is prescribed on
$\Omega_{\mathrm{in}}(\phi)\cap\partial\Omega$ and a zero stress condition
($\bm{n}\cdot\bm{\sigma}(\bm{d})=0$) is prescribed on $\Gamma(\phi)$.

The FCM approximation of this problem without a forcing term reads as: find
$\bm{d}_{h}\in\bm{V}_{h,0}^{1}\doteq[V_{h,0}^{1}]^{D}$ such that

$\int_{\Omega_{\mathrm{in}}(\phi)}\bm{\sigma}(\bm{d}_{h}):\bm{\varepsilon}(\bm{v}_{h})\mathrm{d}\Omega+\int_{\Omega\setminus\Omega_{\mathrm{in}}(\phi)}\alpha_{\mathrm{out}}\bm{\sigma}(\bm{d}_{h}):\bm{\varepsilon}(\bm{v}_{h})\mathrm{d}\Omega=0,\quad\forall\bm{v}_{h}\in\bm{V}_{h,0}^{1}.$
(17)

It is easy to check that the zero-stress condition on $\Gamma(\phi)$ is
recovered as $\alpha_{\mathrm{out}}\rightarrow 0$. For the linear elasticity
equation, we again use the same differential operator in the artificial domain
for stabilization purposes.

A typical TO problem in solid mechanics is the minimization of the strain
energy. In this case, we aim at finding a level-set $\phi$ that minimizes the
cost function

$J(\phi,\bm{d}_{h}(\phi))=\int_{\Omega(\phi)}\bm{\sigma}(\bm{d}_{h}):\bm{\varepsilon}(\bm{d}_{h})\
\mathrm{d}\Omega,$ (18)

where $\bm{d}_{h}(\phi)$ is the solution of (17) given $\phi$.

###  5.3. Linear Elasticity with Fluid Forcing Terms

Once again, we want to obtain the displacement $\bm{d}_{h}(\bm{x})$ that
satisfies the linear elasticity formulation in (16). In this case, however, we
consider the surface traction exerted by the fluid:

$\displaystyle\int_{\Gamma(\phi)}(\bm{n}\cdot\bm{\nabla}\bm{u}_{h}-{p}_{h}\
\bm{n})\cdot\bm{v}{\rm d}x,$ (19)

where the fluid velocity $\bm{u}_{h}$ and pressure field ${p}_{h}$ are
obtained by solving a fluid problem in the domain $\Omega_{\rm out}$. These
fields are obtained by solving the Stokes equations with a Brinkmann
penalization as in [Borrvall2002] but without intermediate interpolation of
permeabilities at the boundary.

In order to approximate the fluid problem, we use a mixed FE method, namely
the equal order pair $\bm{V}_{h,0}^{1}\times{V}_{h}^{1}$.

We find  $\ignorespaces(\bm{u}_{h},{p}_{h})\in\bm{V}_{h,0}^{1}\times
V_{h,0}^{1}$ such that:

$\displaystyle\int_{\Omega}$
$\displaystyle\left[\alpha\bm{u}_{h}\cdot\bm{\psi}_{h}+\mu\bm{\nabla}\bm{u}_{h}\cdot\bm{\nabla}\bm{\psi}_{h}-{p}_{h}(\bm{\nabla}\cdot\bm{\psi}_{h})-(\bm{\nabla}\cdot\bm{u}_{h}){q}_{h}-h^{2}\bm{\nabla}p_{h}\cdot\bm{\nabla}{q_{h}}\right]{\rm
d}x=0,\quad$ (20)
$\displaystyle\forall\bm{\psi}_{h},{q}_{h}\in\bm{V}_{h,0}^{1}\times{V}_{h}^{1},$

where

$\left\\{\begin{aligned} \alpha&=0&&\text{in }\Omega_{\rm out}\\\
\alpha&=\alpha_{u}&&\text{in }\Omega_{\rm in}\\\ \end{aligned}\right.$ (21)

using an artificial porosity $\alpha_{u}$ to make the fluid problem well-posed
in the solid domain $\Omega_{\rm in}$ and enforce the no-slip boundary
condition. In the fluid domain $\Omega_{\rm out}$ we recover the Stokes
equations.

The TO problem once again involves finding a level-set $\phi$ that minimizes
the elastic strain using (18).

###  5.4. Differentiability of the unfitted FE solver

An important property for the convergence of a TO strategy is the notion of
shape differentiability of the cost function. A functional under a PDE
constraint is considered shape differentiable if the mapping $\phi\rightarrow
J(\bm{u_{h}}(\phi),\phi)$ is differentiable at the admissible set of domains
in $\Omega$ defined by $\phi$. In this section, we discuss how the choice of
FE stabilization can affect this property.

The model problem in (14) with a FCM stabilization is equivalent to that of a
typical two-phase conductivity problem. With a solution $\bm{u}_{h}\in
H^{1}(\Omega)$, the functional $J(\bm{u_{h}},\phi)$ for this problem can be
proven to be shape differentiable, see [Allaire2021, Theorem 4.9].

Conversely, unfitted techniques involving stabilization only in the vicinity
of the boundary are in general not shape differentiable. Regions of non-
differentiability arise (typically when the boundary crosses over mesh nodes)
harming the convergence of the geometry to optimized solutions [Sharma2016].
To see why this is the case, we investigate shape perturbations under the
CutFEM and AgFEM formulations.

$\Omega_{\rm in}^{\epsilon}$
---
$\epsilon$
---
Figure 4. A small perturbation with size $\epsilon$ to the LS to form a new
domain $\Omega_{\rm in}^{\epsilon}$.

If we were to use a restricted space for the solution and add ghost penalty
terms in the vicinity of the interface following the CutFEM method
[CutFEM2015], the problem is to find $u_{h}\in W_{h,0}^{1}.$ such that:

$\displaystyle\int_{\Omega_{\mathrm{in}}(\phi)}\mathcal{R}(\bm{u}_{h},\bm{v}_{h})\mathrm{d}\Omega_{\Omega_{\rm
in}}+j(\bm{u}_{h},\bm{v}_{h})_{\Gamma_{G}}=0,\quad\forall v_{h}\in
W_{h,0}^{1},$ (22)

for a ghost penalty term $j$ on a ghost skeleton triangulation $\Gamma_{G}$
using the space $W_{h,0}^{1}$ as in [CutFEM2015]. Consider the change to the
domain $\Omega_{\rm in}$ from Figure 3 caused by a perturbation
$\epsilon\delta$, where $\epsilon\in\mathbb{R}$ and $\delta\in V^{1}_{h}$, to
the LS function $\phi$. The resulting domain $\Omega_{\rm in}^{\epsilon}$ may
be as in Figure 4. The ghost penalty term in this formulation does not depend
on $\epsilon$ and instead changes depending on which cells are cut.
Specifically, if $\Gamma^{\epsilon}$ crosses over a mesh node as in Figure 5,
the ghost skeleton triangulation includes the faces of a new element. Non-zero
terms are integrated on this triangulation introducing discontinuity to the
problem with respect to the shape, since:

$\lim_{\epsilon\to 0}\
[j(\bm{u}_{h},\bm{v}_{h})_{\Gamma_{G}^{\epsilon}}-j(\bm{u}_{h},\bm{v}_{h})_{\Gamma_{G}}]\
\neq 0.$ (23)

With different terms added to the linear system which do not go to zero with
$\epsilon$, we can see that the derivative of the solutions with respect to
the shape can be ill-defined. A cost function operating on the solution could
not, in general, be shape differentiable at these points.

$\Omega_{\rm in}$
---
$\Omega_{\rm out}$
---
$\Gamma$
---
$\Gamma_{G}$
---
$\Gamma_{G}^{\epsilon}$
---
$\Omega_{\rm in}^{\epsilon}$
---
$\Omega_{\rm out}^{\epsilon}$
---
$\Gamma^{\epsilon}$
---
Figure 5. The change in a portion of the ghost skeleton triangulation
$\Gamma_{G}$, depicted by the red faces, before and after a perturbation to
the boundary.

In the AgFEM, the problem is to find $u_{h}\in V^{agg}_{h,0}$ such that:

$\displaystyle\mathscr{R}(\mathcal{E}(\bm{u}_{h}),\mathcal{E}(\bm{v}_{h}))_{\Omega_{\rm
in}}\quad\forall v_{h}\in V^{agg}_{h,0},$ (24)

for the extension operator $\mathcal{E}$ and space $V^{agg}_{h,0}$ as defined
in [Badia2018]. Similar to the CutFEM, regions of non-differentiability in the
problem exist when the zero iso-surface of the LS crosses over mesh nodes:

$\displaystyle\lim_{\epsilon\to 0}\
[\mathscr{R}(\mathcal{E}(\bm{u}_{h}),\mathcal{E}(\bm{v}_{h}))_{\Omega_{\rm
in}}-\mathscr{R}(\mathcal{E}^{\epsilon}(\bm{u}_{h}),\mathcal{E}(\bm{v}_{h}))_{\Omega_{\rm
in}^{\epsilon}}]\neq 0,$ (25)

since $\mathcal{E}^{\epsilon}(\bm{u})\neq\mathcal{E}(\bm{u}_{h})$ in general
because the support for cut cells potentially changes depending on which cells
are cut. Following the same reasoning as above, the method is therefore not
shape differentiable. The same is true for other methods which use
stabilization approaches that act only in the vicinity of cut cells, e.g.
[Lang2014], using similar branching strategies when crossing a node or
reaching a certain threshold.  $\square$

##  6\. Gradient Implementation

To the best of our knowledge, there is no existing implementation of an
unfitted LS TO method that accepts arbitrary residuals defining the PDE and
computes the entire gradient $\frac{dJ}{d\mathbf{p}}$ by automatic
differentiation. Making use of a backward pass, we do this efficiently by
defining differentiation rules for each of the steps in the method.

###  6.1. Integral Differentiation Operator

The backward pass is mainly composed of gradients of integrals with respect to
the DOF of FE functions. To make the derivative computation efficient, we
exploit the fact that the DOF only have an effect on surrounding cells and
utilize the optimizations exploiting sparsity in the FE library Gridap
[Verdugo2021]. Integrals in the domain can be divided into cell-wise
components:

$\mathscr{I}(\bm{u},\bm{v},\phi)=\sum_{K\in\mathcal{T}_{h}}\mathscr{I}^{K}(\mathbf{u}^{K},\mathbf{v}^{K},\mathbf{\phi}^{K}),$
(26)

where
$\mathbf{u}^{K}\in\mathbb{R}^{\Sigma_{u}},\mathbf{v}^{K}\in\mathbb{R}^{\Sigma_{v}}$
and $\phi^{K}\in\mathbb{R}^{\Sigma_{\phi}}$ are the DOF parameterizing the
restrictions of $u,v$ and $\phi$ to the cell $K$ and $\Sigma u$, $\Sigma v$
and $\Sigma\phi$ are the number of DOF in $K$ for the respective functions.
Gradients can then be computed at roughly the cost of an integral evaluation
for each cell $K$:

$\frac{\partial\mathscr{I}}{\partial\phi}^{K}=\nabla^{F}_{\phi}\mathscr{I}^{K}(\mathbf{u}^{K},\mathbf{v}^{K},\phi^{K}),$
(27)

where the operator $\nabla^{F}_{\phi}$ represents taking the gradient with
respect to $\phi^{K}$ using a vectorized forward propagation of dual numbers
[ForwardDiff]. To make taking derivatives in this way possible for the LS, we
implement the integrals so that each $\phi^{K}$ is accepted as the argument to
compute the contribution $\mathscr{I}^{K}$:

$\mathscr{I}^{K}:\phi^{K}\in\mathbb{R}^{\Sigma_{\phi}}\mapsto\mathscr{I}_{K}(\mathbf{u}^{K},\mathbf{v}^{K},\phi^{K})\in\mathbb{R}.$
(28)

where

$\mathscr{I}^{K}(\mathbf{u}^{K},\mathbf{v}^{K},\phi^{K})=\int_{K(\phi^{K})}\mathcal{I}(\mathbf{u}^{K},\mathbf{v}^{K})dK$
(29)

A key point is that the integral function subroutines, including all the
unfitted FE tools, are implemented in such a way as to allow the propagation
of dual numbers through the code. We also make use of a reverse mode operator
$\nabla^{R}$ for the backwards propagation of derivatives used where
appropriate, e.g., for the NN.

###  6.2. Backwards Pass Routine

We now present the backward pass in detail. To compute the sensitivity of the
objective with respect to the parameters, we start with the seed
$\frac{dJ}{dJ}=1$ and propagate derivatives in reverse mode using the chain
rule:

$\frac{dJ}{d\mathbf{p}}=\frac{dJ}{dJ}\left(\frac{\partial{J}}{\partial{\phi}}+\frac{\partial{J}}{\partial{\bm{u}}}\frac{d\bm{u}}{d\phi}\right)\frac{d\phi}{d\mathbf{\varphi}}\frac{d\mathbf{\varphi}}{d{\mathbf{p}}}$
(30)

where an adjoint method on the problem residual $\mathscr{R}$ is used to
differentiate through the PDE:

$\frac{\partial{J}}{\partial{\bm{u}}}\frac{d\bm{u}}{d\phi}={-\lambda^{T}\frac{d\mathscr{R}}{d\phi}}\
(\text{ here we solved
}\frac{d\mathscr{R}}{d\bm{u}}^{T}\lambda=\frac{dJ}{d\bm{u}}^{T}).$ (31)

We then use the chain rule to differentiate through the LS function processing
steps:

$\frac{d\phi}{d\mathbf{\varphi}}=\frac{d\phi}{d{\phi_{s(3)}}}\frac{d{\phi_{s(3)}}}{d{\phi_{f(2)}}}\frac{d{\phi_{f(2)}}}{d\mathbf{\varphi}}.$
(32)

The volume constraint here involved a root-finding method. To differentiate
through this step, we utilize the implicit function theorem:

$\frac{d\phi}{d{\phi_{s(3)}}}=\frac{\partial\phi}{\partial{\phi_{s(3)}}}-\frac{\partial\phi}{\partial
b}\frac{\partial\mathscr{V}}{\partial
b}^{-1}\frac{\partial\mathscr{V}}{\partial{\phi_{s(3)}}},$ (33)

and to differentiate through the signed distance map, we use the adjoint
method once again for the residual $\mathscr{R}_{s}$ equal to the integral in
(7). Finally, we use standard backpropagation to compute the derivative with
respect to the parameters of the NN. The steps of the backward pass are
presented explicitly in Algorithm 1.

Algorithm 1 Backwards Pass

Initialize $\frac{dJ}{dJ}\leftarrow 1$

Extract $\phi^{K}\in\mathbb{R}^{\Sigma\phi}$,
$\mathbf{u}^{K}\in\mathbb{R}^{\Sigma u}$ from $\phi$,$\bm{u}$
$\forall{K}\in\mathcal{T}_{h}$.

for $K\in\mathcal{T}_{h}$ do

$\frac{\partial
J}{\partial\bm{u}}^{K}\leftarrow\nabla^{F}_{u}J(\mathbf{u}^{K},\phi^{K})$

$\frac{\partial
J}{\partial\phi}^{K}\leftarrow\nabla^{F}_{\phi}J(\mathbf{u}^{K},\phi^{K})$

end for

Assemble the gradients $\frac{\partial
J}{\partial\bm{u}}\in\mathbb{R}^{N_{u}}$ and $\frac{\partial
J}{\partial\phi}\in\mathbb{R}^{N}$

Assemble the sparse jacobian associated with the residual
$\frac{\partial\mathscr{R}}{\partial\bm{u}}\in\mathbb{R}^{N_{u},N_{u}}$

Solve the adjoint equation
$\frac{\partial\mathscr{R}}{\partial\bm{u}}\lambda=\frac{\partial
J}{\partial\bm{u}}$ for $\lambda\in\mathbb{R}^{N_{u}}$

Extract $\lambda^{K}\in\mathbb{R}^{\Sigma u}$ from $\lambda$ $\forall
K\in\mathcal{T}_{h}$

for $K\in\mathcal{T}_{h}$ do

$\frac{\partial
J}{\partial\bm{u}}\frac{\partial\bm{u}}{\partial\phi}^{K}\leftarrow\nabla^{F}_{u}\mathscr{R}(\mathbf{u}^{K},\lambda^{K},\phi^{K})$

end for

Assemble the gradient $\frac{\partial
J}{\partial\bm{u}}\frac{\partial\bm{u}}{\partial\phi}\in\mathbb{R}^{N}$

$\frac{dJ}{d\phi}\leftarrow\frac{\partial J}{\partial\phi}+\frac{\partial
J}{\partial\bm{u}}\frac{\partial\bm{u}}{\partial\phi}$

Compute the vector-jacobian-products:

$\frac{\partial J}{\partial
b}\leftarrow\frac{dJ}{d\phi}\nabla^{R}_{b}\phi(\phi_{s(3)},b)$

$\frac{\partial
J}{\partial\phi}\leftarrow\frac{dJ}{d\phi}\nabla^{R}_{\phi}\phi(\phi_{s(3)},b)$

Compute the gradients:

$\frac{\partial V}{\partial\phi}\leftarrow\nabla^{R}_{u}V(\phi_{s(3)},b)$

$\frac{\partial V}{\partial b}\leftarrow\nabla^{F}_{u}V(\phi_{s(3)},b)$

$\frac{dJ}{d{\phi_{s(3)}}}\leftarrow\frac{\partial
J}{\partial{\phi_{s(3)}}}-\frac{\partial J}{\partial
b}\frac{\partial\mathscr{V}}{\partial
b}^{-1}\frac{\partial\mathscr{V}}{\partial{\phi_{s(3)}}}$

Assemble the sparse jacobian associated with the residual
$\frac{\partial\mathscr{R}_{s}}{\partial\phi_{s(3)}}\in\mathbb{R}^{N,N}$

Solve the adjoint equation
$\frac{\partial\mathscr{R}_{s}}{\partial\phi_{s(3)}}\lambda_{s}=\frac{\partial
J}{\partial\phi_{s(3)}}$ for $\lambda_{s}\in\mathbb{R}^{N}$

for $K\in\mathcal{T}_{h}$ do

$\frac{d\mathscr{R}_{s}}{d\phi_{f(2)}}^{K}\leftarrow\nabla^{F}_{\phi_{f(2)}}\mathscr{R}_{s}(\phi_{s(3)}^{K},\lambda_{s}^{K},\phi_{f(2)}^{K})$

end for

Assemble the gradient $\frac{dJ}{d\phi_{f(2)}}\in\mathbb{R}^{N}$

Compute the vector-jacobian-products:

$\frac{dJ}{d\varphi}\leftarrow\frac{dJ}{d\phi_{f(2)}}\nabla^{R}_{\varphi}(\phi_{f(2)}(\varphi))$

$\frac{dJ}{dp}\leftarrow\frac{dJ}{d\varphi}\nabla^{R}_{p}(N(\mathbf{p}))$

##  7\. Numerical Experiments

###  7.1. Benchmark Results

We first compare the optimized results obtained for benchmark problems against
baseline methods using the model problems presented in Section 5.

The method presented in this work, the NN-LS method, is compared against its
non-neural counterpart the Pixel LS (Pixel-LS) method, that is, the same
method without a neural prior where instead the nodal values of the level set
function are taken as the optimization parameters. We also compare the method
against the SIMP method of TO again using both a neural prior (NN-SIMP) and
the standard approach (Pixel-SIMP). Following a standard SIMP implementation
of the heat conduction problem, we use a conductivity based on the power law
$k=\alpha_{T}+(1-\alpha_{T})\rho^{\gamma}$ where $\rho$ is a design variable
given by a FE function constructed on the space $V_{h}^{1}$ where the DOF
values are the optimization parameters or output vector of the NN in the
Pixel-SIMP and NN-SIMP cases, respectively, and $\gamma$ is the penalization
parameter, taken to be equal to $3$. Similarly, for the SIMP implementation of
the structural problem, we use a Youngs Modulus
$E=\alpha_{d}+(1-\alpha_{T})\rho^{\gamma}$. To maintain a fair comparison
between methods, we use the optimizers most widely regarded as suitable for
the particular parameterization. Namely, we take the most commonly used MMA
optimization strategy [Svanberg1987] for the pixel parameterization and the
ADAM strategy [Kingma2014] for the NN.

####  7.1.1. Benchmark Problems

$\Omega$
---
$\Gamma_{N}$
---
$\Gamma_{D}$
---
$\Gamma_{D}$
---
$\Gamma_{N}$
---
(a) Heat conduction problem setup. The top and left sides $\Gamma_{D}$ are
given a Dirichlet condition and the bottom and left sides $\Gamma_{N}$ are
given a Neumann condition.

$\Omega$
---
$F$
---
(b) The right half of the MBB problem exploiting symmetry. The roller supports
provide vertical restraint on the left-hand side and horizontal restraint in
the bottom right corner. A downward force $F$ is prescribed on the top left
corner.

Figure 6. Benchmark problems

The first problem studied is the Poisson equation to model heat conduction
with the setup in Figure 6(a) selected from [GersborgHansen2006]. For this
problem, we set $k_{0}=1m^{2}s^{-1}$, $\alpha_{T}=0.01$, $f=0.01Ks^{-1}$, use
homogenous Dirichlet and Neumann conditions and use a mesh of $95\times 95$
with a $0.4$ volume fraction. For the NN, we set the number of convolutional
layers to $5$, $N_{\bm{\Theta}}=64$, $w=(12,12,24,46,96,96)$,
$l=(12,12,24,46,96,96)$ and $c=(16,128,64,32,16,1)$.

The second problem studied is the typical MBB problem described in
[Sigmund2001] with the setup as in Figure 6(b). For this problem, we set
$\nu=0.3$, $E=1Pa$, $\alpha_{d}=0.001$ and $F=1N$ and use a mesh of $287\times
95$ with a $0.4$ volume fraction. We use the same network as in the Poisson
problem but set $w=(36,36,64,128,256,256)$

The hyperparameters of the network are mainly selected to match those
presented in [Hoyer2019] with minor adjustments for the specific problem. The
size of the convolutional filter chosen is appropriate for this scale of
problem since we have achieved the required level of expressivity, as can be
seen in the subsequent results, through learning complex enough templates
whilst retaining a reasonable number of parameters.

####  7.1.2. Optimized Structures

To compare the SIMP and LS results, the optimized densities from the SIMP
procedure are converted to an LS by taking the 0.5 iso-surface of the density
and recomputing the objective function value using the unfitted machinery
which integrates exactly on the cut cells. To evaluate the effect of the NN
parameterization, we also include the results obtained through the more
standard methods of parameterization by taking the coefficients of the FE
function defining the geometry directly as the design variables.

| NN-LS | Pixel-LS | NN-SIMP | Pixel-SIMP
---|---|---|---|---
Heat |  |  |  |
| 0.0% | 20.3% | 5.8% | 15.7%
MBB |  |  |  |
| 2.6% | 1.7% | 0.2% | 0.0%

Figure 7. Optimized geometries for the various methods. The percentage under
each geometry represents its performance relative to the best performing
geometry in the row as measured by the objective function. The 0.0% geometry
is the best performing in the given row and the remainder are computed as the
respective difference as compared to the best performing objective value.

The optimized geometries using the various methods are seen in Figure 7. In
all cases, the neural parameterization results in more regular geometries. The
pixel-based methods could be regularized by augmenting the objective function
with a penalization term although this would require manual tuning of a
penalization parameter and may have an impact on the convergence. Simplistic
structures could also be obtained for the pixel-based cases by controlling the
filter radius and mesh resolution although this would prevent fine-scale
structure. The NN-LS method instead allows for fine-scale features resolved in
the final layer of the network but still produces performant regular
geometries as the multi-scale influence of the parameters in the neural
parameterization encourages globally performant structures to emerge. The use
of the NN is also seen to suppress numerical artifacts observed in the pixel-
LS solutions. For the MBB problem, the compliance measurements for all methods
fall within a few percent of each other. For this problem, the optimized
solutions roughly share the location of the major members and are only slight
variations from the optimal solution [Rozvany1998].

For the heat conduction problem, however, the NN-LS method has the best
performance by a fair amount followed by the NN-SIMP baseline. The regular
geometries produced by the NNs in these cases outperform their pixel
counterparts by a significant amount in both the LS and SIMP cases. The U-Net
here seems able to find better minima because of its ability to focus on
larger scale structures while the pixel-based parameterization makes
improvements locally by adding finer scale branches to the geometry. Without
handcrafting an initialisation close to the optimal solution, the addition of
the neural network avoids sub-optimal branching structures and instead leads
to the emergence of lamellar needles, the salient feature of optimal
geometries for this problem [Yan2018].

####  7.1.3. Convergence Plots

The convergence of the methods is plotted for the heat conduction and MBB
problems in Figure 8. Since the use of the word iteration in the context of
optimization is somewhat ambiguous, we plot the compliance against the number
of objective function value calls for the process. For the ADAM and MMA
optimizers, the ratio of function and gradient calls is 1:1. To compare the
SIMP and LS methods against each other, we bias the SIMP method using the
optimized structures converted LS once again. The bias is computed as the
difference between the SIMP method’s final objective value using the
interpolated material and the SIMP method’s final objective value using the
converted LS. This bias is applied to all of the series data in Figures 8 for
the SIMP methods.

(a) Heat conduction

(b) MBB

Figure 8. Convergence plot comparison for the benchmark problems. The
objective value at each objective function call is plotted.

For the heat conduction problem in Figure 8(a), the pixel-SIMP method
converges the fastest, albeit to an inferior solution. The rest of the methods
converge at similar rates. For the MBB problem in Figure 8(b), the SIMP
methods converge much faster than the LS methods although the NN-SIMP method
jumps out of the minimum at later iterations and stabilizes later on.

###  7.2. Interface Coupled Multiphysics Problem

The design of a support in a fluid-structure problem is optimized in this
section with the setup in [Jenkins2016]. The structure is optimized under a
forcing term from the fluid integrated on the evolving interface. The goal is
to demonstrate the generality of the method which is shown here by its
capacity to solve a multiphysics problem with interface coupling. Density
methods do not extend naturally to handle such problems and are faced with
difficulty in obtaining accurate coupling between the fluid and structure
since the representation of the interface is spread across cells in the
vicinity of the boundary.

$\Omega_{\mathrm{out}}$
---
$\Omega_{D}$
---
$\Omega_{D}$
---
$\Omega_{D}$
---
$\Omega_{\mathrm{in}}$
---
$\bm{u}_{in}$
---
$\Gamma_{N}$
---
Figure 9. The beam support problem setup. An inlet velocity is prescribed on
the left side, a homogenous Dirichlet condition on the top and bottom walls
and a homogenous Neumann condition on the right wall $\Omega_{N}$. The design
region supporting the beam in light grey can be either fluid or solid. The
solid domain $\Omega_{\mathrm{in}}$ is composed of the dark grey beam and the
solid part of the design region. The fluid region $\Omega_{\mathrm{out}}$ is
composed of the remainder of the channel, including the non-solid part of the
design region.

The problem setup is seen in Figure 9. For this problem, we set the fluid
parameters as $\mu=1m^{2}s^{-1}$ and $\alpha_{u}=2.5\mu/0.01^{2}$ and use a
parabolic velocity profile on the inlet with an average velocity of
$0.01ms^{-1}$. For the structural parameters, we set $E=1Pa$ and $\nu=0.3$ and
use a $0.45$ volume fraction. We use the same network as in the poisson
problem but set $w=(12,12,24,46,96,96)$ and $l=(24,24,48,96,192,192)$. The
optimized geometry for problem is seen in Figure 10.

0.0
---
0.5
---
0.1
---
0.2
---
0.3
---
0.4
---
$p_{h}$[Pa]
---
Figure 10. Optimized geometry for the beam support problem. Streamlines are
plotted to represent the velocity field and pressure is indicated in the
legend.

##  8\. Conclusion

Here we propose a neural LS method as a means of coupling convolutional NNs
and the unfitted LS method. The neural parameterization provides more regular
geometries and similar, in some cases better, performance compared to the
well-known SIMP method of TO and the pixel-based counterpart LS method. The
neural parameterization learns features at multiple scales during the
optimization allowing for the emergence of regular structures whilst
maintaining performance. The method takes longer to converge than the SIMP
method and may not be suitable for simple problems such as those in linear
elasticity. In this case, the network can be used in conjunction with the SIMP
method. In contrast to SIMP, however, the method extends naturally to
interface multiphysics problems and is more suitable in the general case. It
should be noted that in our approach, the expressivity of geometries for the
neural parameterization is the same as pixel counterpart methods since we use
the same space $V_{h}^{1}$ in all cases to describe the geometry. A direction
for development in the future would be to relax this constraint and also allow
for the NN to control the space in which it expresses itself, potentially
allowing for more efficient expressions of better topologies. Further
investigation into architecture design and optimization strategy may also lead
to a greater and more general improvement in performance.

##  9\. Replication of Results

All the codes being used in this paper are distributed as open-source
software. For reproducibility purposes, the implementation of the proposed TO
methodology as well as the drivers being used to compare the different
numerical methods are publicly available in the following repository:
https://github.com/ConnorMallon/NLSTO

## Acknowledgements

This research was partially funded by the Australian Government through the
Australian Research Council (project number DP220103160).

## Conflict of Interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

  *[TO]: topology optimization
  *[LS]: level set
  *[FE]: finite element
  *[FCM]: finite cell method
  *[DOF]: degrees of freedom
  *[FEs]: finite element
  *[NN]: Neural network
  *[NNs]: Neural network
  *[CutFEM]: cut finite element method
  *[AgFEM]: aggregated finite element method