Another source of mesh dependence in topology optimization.
LLNL]Miguel A. Salazar de Troyacorrespondingauthor
LLNL]Geoffrey M. Oxberry
LLNL]Cosmin G. Petra
LLNL]Daniel A Tortorelli
[LLNL]Lawrence Livermore National Laboratory, Livemore, CA, USA, 94550
[correspondingauthor]Corresponding author.

The topology optimization community has regularly employed nonlinear programming (NLP) algorithms  from the operations research community.
However, these algorithms are implemented in the real vector space $\mathbb{R}^n$ instead of the proper function space where the design variable resides.
In this article, we show how the volume fraction variable discretization on non-uniform meshes affects the convergence of $\mathbb{R}^n$ based NLP algorithms.
We do so by first summarizing the functional analysis tools necessary to understand why convergence is affected by the mesh.
Namely, the distinction between derivative and gradient definitions and the role of the mesh-dependent inner product within the NLP algorithm.
These tools are then used to make the Globally Convergent Method of Moving Asymptotes (GCMMA), a popular NLP algorithm in the topology optimization community, converge in a mesh independent fashion when starting from the same initial design.
We then benchmark our algorithm with different examples of topology optimization problems. Namely, ill-conditioned and large scale problems and designs problems with adaptive mesh refinement.

§ INTRODUCTION

Topology optimization finds the optimal distribution of material in a given design domain $D$ to minimize a cost function and satisfy constraint function inequalities.
In the classic element-wise uniform density based method, which in this paper we refer to as element-wise uniform volume fraction method[We opt to use volume fraction over density to avoid confusion with the physical quantity density, which is often used in topology optimization applications such as elastodynamics, fluids, etc.], the optimization algorithm places material in individual elements of a background mesh to define the geometry of the optimal design.
Regions devoid of material are meaningless, hence the motivation to coarsen the mesh in these regions.
By the same token, regions which contain material require a higher mesh resolution.
This cost saving strategy that distributes elements with different sizes within the mesh is known as the Adaptive Mesh Refinement (AMR).

If the elements have different sizes, as in AMR, it is intuitively wrong to think that all design variables have the same contribution to the design.
In particular, we explain why it is necessary to accommodate the element size when calculating the inner products involving the design variables within the NLP algorithms.
However, this is not done within most NLP algorithms used in the topology optimization community as they assume the design is a vector
in the real vector space $\mathbb{R}^n$, i.e., simply a vector of length equal to the number of elements in the mesh, $n$, cf. IPOPT <cit.>, SNOPT <cit.>,  MMA <cit.>, FMINCON <cit.> and Optimality Criteria <cit.>.
On the other hand, the NLP libraries Optizelle <cit.>, Moola <cit.>, ROL <cit.> and TAO <cit.> contain to various extents the capacity of treating design fields as elements of their underlying function spaces.
Related work by <cit.> compares mesh-independent and dependent versions of the steepest descent algorithms for an unconstrained problem and estimates their rates of convergence.
It is also necessary to mention that in the PDE-constrained optimization community, NLP algorithms are inherently mesh-independent as they are implemented in the corresponding function space.
For instance, <cit.> implements an infinite-dimensional (inf-dim) primal-dual interior-point method with a Newton solver,
<cit.> implements an inexact sequential quadratic programming method with an adaptive multilevel mesh refinement scheme and
<cit.> solves a phase field based topology optimization with a projected gradient method for cases where the cost function is only differentiable in $L^{\infty}$.

We choose the $L^2$ function space to represent the set of possible designs on the design domain $D$.
It is equipped with an inner product and discretized to be piecewise uniform[We use uniform to describe functions that do not change in space and constant to describe functions that do not change in time.] over the finite elements.
Other infinite-dimensional function spaces choices, e.g $H^1$, are possible and should be addressed in the future.
The inconsistency of using a design field in $L^2$ with an NLP algorithm formulated in $\mathbb{R}^n$ is generally not a problem because most of the topology optimization studies use uniform meshes.
However, when using meshes with different element sizes, as is the case in AMR, the $\mathbb{R}^n$ viewpoint yields mesh-dependent designs, whereas the $L^2$ approach does not.
An immediate corollary is that restriction by filtration alone does not ensure mesh-independent designs, as is commonly accepted in the topology optimization community.

This work is laid out as follows: Section <ref> presents the mathematical tools we need to implement mesh-independent NLP algorithms.
We use these tools in Section <ref> to make one of the most popular NLP algorithms in the topology optimization community, the Globally Convergent Method of Moving Asymptotes (GCMMA/MMA) algorithm, mesh-independent.
In Section <ref>, we first validate the NLP algorithm by solving three common problems in topology optimization with contrived meshes specifically built to increase the ill-conditioning of the optimization problem.
We then apply the algorithm to a three dimensional problem with a uniform mesh and different levels of refinement to showcase the importance of our algorithm for large scale problems.
Finally, we solve two design problems with different physics and AMR applied during the optimization.
Section <ref> briefly summarizes our findings and presents conclusions.

The function spaces concepts used in this article require a rigorous mathematical discussion to be absolutely precise.
Namely, it is necessary to show the differentiability of the convex approximation within the GCMMA to ensure that applying the Newton's method is mathematically sound.
These details would quickly obscure our main focus.
Our intent is to convey the differences between the $L^2$ and $\mathbb{R}^n$ NLP algorithms as simply as possible and to demonstrate their differences.
We therefore opt to take a more pragmatic approach in our discussions at the expense of glossing over important mathematical details.

§ MATHEMATICAL PRELIMINARIES

A topology optimization algorithm converges in a mesh-independent fashion by treating the design as a field, here a field in the $L^2$ space, using concepts from functional analysis.
The design field is then discretized in a consistent manner using the finite element basis, i.e. resulting in the widely used element-wise uniform volume fraction field[
In this work, for simplicity, we focus on the most popular topology optimization approach, where the design is defined
by an element-wise uniform material volume fraction in the Hilbert space $L^2$, but it could be extended to other parametrizations
such as an $H^1$ nodal-based material volume fraction.].
Notably, the norms in the NLP algorithm that check for convergence are discretized in this finite element space.

To illustrate the proper discretization, consider the unconstrained minimization problem

\begin{equation}
\begin{aligned}
\begin{split}
& \underset{\nu \in V}{\text{min}}
& & \theta(\nu)\,, \\
\end{split}
\label{eq:orig_problem}
\end{aligned}
\end{equation}

with the functional

\begin{align}
\theta & :V\rightarrow \mathbb{R}\,,
\end{align}

where $\nu: V \rightarrow \mathbb{R}$ is our volume fraction design field that belongs to $V$, a Hilbert space on domain $D$, equipped with an inner product $(\cdot, \cdot)_V$, which induces the primal norm $\norm{\cdot}_V$.
For our topology optimization, $V=L^2(D)$, which is equipped with the norm

\begin{align}
\norm{\nu}_{L^2} = \sqrt{(\nu, \nu)_{L^2}} =  \left( \int_{D} \nu^2 ~ dV  \right)^{1/2}\,.
\end{align}

The space of all bounded linear functionals that map $V$ to $\mathbb{R}$ is the dual space $V^*$, which is a subset of $\mathscr{L}(V, \mathbb{R})$, i.e. the space of linear operators from $V$ to $\mathbb{R}$, i.e. $V^* \subset \mathscr{L}(V, \mathbb{R})$.
Both the primal $\norm{\cdot}_V$ and dual $\norm{\cdot}_{V^*}$ norms can be used to check for convergence in NLP algorithms.

To formulate NLP algorithms on the function space $V$, we need the Riesz map from the Riesz representation theorem:
Let $V$ be a Hilbert space with inner product $(\cdot, \cdot)_V$ and dual space $V^*$.
For every $\varphi \in V^*$ there is a unique element $u\in V$ such that $\varphi (v) = \left( u, v \right)_V$ for all $v\in V$.
This one-to-one map is the Riesz map $\Phi:V \rightarrow V^*$ defined such that $\Phi(u) = \varphi$; it is an isometry between $V$ and $V^*$.

The discretization of the primal and dual spaces follow from <cit.>.
We approximate the volume fraction field $\nu \in V$ with $\nu_h \in V_h$, where $V_h$ is the span of basis functions $\mathscr{P} = \left\{\phi_1, ..., \phi_n \right\}$, $\phi_i \in \rojo{V_h}$ and $n$ is the dimension of $V_h$.
Our approximation now reads

\begin{align}
\nu(\mathbf{x}) \approx \nu_h(\mathbf{x}) & = \sumonn{i} \nu^i \phi_i(\mathbf{x}) = \boldsymbol{\nu}^T \boldsymbol{\phi}(\mathbf{x})\,.
\label{eq:l2_disc}
\end{align}

We similarly approximate $\iota\in V$ with $\iota_h \in V_h$ so that the inner product definition yields

\begin{equation}
\begin{aligned}
(\nu_h,\iota_h)_{V_h} & = \int_{D} (\boldsymbol{\nu}^T \boldsymbol{\phi}) (\boldsymbol{\iota}^T \boldsymbol{\phi}) dV \\
& = \boldsymbol{\nu}^T \int_{D} \boldsymbol{\phi} \boldsymbol{\phi}^T ~dV ~ \boldsymbol{\iota}  \\
& = \boldsymbol{\nu}^T \mathbf{M}  \boldsymbol{\iota} \,,\label{eq:disc_inner}
\end{aligned}
\end{equation}

\begin{align}
\mathbf{M} & = \int_{D} \boldsymbol{\phi} \boldsymbol{\phi}^T ~dV
\label{eq:mass_matrix}
\end{align}

is the mass matrix that reflects the mesh discretization.
By construction, $\mathbf{M}$ is symmetric and invertible.

The discretized design field $\nu_h$ is in the Hilbert space $V_h=\left(\mathbb{R}^n, (\cdot,\cdot)_{\mathbf{M}}\right)$, i.e. it is a vector in $\mathbb{R}^n$ of dimension $n$ with an $\mathbf{M}$ inner product.
This inner product induces the norm $\norm{\nu_h}_{V_h} = \norm{\boldsymbol{\nu}}_{\mathbf{M}}= (\boldsymbol{\nu}^T \mathbf{M} \boldsymbol{\nu})^{1/2}$.

Clearly the $L^2$ norm $\norm{\nu_h}_{L^2} = (\boldsymbol{\nu}^T \mathbf{M} \boldsymbol \nu)^{1/2}$ differs from the $\mathbb{R}^n$ norm $\norm{\nu_h}_{\mathbb{R}^n} = (\boldsymbol{\nu}^T  \boldsymbol \nu)^{1/2}$.
In topology optimization, $\nu_h$ is usually discretized via piecewise uniform functions over the individual elements so $\mathbf{M} = diag\left( |\Omega_1|, ..., |\Omega_n|\right)$  where $|\Omega_e|$ is the volume of the element $\Omega_e$.
So if the mesh is uniform, $\mathbf{M} = |\Omega_e| \mathbf{I}$ and hence $\norm{\nu_h}_{L^2} = \sqrt{|\Omega_e|} \norm{\nu_h}_{\mathbb{R}^n}$.

The basis $\mathscr{P}=\left\{\phi_1, ..., \phi_n \right\}$ induces a unique dual basis $\mathscr{P}^* = \left\{\phi^{*1}, ..., \phi^{*n} \right\}$ for $V_h^*$ defined such that $\phi^{*i} \in \rojo{V_h^*}$ and $\phi^{*i}(\phi_j)= \delta_{j}^i ~\forall ~i,j=1,...,n$.
This $\mathscr{P}^*$ basis is used to discretize any $F\in V^*$ as $F_h \in V_h^*$ such that for all $\iota_h \in V_h$

\begin{align}
F_h(\iota_h) = \sum^{n}_{i=1} F_i \phi^{*i} (\iota_h)\,,
\end{align}

where $F_i=F(\phi_i)$ for $i=1,...,n$, i.e. the vector components $F_i$ are interpolated from $F\in V^*$.
In this way, $F_h(\iota_h)$ is computed as

\begin{equation}
\begin{aligned}
F_h(\iota_h) & = \sumonn{i} F_i \phi^{*i} \left( \sumonn{j} \iota^j \phi_j \right)  \,, \\
& = \sumonn{i} F_i   \sumonn{j} \iota^j \phi^{*i}(\phi_j)   \,,            \\
& = \sumonn{i} F_i   \iota^i    \,,                                        \\
& =  \mathbf{F}^T   \boldsymbol{\iota}   \,,\label{eq:dual_norm_disc}
\end{aligned}
\end{equation}

where we used the orthonormal property between the bases $\mathscr{P}$ and $\mathscr{P}^*$ and the linearity of $\phi^{*i}$.
From the Riesz representation theorem, there exists $\nu_h\in V_h$ such that $\Phi(\nu_h) = F_h$ or $\Phi^{-1}(F_h) = \nu_h$, where $\Phi^{-1}: V_h^* \rightarrow V_h$.

\begin{equation}
\begin{aligned}
F_h(\iota_h) & = (\overbrace{\Phi^{-1}(F_h)}^{\nu_h}, \iota_h)_{V} \,,                           \\
& = (\nu_h, \iota_h)_{V} \,,                                                        \\
& = \boldsymbol{\nu}^T \mathbf{M}  \boldsymbol{\iota} \,.\label{eq:dual_norm_disc2}
\end{aligned}
\end{equation}

From Equations (<ref>) and (<ref>) we can see that

\begin{align}
\mathbf{F}^T   \boldsymbol{\iota}  =  \boldsymbol{\nu}^T \mathbf{M}  \boldsymbol{\iota}\,.
\end{align}

\begin{align}
\mathbf{F} = \mathbf{M} \boldsymbol{\nu} \,,
\end{align}

and the discrete Riesz map and its inverse are defined such that

\begin{align}
\Phi_h(\boldsymbol{\nu})=\mathbf{M}\boldsymbol{\nu}
\end{align}

\begin{align}
\Phi_h^{-1}(\mathbf{F}) & = \mathbf{M}^{-1} \mathbf{F} \,.
\label{eq:rieszmap}
\end{align}

Recalling that the Riesz map is an isometry between the spaces $V_h$ and $V_h^*$, we can now define and calculate the norm of an object $F_h \in V_h^*$ as

\begin{equation}
\begin{aligned}
\norm{F_h}_{V^*_h} & = \norm{\Phi^{-1} (F_h)}_{V_h} \,,                                               \\
& = \norm{\mathbf{M}^{-1} \mathbf{F}}_{\mathbf{M}} \,,                             \\
& = \sqrt{\mathbf{F}^T \mathbf{M}^{-1} \mathbf{M} \mathbf{M}^{-1} \mathbf{F} } \,, \\
& = \norm{\mathbf{F}}_{\mathbf{M}^{-1}} \,,
\end{aligned}
\end{equation}

where we used the definition of the discrete inner product in Equation (<ref>).

In this work, we use the Fréchet derivative $D\theta(\nu) \in V^*$ of the function $\theta : V \rightarrow \mathbb{R}$ at $\nu$.
If it exists, this derivative is defined such that[The “little-$o$ notation” $o(\norm{h}_V)$ for a functional $q: V \rightarrow \mathbb{R}$ means
$\lim_{\norm{h}_V \to 0} \frac{q(h)}{ \norm{h}_V} =0$]

\begin{align}
\theta(\nu + h) - \theta(\nu) - D\theta(\nu)[h] = o( \norm{h}_V)
\label{eq:frechet_deriv}
\end{align}

for all $h\in V$.
By definition, $D\theta(\nu)\in V^*$, and hence the Riesz representation theorem tells us there is an object in $V$ that we will denote $\nabla \theta (\nu) \in V$, i.e. the gradient of $\theta$ at $\nu$ such that
$D\theta(\nu)[h]=(\nabla \theta(\nu), h)_V$ for all $h\in V$.
Using the Riesz map, $\Phi(\nabla \theta(\nu)) = D\theta(\nu)$ and because the Riesz map depends on the inner product $(\cdot, \cdot)_V$, so does $\nabla \theta(\nu) \in V$.
This inner product dependence is crucial in our NLP algorithm as the inner product depends on the mesh discretization, notably from (<ref>) we have

\begin{align}
\boldsymbol{\nabla \theta} = \mathbf{M}^{-1} \boldsymbol{D \theta}
\label{eq:rieszmap_grad}
\end{align}

where $\boldsymbol{\nabla \theta}$ and $\boldsymbol{D \theta}$ are the discrete counterparts of $\nabla \theta(\nu)$ and $D\theta(\nu)$.

We are now in position to show how these functional analysis concepts apply to NLP algorithms in the inf-dim space.
We start by examining the most basic NLP algorithm for the solution of the simple unconstrained minimization problem of (<ref>).
i.e, the steepest descent algorithm, for which the iterate $\nu^{(k)}$ is updated as

\begin{align}
\nu^{(k+1)} = \nu^{(k)} - \gamma \nabla \theta(\nu^{(k)}) \,.
\label{eq:steepestdesc}
\end{align}

where $\gamma \geq 0$ is the step length.
The discretized Equation (<ref>) becomes

\begin{equation}
\begin{aligned}
\boldsymbol{\nu}^{(k+1)} & = \boldsymbol{\nu}^{(k)} - \gamma \boldsymbol{\nabla \theta} \,.            \\
& = \boldsymbol{\nu}^{(k)} - \gamma \mathbf{M}^{-1} \boldsymbol{D \theta} \,.
\label{eq:steepestdescdisc}
\end{aligned}
\end{equation}

When calculating the norm to check for convergence, we use $\norm{\nabla \theta (\nu)}_V$, which upon discretization is $\norm{\boldsymbol{\nabla \theta}}_{\mathbf{M}}$.

It seems intuitive that $\nu$ and $\nabla \theta(\nu)$ must be in the same function space since they are added together.

This motivates us to use the gradient $\nabla \theta(\nu)$ and not the derivative $D\theta(\nu)$ in (<ref>), which is contrary to most topology optimization algorithms.
For a uniform mesh $\mathbf{M}=|\Omega_e|\mathbf{I}$ so $\boldsymbol{\nabla \theta} = \frac{1}{|\Omega_e|} \mathbf{I} \boldsymbol{D\theta} $
and hence $\boldsymbol{\nabla \theta}$
and $\boldsymbol{D\theta}$ are
parallel and there is no difference in the search direction.
However, the number of iterations to convergence will be different due to the difference in the inner product.

The second NLP algorithm uses Newton's method, wherein we iterate to find $\nu$ such that

\begin{align}
D\theta(\nu)[\delta \nu]= 0~ \forall \delta \nu \in V.
\end{align}

To do so, we linearize around $\nu^{(k)}$ and solve for the update $\Delta \nu^{(k)}$ via

\begin{align}
D^2\theta(\nu^{(k)})[\Delta \nu^{(k)},\delta \nu] =-D\theta(\nu^{(k)})[\delta \nu]~~ \forall \delta \nu \in V \,,
\label{eq:newtonstep}
\end{align}

where $D^2\theta(\nu^{(k)})[ \cdot, \cdot]: V \rightarrow  \mathscr{L} (V \times V, \mathbb{R})$ is the Hessian, i.e. second derivative of $\theta$ at $\nu^{(k)}$; it is a bilinear map from $(V \times V)$ to $\mathbb{R}$.
The difference here is that we need to supply the NLP algorithm with the derivative $D\theta(\nu)$ (and the Hessian $D^2\theta(\nu)$) and not the gradient $\nabla \theta(\nu)$ as in the steepest descent algorithm.
Upon discretization, Equation (<ref>) becomes

\begin{align}
\boldsymbol{D}^2\boldsymbol{\theta} \boldsymbol{\Delta \nu}^{(k)}  =-\boldsymbol{D\theta}
\end{align}

When calculating the norm to check for convergence, we use $\norm{D \theta (\nu)}_{V^*}$, whose discretization is $\norm{\boldsymbol{D \theta}}_{\mathbf{M}^{-1}}$

Inspired by <cit.>, we showcase the difference between $\nabla \theta(\nu)$ and $D\theta(\nu)$ with the following one dimensional convex unconstrained optimization problem

\begin{align}
\underset{\nu\in V}{\text{min}} ~ \theta(\nu) = \frac{1}{2}( \nu, c \, \nu )_V - ( \nu, b )_V \,,
\end{align}

where $b$ and $c$ are given functions on $V=L^2(D)$ with $D=[1, 10]$, $c(x) = \sinfunc$ and $b(x) = x$.
The solution is trivially calculated by the stationary condition

\begin{equation}
\begin{aligned}
D\theta(\nu) [v] = 0 & =  ( \nabla \theta, v )_V  \,,  \\
& =      ( c \, \nu - b,v )_V \,.
\end{aligned}
\label{eq:derivexample}
\end{equation}

which must hold for all $\nu \in V$ and hence $\nu(x)=\frac{b(x)}{c(x)}$.
We proceed to discretize the function $\theta$ with piecewise uniform elements, resulting in the expression

\begin{align}
\theta(\nubold) =  \frac{1}{2} \nubold^T \mathbf{H} \nubold - \nubold^T \mathbf{M b}  \,,
\end{align}

whose derivative is

\begin{align}
\boldsymbol{D \theta}=  \gradient \,,
\label{eq:discderiv}
\end{align}

where the Hessian matrix $\mathbf{H}$ is calculated as

\begin{align}
\mathbf{H} = \int_D c(x) \boldsymbol \phi(x)  \boldsymbol \phi(x)^T ~dV \,,
\end{align}

using one point quadrature per element.
The vectors $\boldsymbol \nu$ and $\bf b$ represent the values of the functions $\nu(x)$ and $b(x)$ at the quadrature points.
Note that Equation (<ref>) is the discretization of $D\theta(\nu)$ in Equation (<ref>).
Applying the steepest descent algorithm in $\mathbb{R}^n$ yields

\begin{equation}
\begin{aligned}
\nubold^{(k+1)} & = \nubold^{(k)} - \gamma \boldsymbol{D \theta} \,, \\
& = \nubold^{(k)} - \gamma(\gradient)\,.
\label{eq:steepestexample}
\end{aligned}
\end{equation}

The optimal step size $\gamma$ is calculated with the closed-form expression

\begin{align}
\gamma = \stepsize \,,
\end{align}

which we use in Equation (<ref>) to obtain the fixed point iteration

\begin{align}
\nubold^{(k+1)} = \nubold^{(k)} - \stepsize (\gradient) \,.
\label{eq:fixedsteepest}
\end{align}

On the other hand, in the $L^2$ reformulation,
we replace the descent direction $\boldsymbol{D \theta}$ with the gradient $\boldsymbol{\nabla \theta} = \mathbf{M}^{-1} \boldsymbol{D\theta}(\nubold) =  \newgradient$ in Equation (<ref>).a and calculate the optimal step size

\begin{align}
\gamma = \newstepsize = 1
\end{align}

to obtain the fixed point iteration

\begin{align}
\nubold^{(k+1)} = \nubold^{(k)} - (\newgradient)\,.
\label{eq:newfixedsteepest}
\end{align}

We run both fixed point iteration Equations (<ref>) and (<ref>) starting from $\boldsymbol \nu^{(0)}=\boldsymbol 0$.
Convergence is declared when the error $e \leq 10^{-7}$, where $e=\norm{\boldsymbol \nu - \mathbf{b}\oslash \mathbf{c}}$ for $\mathbb{R}^n$ (the operator $\oslash$ is the Hadamard division) and $e=\norm{\nu - b / c}_{L^2}$ for $L^2$.

Table <ref> denotes the iteration history of both methods with $n$ design variables over a one-dimensional mesh with nodes at positions $x_r=10^{y_r}$ where $y_r=\{\frac{r}{n+1}\}_{r=0}^{n+1}$.
The convergence in the $\mathbb{R}^n$ NLP algorithm deteriorates with the number of elements, whereas the number of iterations in the $L^2$ NLP algorithm remains nearly constant.
The slower convergence of the $\mathbb{R}^n$ NLP algorithm is attributed to the fact that we are adding members from different spaces ($\boldsymbol\nu \in V_h$ and $\boldsymbol{D\theta}\in V_h^*$).

$10^1$     $10^3$     $10^5$

$\mathbb{R}^n$     195        270        302

$L^2$              52         56         56

Iteration count for the discrete steepest descent of algorithms Equations (<ref>) and (<ref>).

§ GCMMA IN FUNCTION SPACE

We are now motivated to formulate the first-order GCMMA algorithm <cit.> in the $L^2$ space.
GCMMA <cit.> and its non globally-convergent version MMA <cit.> are widely used NLP algorithms in the topology optimization community.
Following the implementation given in <cit.>, we highlight here the necessary changes to make the GCMMA algorithm converge in a mesh-independent fashion.
To begin, we consider the optimization problem

\begin{equation}
\begin{aligned}
\label{eq:minproblem}
\begin{split}
& \underset{\nu \in V}{\text{min}}
& & \theta_0(\nu) \,,\\
& \text{s.t.}
& & \theta_i(\nu) \leq 0, \; i = 1, \ldots, m\,, \\
&&& \nu_{\text{min}} \leq \nu \leq \nu_{\text{max}} \; \text{a.e}\,.
\end{split}
\end{aligned}
\end{equation}

First, the artificial optimization variables $\mathbf{y} = (y_1,...,y_m)$ are added to ensure feasibility and $z$ is added make certain subclasses of problems, like least squares or minmax problems, easier to formulate, i.e.

\begin{equation}
\begin{aligned}
& \underset{
\subalign{
\nu            & \in V                                                                                \\
\boldsymbol{y} & \in \mathbb{R}^m                                                                                \\
z              & \in \mathbb{R}
&                  & \theta_0(\nu) + a_0z + \sum\limits_{i=1}^m \left( c_i y_i + \frac{1}{2}d_i y_i^2 \right) \,,                                                                     \\
& \text{s.t.}
&                  & \theta_i(\nu) - a_iz - y_i\leq 0, \; i = 1, \ldots, m\,,                                                                                \\
&                  &                                                                                & \nu_{\text{min}} \leq \nu  \leq \nu_{\text{max}}\; \text{a.e.}\,, \\
&                  &                                                                                & \boldsymbol{y} \geq 0 \,,                                         \\
&                  &                                                                                & z\geq 0 \,.
\end{aligned}
\label{eq:mmaproblem}
\end{equation}

where $a_0, a_i, c_i$ and $d_i$ are real numbers which satisfy $a_0 > 0, a_i \geq 0, c_i \geq 0, d_i \geq 0$ and $c_i + d_i > 0$ for all $i$, and also $a_i c_i > a_0$ for all $i$ <cit.>.
Note that we recover the original NLP algorithm for $z=0$ and $\boldsymbol{y}=0$.

For each optimization iteration $k$, we solve the following convex approximate subproblem based on Equation (<ref>), the cost and constraint functions and
their derivatives and the values at the current iterate $(\nu^{(k)},\boldsymbol{y}^{(k)},z^{(k)})$.
Ultimately, we iterate by solving

\begin{equation}
\begin{aligned}
(\nu^{(k+1)}, \boldsymbol{y}^{(k+1)}, & z^{(k+1)}) =
\underset{
\subalign{
\nu                                   & \in V                                                                                \\
\boldsymbol{y}                        & \in \mathbb{R}^m                                                                                \\
z                                     & \in \mathbb{R}
}{\text{arg min}}
&                  & \tilde{\theta}_0(\nu) + a_0z + \sum\limits_{i=1}^m \left( c_i y_i + \frac{1}{2}d_i y_i^2 \right) \,,                                               \\
& \text{s.t.}
&                  & \tilde{\theta}_i(\nu) - a_iz - y_i\leq 0, \; i = 1, \ldots, m. \,,                                                                                \\
&                  &                                                                                & \alpha \leq \nu \leq \beta\; \text{a.e} \,, \\
&                  &                                                                                & \boldsymbol{y} \geq 0 \,,                   \\
&                  &                                                                                & z\geq 0 \,.
\end{aligned}
\label{eq:mmasubproblem}
\end{equation}

where the newly introduced functions $\tilde{\theta}_0$ and  $\tilde{\theta}_i$ and bounds $\alpha$ and $\beta$ are defined momentarily.

In our formulation of the above subproblem, we replace the summations of the approximating functionals $\tilde{\theta}_i$ in <cit.> with integrals over the domain.

\begin{align}
\tilde{\theta}_i(\nu) & =\int_D \left(\frac{p_{i}}{U^{(k)} - \nu}+\frac{q_{i}}{\nu - L^{(k)}}\right)~dV + r_i , \; i = 0,1, \ldots, m \,,      \label{eq:convex_approx} \\
r_i                   & = \theta_i(\nu^{(k)}) - \int_D \left(\frac{p_{i}}{U^{(k)} - \nu^{(k)}}+\frac{q_{i}}{\nu^{(k)} - L^{(k)}}\right)~dV  \,,
\label{eq:mma_integrals}
\end{align}

\begin{align}
p_{i} & = (U^{(k)} - \nu^{(k)})^2\left(1.001\left(\nabla \theta_i(\nu^{(k)})\right)^+
+0.001\left(\nabla \theta_i(\nu^{(k)})\right)^- +\frac{\rho_i^{(k,j)}}{\nu_{\text{max}} - \nu_{\text{min}}}\right) \label{eq:pfunc} \,, \\
q_{i} & = (\nu^{(k)} - L^{(k)})^2\left(0.001\left(\nabla \theta_i(\nu^{(k)})\right)^+
+1.001\left(\nabla \theta_i(\nu^{(k)})\right)^- +\frac{\rho_i^{(k,j)}}{\nu_{\text{max}} - \nu_{\text{min}}}\right) \label{eq:qfunc} \,.
\end{align}

and $U^{(k)}$ and $L^{(k)}$ are the soon to be defined moving upper and lower asymptotes; they are all elements of $V$.
We emphasize here that the original GCMMA implementation does not make a distinction between gradients and derivatives when building the convex approximation in Equation (<ref>).
As shown in this paper, it is vital to use gradients $\nabla \theta_i(\nu^{(k)})$ for $i=0,1, \ldots,m$.
It is important to warn readers that the convex approximations in Equation (<ref>) are not Frechét differentiable where $\nu=0$ in $L^2$.
We do not allow this $\nu=0$ situation, but a more mathematically rigorous rederivation of the method to accept such cases should be considered in the future.
To ensure the subproblem is convex, we use the ramp like functions $\left(a\right)^+ = \text{max} (0, a)$, and $\left(a\right)^-= \text{max} (0 ,-a)$.

The bounds (now fields in $V$) $\alpha$ and $\beta$ are taken as

\begin{equation}
\begin{aligned}
\alpha & = \text{max} \{\nu_{\text{min}}, L^{(k)} + 0.1 (\nu^{(k)} - L^{(k)}), \nu^{(k)} - 0.5 (\nu_{\text{max}} -\nu_{\text{min}} )\} \\
\beta  & = \text{min} \{\nu_{\text{max}}, U^{(k)} - 0.1 (U^{(k)} - \nu^{(k)}), \nu^{(k)} + 0.5 (\nu_{\text{max}} -\nu_{\text{min}} )\}
\end{aligned}
\label{eq:alphabeta}
\end{equation}

The GCMMA differs from the MMA in its attempt to achieve global convergence by controlling the parameter $\rho_i^{(k,j)}$ (which in the MMA is a fixed small positive value, usually lower than $10^{-5}$) in Equations (<ref>) and (<ref>).
Here, the added superscript $j$ corresponds to the inner iteration within the GCMMA.
For the initial $j=0$ inner iteration, the solution $(\nu^{(k,0)},\boldsymbol{y}^{(k,0)},z^{(k,0)})$ of the Equation (<ref>) subproblem, whose details are explained later, is accepted if

\begin{align}
\tilde \theta_i(\nu^{(k,j)}) \geq \theta_i(\nu^{(k,j)})  ;~ ~ i=0,...,m  \,,
\label{eq:condition_check}
\end{align}

whereupon the outer iteration $k+1$ commences with the initial iterate $(\nu^{(k+1)},\boldsymbol{y}^{(k+1)},z^{(k+1)}) = (\nu^{(k,0)},\boldsymbol{y}^{(k,0)},z^{(k,0)})$.
Otherwise, the $j+1$ subproblem (<ref>) is solved with a more conservative convex approximation by replacing $\rho_i^{(k,j)}$ with $\rho_i^{(k,j+1)} > \rho_i^{(k,j)}$ and the Equation (<ref>) inequality is rexamined.
If Equation (<ref>) is satisfied, we begin the outer iteration $k+1$ with the initial iterate $(\nu^{(k+1,0)},\boldsymbol{y}^{(k+1,0)},z^{(k+1,0)}) = (\nu^{(k, j+1)}, \boldsymbol{y}^{(k, j+1)}, z^{(k, j+1)})$, otherwise, the inner $j+2$ subproblem is solved and so on, cf. Figure <ref>.
The termination criteria will be explained in detail later.

block = [rectangle, draw, fill=blue!20,
text width=3cm, text centered, rounded corners, minimum height=0.5cm]
input = [draw, ellipse, fill=green!20,
text width=3cm, text centered, minimum height=0.8cm]
output = [draw, ellipse, fill=red!20,
text width=3cm, text centered, minimum height=0.8cm]
decision = [diamond, aspect= 3, draw, fill=orange!20, inner sep=-6pt,
text width=5cm, text centered, minimum height=0.5cm]
decision2 = [diamond, aspect= 3, draw, fill=orange!20, inner sep=-6pt,
text width=5cm, text centered, minimum height=0.4cm]
line = [draw, -latex]
[node distance = 1.2cm, auto]
[input] (nu_k) $\nu^{(0)}$;
[block, below of = nu_k] (obj_grad) Calculate $\theta_i(\nu^{(k)}), ~\nabla \theta_i(\nu^{(k)})$;
[block, below of = obj_grad] (param_gc) Obtain $\rho^{(k,0)}$;
[block, below of = param_gc, node distance = 1.2cm] (fapp) Build $\tilde \theta_i^{(k,j)}(\nu^{(k)})$;
[block, text width=5cm, below of = fapp] (nu_new) Solve for $\nu^{(k,j)}$ in subproblem Equation (<ref>);
[block, text width=5cm, below of = nu_new] (new_f) Calculate $\theta_i (\nu^{(k,j)})

[decision, below of = new_f, node distance = 2.5cm] (f_condition)

\begin{align*}
\text{If} ~ \tilde \theta_i(\nu^{(k,j)})
& \geq \theta_i(\nu^{(k,j)}) \\
\forall i=0,...,m &
\end{align*}
[output, below of = f_condition, node distance = 3cm] (new_opti) $\nu^{(k+1, 0)} = \nu^{(k,j)}$;
[block, right of = f_condition, node distance = 6cm] (repeat) Obtain $\rho^{(k,j+1)}$;
[line] (repeat) |- node[near start]$j=j+1$ (fapp);
[line] (f_condition) – nodeNo (repeat);
[line] (f_condition) – nodeYes (new_opti);
[line] (nu_k) – (obj_grad);
[line] (obj_grad) – (param_gc);
[line] (param_gc) – node $j=0$ (fapp);
[line] (fapp) – (nu_new);
[line] (nu_new) – (new_f);
[line] (new_f) – (f_condition);
[decision2, below of = new_opti, node distance = 1.5cm] (stopping) Termination criteria;
[line] (stopping) – nodeNo ($(stopping) + (-5,0)$) |- node[near start]$k=k+1$  (obj_grad);
[line] (new_opti) – (stopping);
[block, below of = stopping, node distance = 1.5cm] (final) Solution;
[line] (stopping) – nodeYes (final);

GCMMA algorithm.

The parameters $\rho_i^{(k,j)}$ are calculated following <cit.>, but with integrals over $D$ replacing summations.
For subproblem $(k,0)$

\begin{align}
\rho^{(k,0)}_i = \frac{0.1}{\hat V} \int_D \lvert \nabla \theta_i (\nu^{(k,0)}) \rvert \left( \nu_{\text{max}}  - \nu_{\text{min}}  \right) ~dV ~\text{for}\; i = 0,1, \ldots, m \,,
\label{eq:rho_mechanism}
\end{align}

where $\hat V$ is the volume (area) of $D$.
For the subsequent $(k,j+1)$ subproblems

\begin{equation}
\begin{aligned}
\rho^{(k, j+1)}_i & = \text{min} \left\{ 1.1 \left(\rho^{(k, j)}_i + \delta_i^{(k,j)} \right) , 10 \rho_i^{(k,j)} \right\} & \text{if} ~\delta_i^{(k,j)} > 0 \,,   \\
\rho^{(k, j+1)}_i & = \rho^{(k, j)}_i                                                                                & \text{if} ~\delta_i^{(k,j)} \leq 0\,,
\end{aligned}
\end{equation}

\begin{align}
\delta^{(k,j)}_i = \frac{\theta_i\left(\nu^{(k,j)}\right) - \tilde{\theta}_i\left(\nu^{(k,j)}\right)}{d\left(\nu^{(k,j)}\right)} \,.
\end{align}

\begin{align}
d (\nu) = \int_D \frac{\left( U^{(k)} - L^{(k)} \right) \left(\nu - \nu^{(k)} \right)^2}{\left(U^{(k)} - \nu\right) \left( \nu - L^{(k)} \right) \left( \nu_{\text{max}}  - \nu_{\text{min}} \right)} ~dV \,,
\label{eq:d_mechanism}
\end{align}

The moving asymptote fields $L \in V$ and $U \in V$ are updated via heuristic rules.
For iterations $k=1$ and $k=2$

\begin{equation}
\begin{aligned}
L^{(k)} = \nu^{(k)} - 0.5(\nu_{\text{max}} - \nu_{\text{min}}) \,, \\
U^{(k)} = \nu^{(k)} + 0.5(\nu_{\text{max}} - \nu_{\text{min}}) \,.
\end{aligned}
\label{eq:LUupdate1}
\end{equation}

For iterations $k \geq 3 $,

\begin{equation}
\begin{aligned}
L^{(k)} = \nu^{(k)} - \gamma^{(k)}(\nu^{(k-1)} - L^{(k-1)}) \,, \\
U^{(k)} = \nu^{(k)} + \gamma^{(k)}(U^{(k-1)} - \nu^{(k-1)}) \,.
\end{aligned}
\label{eq:asympt}
\end{equation}

The field $\gamma$ that appears in Equations (<ref>) is determined by the values of $\nu$ in the last three outer iterations.
When there is no oscillation in $\nu$, we reduce the convexity by pushing the asymptotes further apart by choosing a larger $\gamma$ to accelerate convergence.
Otherwise, we use a smaller value to move the asymptotes closer together.
Specifically, we assign

\begin{equation}
\begin{aligned}
\gamma^{(k)} =
\begin{cases}
0.7 & \text{if}~(\nu^{(k)} - \nu^{(k-1)})(\nu^{(k-1)} - \nu^{(k-2)}) < 0 \,,  \\
1.2 & \text{if} ~(\nu^{(k)} - \nu^{(k-1)})(\nu^{(k-1)} - \nu^{(k-2)}) > 0 \,, \\
1   & \text{if} ~(\nu^{(k)} - \nu^{(k-1)})(\nu^{(k-1)} - \nu^{(k-2)}) = 0\,,
\end{cases}
\end{aligned}
\label{eq:gammaupdate}
\end{equation}

subject to the inequalities

\begin{equation}
\begin{aligned}
L^{(k)} & \leq  \nu^{(k)} - 0.01(\nu_{\text{max}} - \nu_{\text{min}}) \,, \\
L^{(k)} & \geq  \nu^{(k)} - 10(\nu_{\text{max}} - \nu_{\text{min}}) \,,   \\
U^{(k)} & \geq  \nu^{(k)} + 0.01(\nu_{\text{max}} - \nu_{\text{min}}) \,, \\
U^{(k)} & \leq  \nu^{(k)} + 10(\nu_{\text{max}} - \nu_{\text{min}}) \,.
\end{aligned}
\label{eq:LUrestrict}
\end{equation}

We remark that the pointwise Equations (<ref>), (<ref>) - (<ref>) are discretized directly by using their element-wise counter parts, which is consistent with our $L^2$ element-wise piecewise uniform parameterization of $\nu$.
We do not worry here about the existence of these discretized counterparts as it is outside of the scope of this paper.

From here, one can solve the MMA subproblem following similar steps to those in the original article <cit.> with the exception of the calculation
of the norms in their corresponding function spaces.
We provide these details in <ref>, where we also summarize all the necessary changes to the original GCMMA.

§ NUMERICAL EXAMPLES

To illustrate the effectiveness of incorporating the $L^2$ function space approach in the GCMMA, we solve six topology optimization problems, five in two (2D) and one in three (3D) dimensions.
For the 2D cases, we use triangular elements, and hexahedral elements for 3D,
both with first order Lagrange basis functions to represent the displacement.
The volume fraction field discretization uses the typical topology optimization
approach with element-wise uniform basis functions.
All of the examples were solved using the finite element library Firedrake <cit.> <cit.> <cit.>,
which uses PETSc <cit.> <cit.> <cit.> as the backend for the linear algebra.
We use the direct solver MUMPS <cit.> in 2D and the PETSc GAMG preconditioner in 3D.
The 2D optimizations ran on a single 2.60 GHz Intel XeonE5-2670 processor.
We employed up to 36 processors for the 3D cases.
The modified GCMMA is an adaptation of a Python implementation of the original MMA algorithm from the GetDP finite element library <cit.>.
It was rewritten for better performance in parallel and to include an interface for the Firedrake-adjoint library <cit.>.
We use the MMA parameters $a_0 = 1$, $c_i = 10000$ and $a_i = d_i = 0$ for all $i \geq 1$.
All results are visualized with ParaView <cit.> and the graphs are plotted with Matplotlib <cit.>.
To launch all the simulations, we use Signac <cit.> and Signac-flow <cit.>.

§.§ Ill-conditioned meshes

In this subsection, we deliberately use meshes with highly refined regions which at first glance,
can be deemed as cherry-picked to validate our approach.
However, these meshes render ill-conditioned optimization problems and it is precisely this issue
that we wish to highlight and resolve.
We first solve three common topology optimization problems in linear elasticity.
The topology optimization problem is formulated as

\begin{equation}
\begin{aligned}
\underset{\nu \in \rojo{V}}{\text{min}}~ \theta_0(\nu)                                      & = \int_{D} \pi(\hat{\nu}, \mathbf{u})~ dV \,,                                 \\
\text{such that} ~ \mathbf{u} \in \rojo{W} \text{satisfies} ~ a(\nu; \mathbf{u},\mathbf{v}) & = L(\mathbf{v}) ~ \text{for all}~ \mathbf{v} \in \rojo{W}  \,,                \\
\theta_i(\nu)                                                                               & =  \int_{D} g_i(\hat{\nu}, \mathbf{u}) dV \leq 0               & i=1,2...m\,,
\end{aligned}
\label{eq:elastic_optimization}
\end{equation}

\begin{equation}
a(\nu; \mathbf{u},\mathbf{v}) =\int_{D} r(\hat\nu) \mathbb{C}[\nabla \mathbf{u}] \cdot \nabla \mathbf{v}  ~dV  \,,
\label{eq:weak_form_simp_mma}
\end{equation}

\begin{equation}
L(\mathbf{v}) =\int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v} ~da  \,.
\end{equation}

The function spaces used are

\begin{equation}
V = \{\nu \in L^2(D) ~|~ 0 \leq \nu \leq 1\}
\end{equation}

\begin{equation}
W = \{\mathbf{u} \in [H^1(D)]^3 ~|~ \mathbf{u}|_{\Gamma_D} = 0\}
\end{equation}

The domain boundary $\Gamma$ is comprised of three complementary regions: $\Gamma_D$, $\Gamma_N$ and
$\Gamma_F$ over which the Dirichlet, non-homogeneous Neumann and the homogeneous Neumann boundary conditions are applied.
The functionals $\theta_i:L^2 \to \mathbb{R}$, $i=0,1,...$ are assumed to be Fréchet differentiable.
The filtered volume fraction $\hat\nu$ in the above is obtained from the PDE-based filter <cit.> to generate a well-posed topology optimization problem

\begin{equation}
\begin{aligned}
-\kappa \nabla^2 \hat{\nu} + \hat\nu     & = \nu ~ & \text{in} ~  D      & \,, \\
\kappa \nabla \hat{\nu} \cdot \mathbf{n} & = 0~    & \text{on} ~  \Gamma &
\end{aligned}
\label{eq:pdefiltermma}
\end{equation}

where $\kappa$ determines the minimum length scale of the design such that a small (large)
$\kappa$ allows for fine (coarse) scale design fluctuations.
We solve Equation (<ref>) with a finite volume scheme to maintain the
volume fraction values between 0 and 1.
We use the SIMP penalization <cit.> to encourage 0-1 designs,
i.e. designs where 0 or $\nu=1$ almost everywhere.
As such,

\begin{align}
r(\hat\nu) = \epsilon_{\nu} + (1 - \epsilon_{\nu})\hat{\nu}^3 \,,
\label{eq:simp}
\end{align}

where $\epsilon_{\nu}=10^{-5}$ ensures the stiffness matrix in the finite element analysis is nonsingular.

Finally, $\mathbb{C}$ is the elasticity tensor corresponding to an isotropic material with Young modulus $E=1$ and Poisson ratio $\upnu=0.3$ and $\mathbf{t}$ is the applied traction on the surface $\Gamma_N$.

As usual, a reduced space approach is taken wherein we account for dependence of $\mathbf{u}$ on $\nu$, i.e. $\mathbf{u} \rightarrow \mathbf{u}(\nu)$ and the adjoint method is used to calculate the derivatives of the cost and constraint functions $\theta_i$.

The first problem we study is the proverbial structural compliance minimization subject to a maximum volume constraint $\hat V = 0.3 |D|$ i.e.

\begin{align}
\theta_0 & = \int_{\Gamma} \mathbf{t} \cdot \mathbf{u} ~da \label{eq:compliance} \,, \\
\theta_1 & = \int_{D} \hat \nu ~dV  - \hat V \label{eq:volume_comp}  \,.
\end{align}

The design domain $D$, cf. Figure <ref>, is subject to the traction
$\mathbf{t}=-1.0\mathbf{e}_2$ on $\Gamma_N$, the length scale parameter is $\kappa = 0.2$ and the initial design is a uniform field $\nu(\mathbf{x}) = 0.1$.
We perform four different optimizations corresponding to
uniform and nonuniform meshes with optimizations in $\mathbb{R}^n$ and $L^2$.
The uniform mesh contains 128,000 elements.
Our non-uniform mesh is illustrated in Figure <ref> and contains 80,577 elements.
It is important to have meshes that are sufficiently refined so the infinite-dimensional response $\mathbf{u}$ is well approximated.
We also include a highly refined arbitrary region on the top to clearly illustrate the deficiency of the NLP algorithm in $\mathbb{R}^n$.
Both meshes are included as beam_uniform.geo and beam_amr.geo in the files to reproduce the results.
In all cases, we run the optimization problems until the number of iterations reaches 200, although some designs converge sooner.

The optimized designs in Table <ref> show that the $\mathbb{R}^n$ NLP algorithm is mesh dependent as opposed to the $L^2$ NLP algorithm.
Figures <ref> - <ref> show the evolution of the cost and constraint functions and the convergence metric, cf. Equations (<ref>).
It is well known that the topology optimization problem is not convex and hence different initial designs might lead to different local minima.
The examples in Table <ref> start from the same initial design.
We therefore attribute the difference in the optimized designs to the mesh dependency of the NLP algorithm in $\mathbb{R}^n$.

[  spring/.style = decorate,decoration=zigzag,amplitude=6pt,segment length=4pt
, scale=0.1, every node/.style=inner sep=0pt]
(SW) at (0,0) ;
(SE) at (100,0) ;
(NE) at (100,40) ;
(NW) at (0,40) ;
[fill=gray!50!white] (0,0) – (100,0) – (100,40) – (0,40) – (0,0);
(D) at (50, 20) $D$;

[pattern=north west lines] (SW) rectangle ($(NW) - (5.0, 0.0)$);
(GammaD) at ($(NW) !.5! (SW) + (7.0, 0)$) $\Gamma_D$;

(loadP) at ($0.5*(NE) + 0.5*(SE)$) ;
[ultra thick] ($(loadP) + (0.0, 2.5)$) – ($(loadP) + (0.0, -2.5)$);

(Gamma) at ($(loadP) - (5, 0)$) $\Gamma$;

[->, ultra thick] (0,0) – (0, 20) node[right]$y$;
[->, ultra thick] (0,0) – (20, 0) node[right]$x$;
(A) at ($(SW) - (0.0, 5.0)$) ;
(B) at ($(SE) - (0.0, 5.0)$) ;
[line style = line width=0.7,
extension start length=-5cm,
extension end length=-5cm](A)(B)$L$;
(C) at ($(NW) - (15.0, 0.0)$) ;
(D) at ($(SW) - (15.0, 0.0)$) ;
[line style = line width=0.7,
extension start length=15cm,
extension end length=15cm](D)(C)$\frac{2L}{5}$;
[line style = line width=0.7,
extension start length=10cm,
extension end length=10cm]($(loadP) + (10, 2.5)$)($(loadP) + (10, -2.5)$)$\frac{4L}{50}$;

Compliance domain. $L=100$

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Uniform mesh

90Non-uniform mesh

Optimized designs for the compliance problem.

Cost function evolution for the compliance problem.

Constraint function evolution for the compliance problem.

Convergence metric evolution for the compliance problem.

We next benchmark our $L^2$ algorithm with the compliant mechanism design problem, cf. Figure <ref> where we minimize the horizontal displacement $u=\mathbf{u}\cdot\mathbf{e}_1$ in the output port $\Gamma_2$, subject to a maximum volume constraint $\hat V = 0.3 |D|$, viz.

\begin{align}
\theta_0 & = \int_{\Gamma_2} u ~da \label{eq:compliant} \,,           \\
\theta_1 & = \int_{D} \hat \nu ~dV - \hat V \,.\label{eq:volume_mech}
\end{align}

Consistent with <cit.>, we introduce Robin boundary conditions into the formulation (<ref>).

\begin{equation}
\begin{aligned}
\mathbf{n} \cdot \mathbb{C}[\nabla \mathbf{u}]  \mathbf{n} & = -k_{in} (\mathbf{u}\cdot \mathbf{n}) + f_x & ~ \text{on} ~\Gamma_1 \,,   \\
\mathbf{n} \cdot \mathbb{C}[\nabla \mathbf{u}]  \mathbf{n} & = -k_{out} (\mathbf{u}\cdot \mathbf{n})      & ~ \text{on} ~\Gamma_2   \,,
\end{aligned}
\end{equation}

where $f_x=10$ and the spring coefficients are $k_{in}=\frac{1}{3}$ and $k_{out}=\frac{0.001}{3}$.
Figures <ref> and <ref> illustrate the design domain and the non-uniform mesh.
The number of elements in the uniform (mechanism_uniform.geo) and non-uniform (mechanism_amr.geo) meshes are 57,600 and 199,404.
The initial design is again $\nu(\mathbf{x}) = 0.1$.
We run all the optimization problems until the number of iterations reaches 200, for which almost all the problems converged to optimized designs.

The optimized designs for the length scale $\kappa=0.8$, summarized and illustrated
in Table <ref> and Figures <ref> - <ref>,
again show that the original $\mathbb{R}^n$ NLP algorithm is mesh dependent.
Notably, the design on the non-uniform mesh with the $\mathbb{R}^n$ algorithm does not even converge in 200 iterations.

[  spring/.style = decorate,decoration=zigzag,amplitude=6pt,segment length=4pt
, scale=0.07, every node/.style=inner sep=0pt]
(NW) at (0,50) ;
[fill=gray!50!white] (0,0) – (100,0) – (100,50) – (0,50) – (0,0);
in 0,4,...,100
[thick] (NW) ++ (,2.0cm) circle (2.0cm);
[thick]  (0,54.5) – (100, 54.5);
[pattern=north west lines]  (0,54.5) rectangle (100, 57);

($(NW) - (0.0, 5.0)$) – ($(NW) - (5.0, 5.0)$);
[spring] ($(NW) - (5.0, 5.0)$) – ($(NW) - (15.0, 5.0)$);
($(NW) - (15.0, 5.0)$) – ($(NW) - (20.0, 5.0)$);
(Gamma1) at ($(NW) - (5, 10)$) $\Gamma_1$;

($(NW) - (20.0, 10.0)$) – ($(NW) - (20.0, 0.0)$);
[very thick] (NW) – ($(NW) + (0.0, -10.0)$);
[pattern=north west lines] ($(NW) - (25.0, 10.0)$) rectangle ($(NW) - (20.0, 0.0)$);

(NE) at (100,50) ;
($(NE) - (0.0, 5.0)$) – ($(NE) + (5.0, -5.0)$);
[spring] ($(NE) + (5.0, -5.0)$) – ($(NE) + (15.0, -5.0)$);
($(NE) + (15.0, -5.0)$) – ($(NE) + (20.0, -5.0)$);
(Gamma2) at ($(NE) + (5, -10)$) $\Gamma_2$;

($(NE) + (20.0, 0.0)$) – ($(NE) + (20.0, -10.0)$);
[very thick] (NE) – ($(NE) + (0.0, -10.0)$);
[pattern=north west lines] ($(NE) + (20.0, 0.0)$) rectangle ($(NE) + (25.0, -10.0)$);

(SW) at (0,0) ;
(SE) at (100,0) ;
(GammaD) at ($(SW) + (5, 2)$) $\Gamma_D$;
[pattern=north west lines] ($(SW) + (0.0, 2.0)$) rectangle ($(SW) - (5.0, 0.0)$);

[line style = line width=0.7,
extension start length=10cm,
extension end length=10cm]($(SW) - (10cm, 0.0)$)($(SW) + (-10cm, 2)$)$2$;

in 0, 2.5, 5.0, 7.5, 10.0
[-Latex[width=1mm]] ($(NW) - (0.0, \y)$) – ($(NW) + (15.0, -\y)$);
at ($(NW) + (20.0, -5.0)$) $f_{x}$;

in 0, 2.5, 5.0, 7.5, 10.0
[color=red, -Latex[width=1mm]] ($(NE) - (0.0, \y)$) – ($(NE) - (15.0, \y)$);
at ($(NE) - (25.0, 5.0)$) $u_{out}$;

(A) at ($(SW) - (0.0, 5.0)$) ;
(B) at ($(SE) - (0.0, 5.0)$) ;
[line style = line width=0.7,
extension start length=-5cm,
extension end length=-5cm](A)(B)120;
(C) at ($(NE) + (50.0, 0.0)$) ;
(D) at ($(SE) + (50.0, 0.0)$) ;
[line style = line width=0.7,
extension start length=-50cm,
extension end length=-50cm](D)(C)60;
(E) at ($(C) + (-10.0, 0.0)$) ;
(F) at ($(E) + (0.0, -10.0)$) ;
[line style = line width=0.7,
extension start length=-40cm,
extension end length=-40cm](F)(E)10;
(G) at ($(NW) + (-30.0, 0.0)$) ;
(H) at ($(G) + (0.0, -10.0)$) ;
[line style = line width=0.7,
extension start length=-30cm,
extension end length=-30cm](G)(H)10;

Design domain and boundary conditions for the compliant mechanism problem.
Domain symmetry is used whereby only the lower half of the structure is analyzed.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Uniform mesh

90Non-uniform mesh

Optimized designs for the compliant mechanism problem.

Cost function evolution for the compliant mechanism problem.

Constraint function evolution for the compliant mechanism problem.

Convergence metric evolution for the compliant mechanism problem.

Our next example is the stress-constrained problem where the goal is to minimize the volume of an L-bracket subject to a maximum
pointwise constraint in the Von Mises stress field $\sigma_{VM} \leq \sigma_y$, cf. Figure <ref>.
We follow the formulation in <cit.>, where the stress constraint is imposed via a penalty method, i.e.

our unconstrained problem uses the cost function

\begin{align}
\theta_0 & = \int_{D} \hat \nu ~da+ \gamma  \norm{\sigma_{VM} - \sigma_y}_+  \,,
\end{align}

where the penalty parameter is $\gamma=10$,

\begin{align}
\norm{\sigma_{VM} - \sigma_y}_+ = \int_{D} R_p\left(\frac{\sigma_{VM}}{\sigma_y}\right)dV \,,
\end{align}

\begin{align}
R_p(x) = (1 + (x)^p)^{\frac{1}{p}} - 1\,,
\label{eq:penalty_func}
\end{align}

$\sigma_y=1.5$ and $p=8$, for the first 300 iterations and $p=20$ for the remaining 100.
The relaxed stress formulation <cit.> uses

\begin{align}
\sigma_{\eta} = \eta_c \mathbb{C}[\nabla \mathbf{u}] \,,
\end{align}

to calculate the Von Mises stress $\sigma_{VM} = \sqrt{\frac{2}{3} \sigma_{\eta}^{'}:\sigma_{\eta}^{'}}$, where $\sigma_{\eta}^{'} = \sigma_{\eta} - \frac{2}{3} \text{tr}(\sigma_{\eta}) \textbf{I}$ and

\begin{align}
\eta_c (\hat\nu) = \hat\nu^{0.5}\,,
\end{align}

The filter parameter is $\kappa=1.2$ and the initial design is $\nu(\mathbf{x}) = 0.5$.
To obtain a better design we extend the domain by adding a region $\Omega_i$ of finite elements near the reentrant corner, cf. Figure <ref>.
This region however, is excluded in the design by enforcing the constraint $\int_{\Omega_i} \hat{\nu}~dV \leq 0$.
This added region lessens the boundary effect of the filter operation in the reentrant corner region which otherwise adversely affects our results <cit.>.
The non-uniform mesh (file lbracket_amr.geo) contains 53,122 elements, whereas the uniform mesh (lbracket_uniform.geo) contains 139,264 elements.

The optimized designs in Table <ref> again illustrate the mesh dependence of the original $\mathbb{R}^n$ NLP algorithm.
This time however, the difference in the cost function values between the $\mathbb{R}^n$ and $L^2$ designs is barely noticeable, cf. the log plot in Figure <ref>.
Due to the high oscillation in the convergence metric, we plot them separately in Figures <ref> - <ref>.

[every node/.style=inner sep=0pt]
[pattern=north west lines]  (-0.05,) rectangle (1.05, 1.05);
[very thick] (-0.05,) – (1.05,);
[fill=gray!50!white] (0,0) – (,0) – (,) – (,) – (,) – (0,) – (0,0);
(GammaD) at (0.5, 0.95) $\Gamma_D$;

(A) at (0, -0.5) ;
(B) at (, -0.5) ;
[line style = line width=0.7,
extension start length=-1cm,
extension end length=-1cm](A)(B)$100$;
(C) at (-0.5, 0) ;
(D) at (-0.5, ) ;
[line style = line width=0.7,
extension start length=1cm,
extension end length=1cm](C)(D)$100$;
(F) at (0, 1.1) ;
(E) at (, 1.1) ;
[line style = line width=0.7,
extension start length=1cm,
extension end length=1cm](F)(E)$40$;
(G) at (1.05, 0) ;
(H) at (1.05, ) ;
[line style = line width=0.7,
extension start length=-1cm,
extension end length=-1cm](G)(H)$40$;
(Gamma) at (0.97, 0.95) $\Gamma_N$;
(L1) at (, ) ;
(L1p) at (, 1.2) ;
(L2) at (0.97, ) ;
(L2p) at (0.97, 1.2) ;
(L3) at (0.94, ) ;
(L3p) at (0.94, 1.2) ;
[thick,->] (L1p) – (L1);
[thick,->] (L2p) – (L2);
[thick,->] (L3p) – (L3);
(L4p) at (0.94, 1.4) ;
(L5p) at (, 1.4) ;
[line style = line width=0.7,
extension start length=1.6cm,
extension end length=1.6cm](L4p)(L5p)$5$;

Intended design domain for the stress-constrained problem.

Cost function evolution for the stress-constrained problem.

Convergence metric evolution for the stress-constrained problem with uniform mesh in $L^2$.

Convergence metric evolution for the stress-constrained problem with uniform mesh in $\mathbb{R}^n$.

Convergence metric evolution for the stress-constrained problem with non-uniform mesh in $\mathbb{R}^n$.

Convergence metric evolution for the stress-constrained problem with non-uniform mesh in $L^2$.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Uniform mesh

90Non-uniform mesh

Optimized designs for the stress-constrained problem.

§.§ Uniform refinement

We return to the compliance problem for our next example, but in three dimensions, cf. Figure <ref>.
The volume constraint is $\hat{V}= 0.15 |D|$,
the surface $\Gamma_N$ is subject to a traction $\mathbf{t} = -1.0 \mathbf{e}_2$, the length scale parameter $\kappa=5 \times 10^{-5}$ and the initial design is $\nu(\mathbf{x})=0.15$.
To prevent the iterative solver from diverging, we set $\epsilon_{\nu}$ to $10^{-4}$.

Design domain for the three dimensional compliance problem.

Instead of comparing two different meshes, we compare the influence of the mesh refinement
starting from the same uniform mesh with 20x10x10 hexahedral elements.
As such, we perform an optimization with the initial mesh uniformly refined once, twice,
three and four times, for a total of 16,000, 128,000, 1,024,000 and 8,192,000 elements.
The construction and refinement of the mesh is done by the utility meshing functions in Firedrake within the code shared along with this paper.
Each hexahedral element is split in eight at each refinement level.
We plot the evolution of the cost function for 1000 iterations corresponding to
the $L^2$ and $\mathbb{R}^n$ NLP algorithms for each refinement level in Figure <ref>.
The behavior of the $\mathbb{R}^n$ NLP algorithm clearly depends on the refinement level, whereas the $L^2$ does not.
Notably, after 1000 iterations, the optimized design with four levels of refinement take compliances values of 10.06 and 2.91 for the $\mathbb{R}^n$ and $L^2$ NLP algorithms.

Cost function evolution of the compliance problem in three dimensions with 1, 2, 3 and 4 uniform refinements.
Results plotted in log-log scale to highlight the differences.

We compare the four level of refinement designs obtained with the $\mathbb{R}^n$ and the $L^2$ NLP algorithms in Figure <ref>.
Upon inspection of the designs, we noticed that the $\mathbb{R}^n$ algorithm fails to reach the lower bound of the volume fraction (0);
it never dips below $\nu=10^{-2}$ for the void phase.
Due to the volume constraint and the “heavier” void phase, the $\mathbb{R}^n$ optimizer cannot add more
mass to the structure and therefore, the compliance is higher.
We conjecture that the reason behind this is due to the difference between the derivative and the gradient.
The $\mathbb{R}^n$ algorithm uses $D\theta$ whereas the $L^2$ uses $\nabla\theta = \mathbf{M}^{-1}D\theta$.
Thus, although $D\theta$ and $\nabla \theta$ are parallel since $\mathbf{M} = |\Omega_e| \mathbf{I}$, the sensitivity of the cost and constraint function of small elements is less influential
in the $\mathbb{R}^n$ vs $L^2$ algorithm.
We further conjecture that using an interpolation scheme with nonzero derivative values for $\nu=0$,
such as RAMP , could alleviate this issue.

Optimized designs for the compliance problem in three dimensions with four levels of refinement, thresholded with volume fraction $\hat\nu$ greater than 0.5

§.§ Adaptive mesh refinement

Our $L^2$ GCMMA implementation is specially suitable for use with AMR strategies during the optimization <cit.>.
In our next two examples, we show the algorithm utility when applying AMR during the optimization and compare it with the original $\mathbb{R}^n$ implementation.
There are three questions to address when applying AMR in a topology optimization problem:
where in the domain we apply AMR, what kind of AMR (coarsening, refinement or both)
and when to apply the AMR during the optimization.
We use an element-based error quantity to determine the regions subject to AMR.
The details are explained in <ref>.
We use two AMR strategies: only refinement or only coarsening.
We apply the AMR after a pre-determined number of iterations
during the optimization[We are not concerned with applying an optimal strategy for AMR during the optimization.
See <cit.> for a more sophisticated scheme.].
Table <ref> summarizes four AMR schemes wherein
it is seen that refinement occurs either during early iterations as the
design evolves or later once the design is well defined.
Coarsening only occurs late in the design evolution as well.
All problem/algorithm combinations are run for 600 iterations.

StrategyAMR type    Coarsening     Refinement
100, 150

10, 80

150, 200

100, 150

Four AMR strategies. Iteration number after which AMR is applied.

Our first AMR example is a repeat of our cantilever example,
but applying coarsening and refinement to the Table <ref> strategies.
The optimized designs, shown in Tables <ref>
and <ref> for refinement and coarsening respectively,
are different when using the $\mathbb{R}^n$ algorithm.
The cost function evolution in Figures <ref>
and <ref> further reflect their difference.
The meshes for the optimized designs illustrated in the Appendix, cf. Figures <ref>,
<ref>, highlight the mesh independence of the
the $L^2$ versus $\mathbb{R}^n$ algorithm.

The cost function evolution for coarsening, cf. Figure <ref>,
is not affected by the different refinement strategies even with the $\mathbb{R}^n$ optimization.
Most likely, this is because both $A$ and $B$ strategies are applied at similar iteration numbers,
after which the designs have almost converged.
It is important to highlight that the algorithm in $L^2$ still performs better.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Optimized designs for the compliance problem with AMR refinement only.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Optimized designs for the compliance problem with AMR coarsening only.

Cost function evolution for the compliance problem with AMR coarsening.

Cost function evolution for the compliance problem with AMR refinement.

Next we solve a coupled thermal flow problem similar than in <cit.>
to demonstrate the algorithm's application to more complex physical phenomena.
The domain in Figure <ref> is the cross section of a heat exchanger, where the design variable $\nu$ represents the volume fraction of a heat generating solid material that needs to be distributed to control the temperature and the fluid flow.
The goal of the optimization is to maximize the heat generated in the solid while keeping the pressure drop in the flow from the inlet $\Gamma_1$ to the outlet $\Gamma_2$ lower than a fixed value $P_{\text{drop}}$.
As such, the optimization problem reads

[ every node/.style=inner sep=0pt, spring/.style = decorate,decoration=zigzag,amplitude=6pt,segment length=4pt, scale=5.0 ]
(A) at (-0.1, 0.0) ;
(B) at (-0.1, -0.1) ;
(C) at (0, -0.1) ;
(D) at (0.0, -0.5) ;
(E) at (1.0, -0.5) ;
(F) at (1.0, -0.1) ;
(G) at (1.1, -0.1) ;
(H) at (1.1, 0.0) ;
(J) at (0.0, 0.0) ;
($(A)$) – ($(B)$) – ($(C)$) – ($(D)$) – ($(E)$) – ($(F)$) – ($(G)$) – ($(H)$) – ($(A)$);
[line style = line width=0.7,
extension start length=-0.1cm,
extension end length=-0.1cm]($(A) + (-0.1, 0)$)($(B) + (-0.1, 0.0)$)$\frac{L}{10}$;
[line style = line width=0.7,
extension start length=0.1cm,
extension end length=0.1cm]($(D) + (-0.2, 0)$)($(B) + (-0.1, 0.0)$)$\frac{2L}{5}$;
[line style = line width=0.7,
extension start length=0.1cm,
extension end length=0.1cm]($(E) + (0.0, -0.1)$)($(D) + (0.0, -0.1)$)$L$;
[line style = line width=0.7,
extension start length=-0.1cm,
extension end length=-0.1cm]($(A) + (0.0, -0.2)$)($(J) + (0.0, -0.2)$)$\frac{L}{10}$;
(GammaD) at ($(A) !.5! (B) + (0.1, 0)$) $\Gamma_1$;
(GammaF) at ($(G) !.5! (H) + (0.1, 0)$) $\Gamma_2$;
(Domain) at (0.5, -0.1) $D$;
[ultra thick, draw=black, fill=black, opacity=0.4] (0.5, -0.3) ellipse (0.2cm and 0.1cm);
(Domain) at (0.5, -0.31) $\nu=1$;

Thermal flow domain. $L$=1.

\begin{align}
\underset{\nu \in V}{\text{max}}~ J(\nu)          & = \int_{D} \nu B (1 - T)~dV \,, \\
\text{s.t.} ~\left(\ww, p, T \right) \in \VV{} \times Q \times W ~\text{satisfy} ~      \\
F(\nu, \ww, p; \vv, q)  & = 0 \label{eq:flow_eq}\,, \\
a_T(\ww; T, m) + c_T(\ww; T, m)                                & = 0     \label{eq:thermal_eq}                                                    \,,     \\
\text{for all} ~\left(\vv, q, m \right) \in \VV{0} \times Q  \times W_0                \\
G(\nu)                                             & = \int_{\Gamma_{1}} p ~dA - \int_{\Gamma_{2}} p ~dA \leq P_{\text{drop}}      \,,
\end{align}

where $(\mathbf{w}, p, T)$ are the flow velocity, pressure and temperature and $(\mathbf{v}, q, m)$ are their admissible counterparts.
Equation (<ref>) is the weak form of the Navier-Stokes equation wherein

\begin{equation}
\begin{aligned}
F(\nu, \ww, p; \vv, q) & = \int_{D}\left( \left(\ww \cdot {\nabla} \ww \right) \cdot \vv + \frac{1}{Re} {\nabla} \ww	: {\nabla} \vv
+ \frac{1}{Da} r(\nu) \ww \cdot \vv \right) ~dV\\
&+ \int_D \left( p \nabla \cdot \vv + q \nabla \cdot \ww \right)~dV                                                    \,,      \\
\end{aligned}
\label{eq:stokes_expanded}
\end{equation}

with Reynolds number $Re=1.0$ and Darcy number $Da=10^{-6}$.
The RAMP function <cit.> with $q_{\text{RAMP}} = 20.0$, i.e.

\begin{align}
r(\nu) = \frac{\nu}{1+q_{\text{RAMP}}(1-\nu)}
\end{align}

is used to obtain discrete 0-1 designs.
At the inlet $\Gamma_1$, the Dirichlet condition $\ww = \ww_1$
is a horizontal parabolic profile with a maximum non-dimensional velocity $W_{\text{max}}=1$.
At the outlet $\Gamma_2$, a traction-free condition is imposed.
The remaining boundary $\Gamma \setminus(\Gamma_1 \cup \Gamma_2)$ has no-slip condition $\ww = 0$.

Equation (<ref>) is the weak form of the advection-diffusion equation wherein

\begin{equation}
\begin{aligned}
a_T(\nu; T, m) & = \int_{D} m \ww \cdot {\nabla} T + \frac{1}{Pe} {\nabla} T	: {\nabla} m  - \nu B (1 - T) m ~dV \,, \\
\label{eq:temp_gls}
\end{aligned}
\end{equation}

We assume the same thermal conductivities in the solid and and the fluid,
a Peclet number $Pe=10^4$.
The heat source in the solid $\nu B (1 - T)$ is proportional to the difference between a reference temperature 1 and the local temperature $T$, and the non-dimensional heat generation coefficient $B=0.01$ (more details on the heat source are in <cit.>).
The Galerkin Least Squares (GLS) stabilization term

\begin{equation}
\begin{aligned}
c_T(\ww; T, m)                & = \int_D \tau_{GLS}\mathcal{L}_T(T) \cdot \mathcal{L}_T(m) ~dV
\end{aligned}
\end{equation}

stabilizes the otherwise highly oscillatory boundary layers due to the fluid velocity field,

\begin{equation}
\begin{aligned}
\tau_{GLS} = \beta_{GLS} \left( \frac{4 \ww \cdot \ww}{h^2} + \left(9 \frac{4 }{h^2Pe} \right)^2 \right)^{-0.5} \,,
\end{aligned}
\end{equation}

$h$ is the element cell size and

\begin{align}
\mathcal{L}_T(T) = \ww \cdot {\nabla} T + \frac{1}{Pe} \nabla T  - \nu B (1 - T)\,.
\end{align}

is the residual.
We use $\beta_{GLS} = 0.9$ in our examples.
We apply a constant temperature $T=0$ on $\Gamma_1$ and adiabatic boundary conditions over all surfaces with the exception of $\Gamma_1$ and $\Gamma_2$.
The finite element discretization of $(\mathbf{w}, p, T)$ uses linear Lagrange elements.
The cost function $J(\nu)$ aims to maximize the heat generation in the design domain
and the constraint $G(\nu) < P_{\text{drop}}=70.0$ limits the pressure drop in the system
and serves to regularize the problem as it imposes an upper bound on the fluid-solid
interface where $\nabla \nu \neq 0$.

Optimized designs are shown in Tables <ref> and <ref> for refinement and coarsening strategies.
Again, the designs obtained using the $\mathbb{R}^n$ algorithm differ.
The cost function evolution in Figures <ref> and <ref> reflect this dependency as well.
The meshes for the optimized designs illustrated in the Appendix, cf. Figures <ref> and <ref> highlight the mesh independence of the $L^2$ algorithm.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Optimized designs for the thermal flow problem with AMR refinement only.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Optimized designs for the thermal flow problem with AMR coarsening only.

Cost function evolution for the thermal flow problem with AMR coarsening.

Cost function evolution for the thermal flow problem with AMR refinement.

§ CONCLUSION

In this work, we presented the necessary mathematical concepts to understand
the relationship between the domain discretization and the NLP algorithm and applied them to the GCMMA algorithm.
Our $L^2$ GCMMA implementation is benchmarked with several problems in
topology optimization and is able to obtain mesh independent designs starting
from the same initial designs, while the original $\mathbb{R}^n$ NLP algorithm is not.
We first showed how the new algorithm solves ill-conditioned optimization problems where the $\mathbb{R}^n$ algorithm fails.
Then we showed its efficiency when solving large scale problems over uniform meshes.
Lastly, we illustrated the effectiveness of the $L^2$ algorithm when applying AMR during the optimization for two problems with different physics.
For future work, the algorithm can be extended to handle design fields in
other common spaces in topology optimization such as $H^1$, i.e. for nodal design variables or for B-splines.
Lastly, a rigorous mathematical proof is necessary to ensure the NLP algorithm is
mathematical sound for all corner cases.

§ REPLICATION OF RESULTS

The scripts used in this article are archived in <cit.> and require the pyMMAopt library <cit.>

§ ACKNOWLEDGEMENTS

This work was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract DE-AC52-07NA27344. The author thanks the Livermore Graduate Scholar Program for its support.
On behalf of all authors, the corresponding author states that there is no conflict of interest. LLNL-JRNL-820905.

§ APPENDIX

§.§ GCMMA

We assume existence and uniqueness of the MMA subproblem (<ref>) and solve it using a primal-dual interior-point NLP algorithm as in <cit.>, where a sequence of relaxed KKT conditions are solved using the Newton-Raphson method.
Omitting the iteration indices $(k, j)$ for clarity we solve the MMA subproblem (<ref>) which is restated here as

\begin{equation}
\begin{aligned}
& \underset{
\subalign{
\nu            & \in V                                                                                \\
\boldsymbol{y} & \in \mathbb{R}^m                                                                                \\
z              & \in \mathbb{R}                                                                                \\
\boldsymbol{s} & \in \mathbb{R}^m
&                  & \int_{D} \left(\frac{p_{o}}{U - \nu}+\frac{q_{o}}{\nu - L}\right) dV +a_0z +  \sum\limits_{i=1}^m \left( c_i y_i + \frac{1}{2}d_i y_i^2 \right) + r_0 \,,                                                \\
& \text{s.t.}
&                  & \int_{D} \left(\frac{p_{i}}{U - \nu}+\frac{q_{i}}{\nu - L}\right) dV - a_i z - y_i + r_i + s_i = 0, \; i = 1, \ldots, m \,,                                                                              \\
&                  &                                                                                & \alpha \leq \nu \leq \beta \; \text{a.e.}\,, \\
&                  &                                                                                & \boldsymbol{y} \geq 0 \,,                    \\
&                  &                                                                                & z\geq 0 \,,                                  \\
&                  &                                                                                & \boldsymbol{s}\geq 0 \,,
\end{aligned}
\label{eq:subprob_cost}
\end{equation}

where we have introduced the slack variable $\boldsymbol{s} \in \mathbb{R}^m$ to transform the  $\theta_i(\nu) \leq 0$ inequality constraints into $m$ equality constraints.

To solve the above, we introduce log-barrier terms for the inequality constraints.
The Lagrangian for the resulting penalized problem is

\begin{equation}
\begin{aligned}
\Lagr (\nu,\mathbf{y},z,\lambda, \mathbf{s}) = &
\int_{D} \left(\frac{p_{o}}{U - \nu}+\frac{q_{o}}{\nu - L}\right) dV +a_0z +  \sum\limits_{i=1}^m \left( c_i y_i +  \frac{1}{2}d_i y_i^2 \right) + r_0                                                                                \\
& +\sum\limits_{i=1}^m \lambda_i \left( \int_{D} \left(\frac{p_{i}}{U - \nu}+\frac{q_{i}}{\nu - L}\right) dV - a_i z - y_i + r_i + s_i \right)                                                   \\
& - \epsilon \int_{D} \text{ln}( \nu - \alpha) ~dV - \epsilon \int_{D}  \text{ln} (\beta - \nu) dV - \epsilon \sumonm \text{ln} y_i - \epsilon \text{ln}(z) - \epsilon \sumonm \text{ln} s_i \,,
\end{aligned}
\label{eq:lagrasub}
\end{equation}

where $\lambda_i;~i=1\dots m$ are the Lagrange multipliers for the inequality constraints.

Stationarity of the Lagrangian is then solved sequentially for different values of the barrier parameter $\epsilon$, starting at $\epsilon=1$ and decreasing by a factor of 0.1 until $\epsilon \leq 10^{-5}$.
The MMA subproblem iteration loop terminates when the norm of the derivative of $\Lagr$ is less than $0.9\epsilon$.
Note that the bound constraints on the volume fraction field $\nu\in V$ are enforced via integrals instead of via summations as in the original NLP algorithm.

To keep the notation clean, we group the variables as $\chi=\left( \nu,\mathbf{y},z,\boldsymbol\lambda, \mathbf{s}\right)$ and express the function $\psi_j (\nu): V \rightarrow V,~ j = 0, \ldots, m $ as

\begin{equation}
\psi_j (\nu) = \frac{p_{j}}{U - \nu} + \frac{q_{j}}{\nu - L},~ j = 0, \ldots, m  \,,
\label{eq:psioperator}
\end{equation}

In this way, the stationary conditions on the Lagrangian become

\begin{equation}
\begin{aligned}
D_{\nu}\Lagr (\chi)[ \delta \nu]             & = \int_D \left( D\psi_0(\nu)[\delta \nu] + \sumonm \lambda_i D\psi_i(\nu)[\delta \nu] - \frac{\epsilon \delta \nu}{\nu - \alpha} + \frac{\epsilon \delta \nu}{\beta - \nu} \right) ~ dV & = 0 & ~~\forall \delta \nu \in V \,,            \\
D_{y_i}\Lagr (\chi)[ \delta y_i]             & = \left( c_i + d_i y_i - \lambda_i - \frac{\epsilon}{y_i} \right) \delta y_i                                                                                & = 0 & ~~\forall \delta y_i  ~i=1...m \,,        \\
D_z\Lagr (\chi)[ \delta z]                   & = (a_0 - \sumonm \lambda_i a_i - \frac{\epsilon}{z}) \delta z                                                                                & = 0 & ~~\forall \delta z \,,                    \\
D_{\lambda_i}\Lagr (\chi)[ \delta \lambda_i] & = \left(\int_D \psi_i(\nu) ~dV - a_i z - y_i +r_i + s_i \right) \delta \lambda_i                                                                                & = 0 & ~~ \forall \delta \lambda_i ~i=1...m  \,, \\
D_{s_i}\Lagr (\chi)[ \delta s_i]             & = (\lambda_i - \frac{\epsilon}{s_i})\delta s_i                                                                                & = 0 & ~~\forall \delta s_i ~i=1...m  \,,
\end{aligned}\label{eq:lagra1}
\end{equation}

Equation (<ref>).b is fulfilled for any $\delta y_i$, therefore the term in the parenthesis equals zero and similarly for Equations (<ref>).c-(<ref>).e.
From the last Equation (<ref>).e, we obtain $\lambda_i s_i = \epsilon\;\text{for}~i=1...m$.

We next introduce the variables

\begin{equation}
\begin{aligned}
\varepsilon & = \frac{\epsilon}{\nu - \alpha}  \,,    \\
\eta        & = \frac{\epsilon}{\beta - \nu}  \,,     \\
\mu_i       & = \frac{\epsilon}{y_i}\,, ~i=1...m  \,, \\
\zeta       & = \frac{\epsilon}{z} \,,
\end{aligned}
\end{equation}

where $\varepsilon, \eta \in V$ and $\mu_i, \zeta \in \mathbb{R}$ to generate the following system of nonlinear equations

\begin{equation}
\begin{aligned}
D_{\nu}\Lagr (\chi)[\delta \nu] & = \int_D \left( D\psi_0(\nu)[\delta \nu] + \sumonm \lambda_i D\psi_i(\nu)[\delta \nu] - \varepsilon\delta \nu + \eta\delta \nu \right) ~ dV & = 0 & ~~\forall \delta \nu \in V \,, \\
D_{y_i}\Lagr (\chi)             & =  c_i + d_i y_i - \lambda_i - \mu_i                                                                                & = 0 & ~i=1...m \,,                   \\
D_z\Lagr (\chi)                 & =  a_0 - \sumonm \lambda_i a_i  - \zeta                                                                                & = 0 & \,,                            \\
D_{\lambda_i}\Lagr (\chi)       & = \int_D \psi_i(\nu) ~dV - a_i z - y_i +r_i + s_i                                                                                & = 0 & ~i=1...m  \,,                  \\
\lambda_i s_i                   & =\epsilon                                                                                &     & ~ ~i=1...m \,,                 \\
(\nu - \alpha) \varepsilon      & =\epsilon   \,,                                                                                \\
(\beta - \nu) \eta              & =\epsilon  \,,                                                                                \\
\mu_i y_i                       & =\epsilon                                                                                &     & ~ ~i=1...m \,,                 \\
\zeta z                         & =\epsilon  \,.
\end{aligned}
\label{eq:kkt_subproblem}
\end{equation}

We express the above in a more compact form $F(\upsilon)=0$, where $\upsilon=\left( \nu,\mathbf{y},z,\boldsymbol\lambda, \mathbf{s}, \varepsilon, \eta, \boldsymbol\mu, \zeta\right) $ and
$F:V \times \mathbb{R}^m \times \mathbb{R}\times \mathbb{R}^m\times \mathbb{R}^m \times V \times V \times \mathbb{R}^m \times \mathbb{R}
\rightarrow V^* \times \mathbb{R}^m \times \mathbb{R}\times \mathbb{R}^m\times\mathbb{R}^m \times V \times V \times \mathbb{R}^m \times \mathbb{R}  $ is defined as

\begin{equation}
F(\upsilon) =
\left[	\begin{array}{l}
\int_D \left( D\psi_0(\nu)[\delta \nu] +\textstyle \sumonm \lambda_i D\psi_i(\nu)[\delta \nu] - \varepsilon \delta \nu + \eta \delta \nu\right) ~dV \\
c_i + d_i y_i - \lambda_i - \mu_i                                                                                \\
a_0 - \sumonm\lambda_i a_i   - \zeta                                                                                \\
\int_D \psi_i(\nu) ~dV - a_i z - y_i +r_i + s_i ~                                                                                \\
\lambda_i s_i -\epsilon                                                                                \\
\int_D ((\nu - \alpha) \varepsilon - \epsilon) \delta \varepsilon ~ dV                                                                              \\
\int_D ((\beta - \nu) \eta - \epsilon ) \delta \eta ~ dV                                                                                \\
\mu_i y_i -\epsilon                                                                                \\
\zeta z -\epsilon
\end{array} \right]
\begin{bmatrix}
\delta_{\nu}         \\
\delta_{y_i}         \\
\delta_{z}           \\
\delta_{\lambda_i}   \\
\delta_{s_i}         \\
\delta_{\varepsilon} \\
\delta_{\eta}        \\
\delta_{\mu_i}       \\
\delta_{\zeta}
\end{bmatrix}
\label{eq:newton_system}
\end{equation}

for all $\delta \nu \in V$, $\delta \varepsilon \in V$ and $\delta \eta \in V$.
As seen above, for discretization purposes, we enforce equations (<ref>).f and (<ref>).g weakly.

We solve $F(\upsilon)=0$ (for all $\delta \upsilon$) via Newton-Raphson by linearizing around the iterate $\upsilon^{(l)}$ and requiring the update $\Delta \upsilon^{(l)}$ to satisfy

\begin{equation}
D F(\upsilon^{(l)}) \Delta \upsilon^{(l)} = - F(\upsilon^{(l)})  \,,
\end{equation}

which we expand into

\begin{equation}
\left[	\begin{array}{l}
\int_D \left( D^2\psi_0(\nu)[\delta \nu,\Delta \nu] + \sumonm \lambda_i D^2\psi_i(\nu)[\delta \nu, \Delta \nu] +\sumonm D\psi_i(\nu)[\delta \nu] \Delta \lambda_i - \delta \nu \Delta \varepsilon + \delta \nu \Delta \eta \right) ~dV \\
d_i \Delta y_i - \Delta \lambda_i - \Delta \mu_i                                                                                \\
- \sumonm a_i \Delta \lambda_i - \Delta \zeta                                                                                \\
\int_D D\psi_i(\nu)[\Delta \nu] ~dV - a_i \Delta z - \Delta y_i + \Delta s_i                                                                                \\
\Delta \lambda_i s_i + \lambda_i \Delta s_i                                                                                \\
\int_D \left( \Delta \nu \varepsilon  + (\nu - \alpha ) \Delta \varepsilon \right) \delta \varepsilon ~dV                                                                                \\
\int_D \left(- \Delta \nu  \eta  + (\beta - \nu) \Delta \eta \right) \delta \eta ~dV                                                                                \\
\Delta \mu_i y_i + \mu_i \Delta y_i                                                                                \\
\Delta \zeta z  +    \zeta \Delta z
\end{array} \right]
\begin{bmatrix}
-\delta_{\nu}         \\
-\delta_{y_i}         \\
-\delta_{z}           \\
-\delta_{\lambda_i}   \\
-\delta_{s_i}         \\
-\delta_{\varepsilon} \\
-\delta_{\eta}        \\
-\delta_{\mu_i}       \\
\end{bmatrix}
\label{eq:newtonlinearsystem}
\end{equation}

for all $\delta \nu, \delta \varepsilon$ and $\delta \eta$ in $V$.

We proceed to discretize the fields $\nu, \varepsilon$ and $\eta$ in $V$ to be piecewise uniform over the finite elements as in Equation (<ref>), i.e. using the piecewise uniform basis functions $\mathscr{P} = \left\{\phi_1, ..., \phi_n \right\}$.
So that, e.g.

\begin{equation}
\begin{aligned}
\Delta \nu(\mathbf{x}) = \boldsymbol{\Delta \mathbf{\nu}}^T \boldsymbol{\phi} (\mathbf{x})\,, \\
\delta \nu (\mathbf{x})= \boldsymbol{\delta \mathbf{\nu}}^T \boldsymbol{\phi}(\mathbf{x})\,.
\end{aligned}
\end{equation}

Substituting these expressions in Equations (<ref>) and (<ref>) and using the arbitrariness of $\delta \nu, \delta \varepsilon$ and $\delta \eta$ yields the discretized residuals

\begin{equation}
\begin{aligned}
\boldsymbol{\delta_{\nu}}         & = \int_D \left( D\psi_0(\nu) +\textstyle \sumonm \lambda_i D\psi_i(\nu) - \varepsilon + \eta\right) \boldsymbol{\phi} ~dV   \,,                     \\
\boldsymbol{\delta_{\varepsilon}} & =     (\boldsymbol{\nu} - \boldsymbol{\alpha}) \circ \boldsymbol{\varepsilon}  - \epsilon \boldsymbol{1} \,, \label{eq:epsi_hadamard} \footnotemark \\
\boldsymbol{\delta_{\eta}}        & =     (\boldsymbol{\beta} - \boldsymbol{\nu} ) \circ \boldsymbol{\eta}  - \epsilon\boldsymbol{1}  \,,
%    \delta_{\boldsymbol\mu} &=     \left( \mu_1 y_1 -\varepsilon, \mu_2 y_2 -\varepsilon,..., \mu_m y_m -\varepsilon \right)^T  & ~i=1...m \,,\\
\end{aligned}
\end{equation}

Note that for the choice of $\boldsymbol\phi$ being piecewise uniform over the elements $\int_D \left( \varepsilon (\nu - \alpha) - \epsilon \right) \boldsymbol\phi ~dV = \text{diag}(|\Omega_e|^1, |\Omega_e|^2, ...,) \left( \boldsymbol{\varepsilon} \circ (\boldsymbol{\nu} - \boldsymbol{\alpha}) - \epsilon \boldsymbol{1} \right)$. Equation (<ref>).b follows since $\text{diag}(|\Omega_e|^1, |\Omega_e|^2, ...,)$ is invertible.
where $\boldsymbol{1}=(1,1,...,1)$ and the operator $\circ$ denotes the Hadamard product, i.e. the component wise multiplication between two vectors.
We also define

\begin{equation}
\begin{aligned}
\boldsymbol{\delta_y}         & =  \boldsymbol{c} + \boldsymbol{d} \circ \boldsymbol{y} - \boldsymbol{\lambda} - \boldsymbol{\mu} \,,                                            \\
\delta_{z}                    & = a_0 - \boldsymbol{\lambda^T a} - \zeta \,,                                                                                \\
\boldsymbol{\delta_{\lambda}} & = \left[\int_D \psi_1(\nu) ~dV, ... ,\int_D \psi_m(\nu) ~dV \right]^T - \boldsymbol{a} z - \boldsymbol{y} + \boldsymbol{r} + \boldsymbol{s}  \,, \\
\boldsymbol{\delta_{s}}       & = \boldsymbol{\lambda} \circ \boldsymbol{s} - \epsilon \boldsymbol{1} \,,                                                                        \\
\boldsymbol{\delta_{\mu}}     & = \boldsymbol{\mu} \circ \boldsymbol{y} - \epsilon \boldsymbol{1} \,,                                                                            \\
\delta_{\zeta}                & = \zeta z - \epsilon \,.
\end{aligned}
\end{equation}

In this way, the discretized update equations read

\begin{equation}
\begin{aligned}
\boldsymbol\Psi \boldsymbol\Delta \boldsymbol{\nu}
+ \ggb \boldsymbol{\Delta \lambda}
- \mathbf{M} ( \boldsymbol{\Delta \varepsilon} - \boldsymbol{\Delta \eta})                                                                          & =- \boldsymbol{\delta_{\nu}}  \,,         \\
\langle \boldsymbol{d} \rangle \boldsymbol{\Delta y} - \boldsymbol{\Delta \lambda} - \boldsymbol{\Delta \mu}                                        & =- \boldsymbol{\delta_{y}}\,,             \\
- \mathbf{a}^T \boldsymbol{\Delta \lambda} - \Delta \zeta                                                                                & =- \delta_z  \,,                          \\
\ggb^T \boldsymbol{\Delta \nu} - \boldsymbol{\Delta y} - \boldsymbol{a} \Delta z + \boldsymbol{\Delta s}                                            & =- \boldsymbol{\delta_{\lambda}} \,,      \\
\langle \boldsymbol s \rangle \Delta \boldsymbol\lambda  + \langle \boldsymbol\lambda \rangle \Delta \boldsymbol s                                  & =- \boldsymbol\delta_{\boldsymbol{s}} \,. \\
\langle \boldsymbol\varepsilon \rangle \boldsymbol{\Delta \nu} + \langle \boldsymbol\nu - \boldsymbol\alpha \rangle \boldsymbol{\Delta \varepsilon} & =- \boldsymbol{\delta_{\varepsilon}} \,,  \\
-\langle \boldsymbol\eta \rangle \boldsymbol{\Delta \nu} + \langle \boldsymbol\beta - \boldsymbol\nu \rangle \boldsymbol{\Delta \eta}               & =- \boldsymbol{\delta_{\eta}} \,,         \\
\langle \boldsymbol\mu \rangle \Delta \boldsymbol y + \langle \boldsymbol y \rangle \Delta \boldsymbol\mu                                           & =-\boldsymbol\delta_{\boldsymbol\mu} \,,  \\
\zeta \Delta z + z\Delta \zeta                                                                                & =- \delta_{\zeta} \,,
\end{aligned}
\label{eq:discretized_residual}
\end{equation}

where the operator $\langle \cdot \rangle$ is a diagonal matrix, e.g. $\langle \boldsymbol \nu \rangle = \text{diag}(\nu_1, \nu_2, ..., \nu_n)$ and
$\boldsymbol \Psi$ and $\mathbf{G}$ are the block matrices

\begin{equation}
\begin{aligned}
\boldsymbol\Psi & = \int_D \left( D^2\psi_0(\nu) + \sumonm \lambda_i D^2\psi_i(\nu) \right) \bsi \bsi^T ~dV \,,                                       \\
\mathbf{G}      & = \left[\int_D D\psi_1(\nu) \bsi ~dV, \int_D D\psi_2(\nu) \bsi ~dV, ...~, \int_D D\psi_m(\nu) \bsi ~dV \right]\label{eq:matrixG} \,
\end{aligned}
\end{equation}

where the integrals are calculated using one quadrature point per element.
Each component of the matrix $\mathbf{G}$ is a column vector of dimensions $n \times 1$ due to the basis functions $\boldsymbol \phi$.
Matrix $\boldsymbol \Psi$ has dimensions $n \times n$ and $\mathbf{G}$, $n \times m$.
Ultimately, we are left with the update equation

\begin{equation}
\begin{pmatrix}
\boldsymbol\Psi                        &                                &                 & \ggb                           &                                      & -\mm                                                & \mm                                                 &                            &    \\
& \langle \mathbf{d} \rangle     &                 & -\mathbf{I}_m                  &                                      &                                                     &                                                     & -\mathbf{I}_m              &    \\
&                                &                 & -\boldsymbol{a}^T              &                                      &                                                     &                                                     &                            & -1 \\
\ggb^T                                 & -\mathbf{I}_m                  & -\boldsymbol{a} &                                & \mathbf{I}_m                         &                                                     &                                                     &                                 \\
&                                &                 & \langle \boldsymbol{s} \rangle & \langle \boldsymbol   \lambda\rangle &                                                     &                                                     &                            &    \\
\langle \boldsymbol\varepsilon \rangle &                                &                 &                                &                                      & \langle \boldsymbol\nu  - \boldsymbol\alpha \rangle &                                                     &                            &    \\
-\langle \boldsymbol\eta \rangle       &                                &                 &                                &                                      &                                                     & \langle \boldsymbol\beta -  \boldsymbol\nu  \rangle &                            &    \\
& \langle \boldsymbol\mu \rangle &                 &                                &                                      &                                                     &                                                     & \langle \mathbf{y} \rangle &    \\
&                                & \zeta           &                                &                                      &                                                     &                                                     &                            & z
\end{pmatrix}
\begin{pmatrix}
\boldsymbol\Delta \boldsymbol \nu         \\
\boldsymbol\Delta \boldsymbol y           \\
\Delta z                                  \\
\boldsymbol\Delta \boldsymbol \lambda     \\
\boldsymbol\Delta \boldsymbol s           \\
\boldsymbol\Delta \boldsymbol \varepsilon \\
\boldsymbol\Delta \boldsymbol \eta        \\
\boldsymbol\Delta \boldsymbol \mu         \\
\Delta \zeta
\end{pmatrix}
\begin{pmatrix}
-\boldsymbol{\delta_{\nu}}         \\
-\boldsymbol{\delta_y}             \\
-{\delta_z}                        \\
-\boldsymbol{\delta_{\lambda}}     \\
-\boldsymbol{\delta_{s}}           \\
-\boldsymbol{\delta_{\varepsilon}} \\
-\boldsymbol{\delta_{\eta}}        \\
-\boldsymbol{\delta_{\mu}}         \\
\end{pmatrix} \,,
\label{eq:hessian_matrix}
\end{equation}

where $\mathbf{I}_m$ is the $m\times m$ identity matrix.
The pointwise operations in Equations (<ref>), (<ref>) - (<ref>) are carried out on the vector components of their corresponding discretized fields.

As in <cit.>, the linear system of equations in (<ref>) can be solved for
$\boldsymbol{\Delta \upsilon} = (\bsdelta{\nu}, \bsdelta{y}, \Delta z, \bsdelta{\lambda}, \bsdelta{s}, \bsdelta{\varepsilon}, \bsdelta{\eta}, \bsdelta{\mu}, \Delta \zeta)$
quickly given that all the block matrices are diagonal, although this only happens for the element-wise uniform discretization.

We limit the Newton step from violating the inequalities in Equations (<ref>) as in <cit.>.
To do this, we first maximize the parameter $t$ such that

\begin{equation}
\begin{aligned}
t                                            & \leq 1 \label{eq:ls1}\,,                                                                                \\
\nu_l + t\Delta \nu_l - \alpha_l             & \geq 0.01\left(\nu_l - \alpha_l \right) ~ \text{for all  $l=1,...,n$} \,,                                                                      \\
\beta_l - \left(\nu_l + t\Delta \nu_l\right) & \geq 0.01\left(\beta_l - \nu_l  \right) ~ \text{for all  $l=1,...,n$} \,,                                                                      \\
\left( \mathbf{y}, z , \boldsymbol\lambda, \boldsymbol s,\boldsymbol\varepsilon, \boldsymbol\eta, \boldsymbol\mu, \zeta, \right)
+ t  \left( \boldsymbol\Delta \mathbf{y}, \Delta z , \boldsymbol\Delta \boldsymbol\lambda, \boldsymbol\Delta \boldsymbol s,\boldsymbol\Delta \boldsymbol\varepsilon, \boldsymbol\Delta \boldsymbol\eta, \boldsymbol\Delta \boldsymbol\mu, \Delta \zeta \right)
& \geq 0.01 \left( \mathbf{y}, z , \boldsymbol\lambda, \boldsymbol s, \boldsymbol\varepsilon, \boldsymbol\eta, \boldsymbol\mu, \zeta \right) \,.
\end{aligned}
\end{equation}

where all inequalities are performed on the vector components and $n$ is the number of mesh elements.
As in <cit.>, we assign $\boldsymbol\upsilon \leftarrow \boldsymbol\upsilon + \tau \bsdelta{\upsilon}$ where $\tau$ is the largest of $t,~t/2,~t/4,...$ such that

\begin{align}
\norm{F(\upsilon + \tau \Delta \upsilon)}_Q < \norm{F(\upsilon)}_Q \,.\label{eq:ls5}
\end{align}

We emphasize that $\norm{\cdot}_Q$ is not the ordinary vector norm, but the norm induced by the space $Q= V^* \times \mathbb{R}^m \times \mathbb{R}\times \mathbb{R}^m \times\mathbb{R}^m \times V \times V \times \mathbb{R}^m \times \mathbb{R} $, i.e.

\begin{equation}
\norm{F(\upsilon)}_Q = \sqrt{ \norm{\delta_{\nu}}_{V^*_h}^2 +
\norm{\delta_y}^2_{\mathbb{R}^m} +
\lVert\delta_z \rVert^2 +
\norm{\delta_{\lambda}}^2_{\mathbb{R}^m} +
\norm{\delta_s}^2_{\mathbb{R}^m} +
\norm{\delta_{\varepsilon}}_{V_h}^2 +
\norm{\delta_{\eta}}_{V_h}^2 +
\norm{\delta_{\mu}}^2_{\mathbb{R}^m} +
\lVert\delta_{\zeta}\rVert^2} \,,
\label{eq:norm_resi_mma}
\end{equation}

where, e.g.

\begin{equation}
\begin{aligned}
\norm{\delta_{\nu}}_{V_h^*}       & = \boldsymbol\delta_{\nu}^T \mathbf{M}^{-1} \boldsymbol\delta_{\nu} \,,             \\
\norm{\delta_{\varepsilon}}_{V_h} & = \boldsymbol\delta_{\varepsilon}^T \mathbf{M} \boldsymbol\delta_{\varepsilon}  \,.
\end{aligned}
\label{eq:norm_partial_resi}
\end{equation}

To recover the original NLP algorithm by <cit.> posed in the sequence space $\mathbb{R}^m$, we merely replace $\mathbf{M}$ with the identity matrix in the norm calculations and Equation (<ref>).

The discretized Fréchet derivative $\mathbf{D}\boldsymbol\theta$ needs to be passed to the NLP algorithm.
However, the GCMMA requires the discretized gradient $\boldsymbol{\nabla\theta}$ cf. Equations (<ref>), (<ref>) and (<ref>).a which is calculated using the Riesz map (<ref>) as

\begin{equation}
\begin{aligned}
\boldsymbol{\nabla \theta} & = \Phi_h^{-1}(\mathbf{D}\boldsymbol{\theta}) \,,    \\
& = \mathbf{M}^{-1} \mathbf{D}\boldsymbol{\theta} \,.
\end{aligned}
\end{equation}

The gradient components are used to evaluate Equations (<ref>) and (<ref>), i.e.

\begin{equation}
\begin{aligned}
p^{(k,j)}_{i,l} = (U_l^{(k)} - \nu_l^{(k,j)})^2\left(1.001\left(\nabla \theta_{i,l}(\nu^{(k,j)})\right)^+ +0.001\left(\nabla \theta_{i,l}(\nu^{(k,j)})\right)^- +\frac{\rho_i^{(k,j)}}{\nu_{\text{max}} - \nu_{\text{min}}}\right) \,, \\
q^{(k,j)}_{i,l} = (\nu_l^{(k)} - L_l^{(k,j)})^2\left(0.001\left(\nabla \theta_{i,l}(\nu^{(k,j)})\right)^+ +1.001\left(\nabla \theta_{i,l}(\nu^{(k,j)})\right)^- +\frac{\rho_i^{(k,j)}}{\nu_{\text{max}} - \nu_{\text{min}}}\right) \,.
\end{aligned}
\end{equation}

where the subscript $l$ in $p^{(k,j)}_{i,l}$ and $q^{(k,j)}_{i,l}$ corresponds to the mesh element $l$ and in $\nabla \theta_{i,l}$ it corresponds to the component of the gradient $\boldsymbol{\nabla \theta}_i$.

The termination criteria of the original problem (<ref>), i.e. not the subproblem (<ref>), is derived from its KKT conditions.
Notably, here we follow the convergence metric of <cit.> and monitor the gradient of the Lagrangian

\begin{align}
\Lagr(\nu, \lambda, \tau^+, \tau^-) & = \theta_0(\nu) + \sumonm \lambda_i \theta_i(\nu) +  \int_D \tau^-(\nu_{\text{min}} - \nu)~dV
+ \int_D \tau^+(\nu - \nu_{\text{max}})~dV \,.
\end{align}

Upon defining $\omega = (\nu, \lambda, \tau^-, \tau^+)$, the KKT conditions read

\begin{equation}
\begin{aligned}
\nabla_\nu    \Lagr(\omega)    & =\nabla \theta_0(\nu) + \sumonm \lambda_i \nabla \theta_i(\nu) - \tau^- + \tau^+ =0 \,,\phantom{i = 1, \ldots, m \,,}                                  \\
\theta_i(\nu)                  & \leq 0                                                                                & \; ~ i = 1, \ldots, m\,,\notag \\
\lambda_i \theta_i(\nu)        & = 0                                                                                & \; ~ i = 1, \ldots, m\,,\notag \\
\lambda_i                      & \geq 0                                                                                & \; ~ i = 1, \ldots, m\,,\notag \\
\nu_{\text{min}} \leq \nu      & \leq \nu_{\text{max}} \,,\notag                                                                                \\
\tau^-(\nu_{\text{min}} - \nu) & = 0 \,,                                                                                \\
\tau^+(\nu - \nu_{\text{max}}) & = 0 \,,                                                                                \\
\tau^+, \tau^-                 & \geq 0\notag \,.
\end{aligned}
\label{eq:kkt_stopping}
\end{equation}

Wherever the bound constraints are active, i.e. wherever $\nu = \nu_{\text{min}}$, or $\nu = \nu_{\text{max}}$, their corresponding Lagrange multipliers are $\tau^- =  \nabla \theta_0(\nu) + \sum_{i=1}^m \lambda_i \nabla \theta_i(\nu) $ or $\tau^+ = - \left( \nabla \theta_0(\nu) + \sum_{i=1}^m \lambda_i \nabla \theta_i(\nu) \right)$ respectively.
And since $\tau^- \geq 0 $, we have $\nabla \theta_0(\nu) + \sum_{i=1}^m \lambda_i \nabla \theta_i(\nu) \geq 0 $ and similarly since $\tau^+ \geq 0 $ we have $\nabla \theta_0(\nu) + \sum_{i=1}^m \lambda_i \nabla \theta_i(\nu) \leq 0 $.
Using these inequalities in the complementary slackness Equations (<ref>).f - (<ref>).g, transforms them to

\begin{equation}
\begin{aligned}
\lambda_i (\theta_i(\nu))^-                                                                           & = 0\,,\; i = 1, \ldots, m  \,,   \\
(\theta_i(\nu))^+                                                                                & = 0\,, \; i = 1, \ldots, m \,,   \\
\lambda_i                                                                                & \geq 0\,,\; i = 1, \ldots, m \,, \\
(\nu_{\text{min}} - \nu)\left( \nabla \theta_0(\nu) + \sumonm\lambda_i \nabla \theta_i(\nu) \right)^+ & = 0 \,,                          \\
(\nu_{\text{max}} - \nu)\left( \nabla \theta_0(\nu) + \sumonm\lambda_i \nabla \theta_i(\nu) \right)^- & = 0 \,,                          \\
\tau^+, \tau^-                                                                                & \geq 0
\end{aligned}
\label{eq:kkt_stopping_final}
\end{equation}

where the Lagrange multipliers $\lambda_i$ are obtained from the solution of the convex approximation problem (<ref>).
The norm of the KKT conditions is given by the norms of the left hand sides of Equations (<ref>).a-(<ref>).b and (<ref>).d-(<ref>).e where each norm is taken in its corresponding space, i.e., Equations (<ref>).a-(<ref>).b in $\mathbb{R}^m$ and Equations (<ref>).d-(<ref>).e in $V$.
The primal-dual NLP algorithm for the subproblem is outlined in Algorithm <ref>
MMA subproblem algorithm outline.
[1]
$\mathbf{Input:}$ Starting point for $k=1$: $\nu^{(k)} = \frac{\alpha + \beta}{2}, y^{(k)}_i = s^{(k)}_i = \lambda^{(k)}_i = 1.0, ~ \mu^{(k)}_i = \text{max}\{1.0, 0.5 c_i\}~ \text{for} ~ i=1,...,m, ~ z^{(k)}=\zeta^{(k)}=1.0, $
$\varepsilon^{(k)} = \text{max} \{ 1.0, \frac{1.0}{\nu^{(k)} - \alpha} \},\eta^{(k)} = \text{max}\{ 1.0, \frac{1.0}{\beta - \nu^{(k)}} \}  $ and $\epsilon = 1$.
$\epsilon > 10^{-5}$
Calculate $\boldsymbol{\Delta \upsilon}$ in Equation (<ref>).
Calculate step length $\tau$ in Equation (<ref>).
Let $\boldsymbol{\upsilon}^{(k+1)} = \boldsymbol{\upsilon}^{(k)} + \tau \boldsymbol{\Delta \upsilon}^{(k)}$.
$\norm{F(\upsilon)}_Q<0.9 \epsilon$
$\epsilon = 0.1 \epsilon$

As seen in the calculation of the search direction, cf. Equation (<ref>), the residual norm, cf. Equation (<ref>) and the convergence metric, cf. Equation (<ref>), the $L^2$ NLP algorithm cannot be obtained by merely scaling the derivatives that feed the $\mathbb{R}^n$ NLP algorithm by the mass matrix, e.g. be replacing $D\theta_i$ with $\nabla \boldsymbol \theta = \mathbf{M}^{-1} \boldsymbol{D\theta}$.
Obtaining a mesh independent NLP algorithm requires starting from the infinite dimensional formulation and using the tools presented in this paper to obtain the proper discretization.
Summarizing all the changes made to the GCMMA algorithm:

* The design variable $\nu$, its lower and upper bounds, the moving asymptotes $L$ and $U$ and the subproblem bounds $\alpha$ and $\beta$ are functions in the Hilbert space $L^2$.

* The gradients $\nabla \theta_i(\nu)$ for $i=0,...,m$ are used in the convex approximation, i.e. Equation (<ref>) instead of derivatives $D\theta_i(\nu)$.
We obtain the gradients $\nabla \theta_i(\nu)$ by applying the Riesz map as in Equation (<ref>).

* The convex approximation in Equation (<ref>) is built with an integral over the domain, instead of a summation over the vector of design variables.

* Similarly, the summations in the global convergence mechanism, i.e. Equations (<ref>) and (<ref>), are replaced with integrals over the domain.

* The norms used to check for convergence of convex approximation subproblem, i.e. Equation (<ref>) and the convergence metric, i.e. (<ref>) are taken in the appropriate spaces

Non-uniform mesh for the compliance problem.

Non-uniform mesh for the mechanism problem.

Non-uniform mesh for the stress constrained problem.

§.§ Adaptive mesh refinement

We apply AMR using an element-based error quantity that measures the jump of the volume fraction $\hat\nu$ ($\nu$ for the thermal problem) across each element $K$ edges:

\begin{align}
e_K = \frac{\sum_{\mathcal{E} \in K} \int_{\mathcal{E}} \llbracket \hat\nu \rrbracket ~dS}{\sum_{\mathcal{E} \in K} \int_{\mathcal{E}} ~dS}\,,
\end{align}

where $\mathcal{E}$ are the element $K$ edges.
Note that $e_K$ is properly weighted to account for the element size.

We perform either only refinement or coarsening during the AMR operation using Triangle
<cit.> and MeshPy <cit.>.
During refinement, we use MeshPy to update the element $\Omega_k$ area $a_K$ to

\begin{align}
a_K = \text{max} \left( \frac{a_K}{\left(1 + 100 \frac{e_K}{\underset{K\in \mathcal{T}}{\text{max}}e_K}\right)^2}, a_{\text{min}} \right) \,.
\end{align}

$\mathcal{T}$ is the mesh discretization.
The lower area limit $a_{\text{min}}$ prevents very small elements and is equal to 0.01 and $6\times 10^{-6}$ for the compliance and thermal-flow problem respectively.

The coarsening requires a different criteria because MeshPy does not have such functionality.
Instead, we remesh the domain from scratch, letting the largest element size to be
$E_{\mathcal{T}} \left(~\underset{K\in \mathcal{T}}{\text{max }} a_K \right)$,
where the expansion factor $E_{\mathcal{T}}$ controls the coarsening.
We use $E_{\mathcal{T}}=10.0$ and 2.0 for the compliance and thermal-flow problem respectively.
Since we want to keep the design geometry refined, we define the high error region

\begin{align}
D_R = \left\{ \mathbf{x} \in D ~|~ e_{K(\mathbf{x})} > e_{50}\right\}
\end{align}

where $e_{50}$ is the 50-th of the element error distribution.
This region is refined until the element areas are smaller than the smallest element
from the previous mesh.
$K(\mathbf{x})$ refers to those elements from the previous mesh that contain $\mathbf{x}$,

We remark here that the beam and thermal flow AMR examples start from the same respective mesh
and either refine it or coarsen it. This ensures that the initial mesh is fine enough
to approximate the infinite dimensional solution of the state variables.
An initially too coarse mesh will introduce a large error in the state variables
and the cost and constraint function derivatives which will affect the optimized design.

To ensure continuity of the optimization algorithm at iteration $k$,
the design variable $\nu_k$, the previous iterations $\nu_{k-1}$, $\nu_{k-2}$ and the MMA moving asymptotes
$L_{k}$ and $U_{k}$ are projected onto the refined or coarsened mesh.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Meshes for optimized designs for the AMR compliance problem with refinement.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Meshes for optimized designs for the AMR compliance problem with coarsening.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Meshes for optimized designs for the AMR thermal flow problem with refinement.

Optimization in $\mathbb{R}^n$     Optimization in $L^2$
90Strategy A

90Strategy B

Meshes for optimized designs for the AMR thermal flow problem with coarsening.