%~Mouliné par MaN_auto v.0.29.1 2023-10-18 10:55:22 \documentclass[CRMATH,Unicode,XML]{cedram} \TopicFR{Analyse numérique, Mécanique} \TopicEN{Numerical analysis, Mechanics} %\usepackage[colorlinks=true]{hyperref} \usepackage[nameinlink,capitalise]{cleveref} %\usepackage{gensymb} \usepackage{algorithm} \usepackage{algorithmic} %\usepackage{bm} \usepackage[abs]{overpic} \usepackage{amssymb} \usepackage{tikz} \usepackage{listings} \usepackage{mathtools} \crefname{equation}{}{} \crefname{figure}{Figure}{Figures} \crefname{remark}{Remark}{Remarks} \crefname{rema}{Remark}{Remarks} \newcommand*\modu[1]{{\bf #1}}%{{\textcolor{blue}{\bf #1}}} \makeatletter \def\@setafterauthor{% \vglue3mm% % \hspace*{0pt}% \begingroup\hsize=12cm\advance\hsize\abstractmarginL\raggedright \noindent %\hspace*{\abstractmarginL}\begin{minipage}[t]{10cm} \leftskip\abstractmarginL \normalfont\Small \@afterauthor\par \endgroup \vskip2pt plus 3pt minus 1pt } \makeatother \usetikzlibrary{calc,positioning} \definecolor{mygray}{rgb}{1.0,0.95,0.90} \definecolor{truegray}{rgb}{0.9,0.9,0.9} \definecolor{mygreen}{rgb}{0,0.6,0} \definecolor{mygreen2}{rgb}{0,0.8,0} \definecolor{mypink}{rgb}{1.0,0.56,0.81} \lstdefinestyle{ff} { language = C++, basicstyle=\ttfamily, backgroundcolor=\color{truegray}, frame=tlbr,framesep=4pt,framerule=1pt, captionpos=b, commentstyle=\color{mygreen}, emph=[1] { macro, border, fespace, }, emphstyle=[1]{\bf\color{blue}}, emph=[2] { varf, matrix, real, string, int, Vh, Vh0, Vh0i, Vh2, Vhi, mesh, }, emphstyle=[2]{\bf\color{black}}, emph=[3] { int2d, int1d, on, readmesh, savemesh, load, include, trunc, }, emphstyle=[3]{\bf\color{violet}}, emph=[4] { EigenValue, buildmesh, adaptmesh, movemesh, checkmovemesh, solve, problem, laplace, elas, }, emphstyle=[4]{\color{red} \bfseries}, emph=[5] { volfrac, readsol, advectRedist, normalvec, curvature }, emphstyle=[5]{\color{mygreen2} \bfseries}, } \lstdefinestyle{shell} { basicstyle=\ttfamily\color{white}, frame=tlbr,framesep=4pt,framerule=1pt, backgroundcolor=\color{black}, captionpos=b, emph=[1] { python, }, emphstyle=[1]{\color{red}\bfseries}, } \lstdefinelanguage{mypython}{ basicstyle=\ttfamily\color{black}, morekeywords = [1]{print,format}, morekeywords = [2]{if,else,and,or,break,continue,def,return}, keywordstyle = [1]\bf\color{lightgray}, keywordstyle = [2]{\bf\color{mypink}}, morestring=[b]", stringstyle = {\color{lightgray}}, morecomment=[l]{\#}, } \lstdefinestyle{python} { language = mypython, commentstyle=\color{mygreen}, frame=tlbr,framesep=4pt,framerule=1pt, showstringspaces=false, backgroundcolor=\color{mygray}, captionpos=b, emph=[1] { mechtools, mshtools, inigeom, inout, lstools, path, subprocess, }, emphstyle=[1]{\bf\color{blue}}, } \DeclareMathOperator{\dv}{div} \DeclareMathOperator{\tr}{tr} \DeclareMathOperator{\I}{I} \DeclareMathOperator{\Id}{Id} \DeclareMathOperator{\Per}{Per} \DeclareMathOperator{\Vol}{Vol} \DeclareMathOperator{\dVol}{dVol} \DeclareMathOperator{\lift}{Lift} \DeclareMathOperator{\drag}{Drag} \DeclareMathOperator{\Reop}{Re} \newcommand{\orm}{\mathrm{o}} \newcommand{\dd}{\mathrm{d}} \newcommand*{\dx}{\dd x} \newcommand*{\ds}{\dd s} \newcommand*{\dt}{\dd t} \newcommand*{\dJ}{\dd J} \newcommand*{\dC}{\dd C} \newcommand{\e}{\varepsilon} \newcommand{\calQ}{{\mathcal Q}} \newcommand{\calS}{{\mathcal S}} \newcommand{\calT}{{\mathcal T}} \newcommand{\calTt}{{\mathcal T}_{\mathrm{temp}}} \newcommand{\caliT}{{\mathcal T}_{\mathrm{int}}} \newcommand{\caleT}{{\mathcal T}_{\mathrm{ext}}} \newcommand{\Winfty}{W^{1,\infty}(\mathbb{R}^d,\mathbb{R}^d)} \newcommand{\R}{{\mathbb R}} \newcommand{\CD}[1]{#1}%{{\color{blue}{#1}}} \newcommand{\D}{\ensuremath \,\dd} \newcommand{\In}{\text{ in }} \newcommand{\On}{\text{ on }} \newcommand{\n}{{\bm n}} \newcommand{\partialn}[1]{\frac{\partial #1}{\partial n}} \graphicspath{{./figures/}} \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} \hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red} \newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}} \newcommand*{\relabelitem}{\renewcommand{\labelenumi}{(\theenumi)$\mkern -18mu$}} \newcommand*{\relabeliitem}{\renewcommand{\labelenumii}{(\theenumii)$\mkern -18mu$}} \newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel} \newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel} \newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel} \newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel} \let\oldtilde\tilde \renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}} \let\oldhat\hat \renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}} \let\oldforall\forall \renewcommand*{\forall}{\mathrel{\oldforall}} \title{Shape optimization using a level set based mesh evolution method: an overview and tutorial} \alttitle{Une strat\'egie d'\'evolution de maillage bas\'ee sur la m\'ethode des lignes de niveaux pour l'optimisation de formes : survol et mise en pratique} \author{\firstname{Charles} \lastname{Dapogny}} \address{Univ. Grenoble Alpes, CNRS, Grenoble INP (Institute of Engineering Univ. Grenoble Alpes), LJK, 38000 Grenoble, France} \email{charles.dapogny@univ-grenoble-alpes.fr} \author{\firstname{Florian} \lastname{Feppon}} \address{Department of Computer Science, KU Leuven, Celestijnenlaan 200A, 3001 Heverlee, Belgium} \email{florian.feppon@kuleuven.be} \thanks{The work of C.D. is partially supported by the project ANR-18-CE40-0013 SHAPO, financed by the French Agence Nationale de la Recherche (ANR).} \subjclass{49M41, 65K05, 65N50, 74P05, 74P10, 90C90} \begin{abstract} This article revolves around a recent numerical framework for shape and topology optimization, which features an exact mesh of the shape at each iteration of the process, while still leaving the room for an arbitrary evolution of the latter (including changes in its topology). In a nutshell, two complementary representations of the shape are combined: on the one hand, it is meshed exactly, which allows for precise mechanical calculations based on the finite element method; on the other hand, it is described implicitly, using the level set method, which makes it possible to track its evolution in a robust way. In the first part of this work, we overview the main aspects of this numerical strategy. After a brief presentation of some necessary background material -- related to shape optimization and meshing, among others -- we describe the numerical schemes involved, notably when it comes to the practice of the level set method, the remeshing algorithms, and the considered optimization solver. This strategy is illustrated with 2d and 3d numerical examples in various physical contexts. In the second part of this article, we propose a simple albeit efficient \texttt{python}-based implementation of this framework. The code is described with a fair amount of details, and it is expected that the reader can easily elaborate upon the presented examples to tackle his own problems. \end{abstract} \begin{altabstract} Cet article traite d'un cadre de travail r\'ecent d\'edi\'e \`a la r\'esolution num\'erique de probl\`emes d'optimisation de formes\,; il s'illustre par une repr\'esentation exacte, maill\'ee, de la forme \`a chaque it\'eration du proc\'ed\'e, tout en laissant la place \`a une \'evolution arbitraire de celle-ci (y compris des changements de sa topologie). L'id\'ee centrale de cette strat\'egie est de combiner deux repr\'esentations compl\'ementaires de la forme\,: d'une part, celle-ci est maill\'ee explicitement, de sorte qu'il est possible d'effectuer des calculs m\'ecaniques pr\'ecis par la m\'ethode des \'el\'ements finis\,; d'autre part, elle est d\'ecrite implicitement, par la m\'ethode des lignes de niveaux, facilitant ainsi le suivi robuste de son \'evolution. Dans la premi\`ere partie de ce travail, on r\'esume les points saillants de cette strat\'egie num\'erique. Apr\`es avoir bri\`evement rappel\'e quelques notions de base -- en lien, entre autres, avec l'optimisation de formes et le maillage -- on d\'ecrit les sch\'emas num\'eriques mis en jeu, notamment pour la pratique de la m\'ethode des lignes de niveaux, les algorithmes de remaillage, et l'algorithme d'optimisation num\'erique. Cette m\'ethodologie est illustr\'ee par des exemples num\'eriques en deux et trois dimensions d'espace, dans diff\'erents contextes physiques. Dans la seconde partie de cet article, on propose une impl\'ementation \texttt{python} open-source, simple mais efficace, de ce cadre de travail. Le code est d\'etaill\'e de sorte que le lecteur puisse facilement intervenir dedans et le modifier pour traiter un probl\`eme de son choix. \end{altabstract} \begin{document} \maketitle \section{Introduction} The scarcity and the cost of resources such as raw materials and energy make it ever more necessary to tailor the design of physical devices from the early stages of design, so that they achieve their purpose with the minimum amount of constituent material and input power. For this reason, shape and topology optimization has aroused a tremendous enthusiasm in both academic and industrial communities, as a set of fully automated algorithms for predicting optimized designs with respect to physical requirements, see e.g.~\cite{allaire2007conception,azegami2020shape,bendsoe2013topology,eschenauer2001topology,mohammadi2010applied}. From the mathematical point of view, shape and topology optimization problems show up under the generic form \[ \min_\Omega J(\Omega) \text{ s.t. } \begin{cases} G(\Omega) = 0,\\ H(\Omega) \leq 0, \end{cases} \] where the objective and the equality and inequality constraint functionals $J(\Omega)$, $G(\Omega)$ and $H(\Omega)$ depend on the optimized shape $\Omega$ in a possibly very intricate way: in concrete applications, they involve a taste of the physical behavior of $\Omega$ through the solution to a partial differential equation posed on $\Omega$. For instance, \begin{itemize} \item When $\Omega$ accounts for a mechanical structure subjected to prescribed loads, $J(\Omega)$, $G(\Omega)$ and $H(\Omega)$ bring into play the elastic displacement, characterized as the solution to the linear elasticity system posed on $\Omega$; \item When $\Omega$ represents a fluid duct, $J(\Omega)$, $G(\Omega)$ and $H(\Omega)$ depend on the velocity of the fluid, which is governed by the Stokes or the Navier--Stokes equations on $\Omega$. \end{itemize} Leaving aside theoretical issues, regarding for instance the existence of an optimal solution (a purpose for which we refer to e.g.~\cite{bucur2002variational,henrot2018shape}), the numerical treatment of such problems is classically plagued by the need to reconcile two antagonistic needs. On the one hand, the evaluation of the objective and constraint functionals $J(\Omega)$, $G(\Omega)$, $H(\Omega)$ and the calculation of their sensitivities with respect to ``small'' variations of the design $\Omega$ require to solve one or several partial differential equations posed on $\Omega$. This task is typically accomplished by using the finite element method, which in turn relies on a high-quality mesh $\calT$ of $\Omega$. On the other hand, the update of the shape $\Omega$ between successive stages of the optimization process, which is typically realized by deforming $\Omega$ according to a velocity field $\theta : \R^d \to \R^d$, is a difficult operation to translate in terms of a mesh of $\Omega$, since the latter very likely becomes degenerate or even invalid in the course of such process. This difficulty of finding a numerical representation of the shape which is amenable to all the operations of a shape and topology optimization workflow has been acknowledged as a genuine bottleneck since the early days of shape optimization. It is also encountered in multiple disciplines such as inverse problems, image processing, etc., see for instance Chap. 23 in~\cite{frey2007mesh} about this point. Various paradigms have been thought of to overcome this issue; let us notably mention two of them: \begin{itemize} \item Relaxation methods replace the parametrization of the shape and topology optimization problem by classical ``black-and-white'' designs $\Omega$ -- or equivalently their characteristic function, taking the values $0$ and $1$ outside and inside $\Omega$ respectively -- by ``grayscale'' density functions $h: D \to [0,1]$ defined on a large ``hold-all'' domain $D$ (and possibly a microstructure tensor, describing the local microscopic structure of the shape). This change of perspectives can be rigorously justified by the mathematical theory of homogenization~\cite{allaire2002shape, allaire1997shape, kohn1986optimal}, which has inspired simplified, heuristic variants of ``classical'' shape optimization problems such as the popular SIMP framework in mechanical engineering, see~\cite{bendsoe1988generating,bendsoe2013topology}. \item Eulerian methods, such as the level set method~\cite{allaire2004structural,osher2001level,sethian2000structural,wang2003level} or the phase-field method~\cite{blank2012phase,bourdin_chambolle_2003,dede2012isogeometric,takezawa2010shape} bring into play a fixed mesh of a large ``hold-all'' domain $D$. The sought shape $\Omega$ is defined implicitly, in terms of quantities defined on the whole domain $D$. For instance, the level set method features shapes defined as negative subdomains $\left\{ x \in D,\:\: \phi(x) < 0 \right\}$ of functions $\phi : D \to \R$, see \cref{sec.lsgen} below for more details. \end{itemize} The aforementioned methods alleviate the need to track the evolution of a mesh of the optimized shape $\Omega$, but this comes at a price: the fact that no proper mesh of $\Omega$ is available raises the need for an approximation of the partial differential equations characterizing its physical behavior, or for an advanced and often intrusive finite element method tailored for implicitly-defined domains, such as the eXtended Finite Element Method (X-FEM) \cite{duysinx2006generalized}, or the cut-FEM method~\cite{burman2018shape}. In this article, we overview a recent mesh evolution framework introduced in~\cite{allaire2011topology,allaire2013mesh,allaire2014shape}, which allows to describe arbitrarily large deformations of the shape throughout the process, and which at the same time features an exact, well-shaped mesh of the latter. This framework leverages two complementary representations of shapes: on the one hand, they are meshed exactly, and so accurate mechanical computations can be conducted by standard finite element methods and solvers used in a ``black-box'' fashion. On the other hand, they are represented by means of the level set method on a larger, fixed computational domain, which makes it possible to account for arbitrarily large deformations. The core of this strategy is a set of numerical algorithms allowing to switch from one of these representations to the other, and so to consistently use that which is most relevant for every operation involved in the shape optimization workflow. Since its introduction in the context of structural mechanics, this idea has been successfully applied in a variety of more challenging physical contexts, such as fluid mechanics and fluid-structure interaction~\cite{feppon2018shape,feppon2020topology}, heat transfer~\cite{feppon2021body}, quantum chemistry~\cite{braida2022shape}, plasticity~\cite{desai2021topology}, fracture mechanics~\cite{desai2022topology}, etc. The aim of this article is twofold. First and foremost, we describe this level set based mesh evolution strategy and its main numerical ingredients, such as the algorithms involved in the practice of the level set method, the meshing aspects of the framework, and the efficient resolution of infinite-dimensional constrained optimization problems. We conclude this presentation with a selection of numerical results obtained by application of this strategy. The second purpose of this work is in line with the multiple educational articles that have been devoted to shape and topology optimization, see e.g.~\cite{allaire2006structural,laurain2018level,sigmund200199} or~\cite{wang2021comprehensive} for an exhaustive list. We propose and detail a \texttt{python} implementation of this level set based mesh evolution framework which is at the same time pedagogical (and thus reasonably simple), and general enough to allow the user to easily build upon this code so as to address his personal shape optimization problems. The remainder of this article is organized as follows. In \cref{sec.th}, we present the model physical setting of our discussion, that of shape optimization of elastic structures. We introduce the chief theoretical concepts at stake, and provide a general sketch of a typical shape optimization procedure. In \cref{sec.meshissue}, we discuss the choice of an adequate representation of the shape in the perspective of accounting for its update between the various iterations of the shape optimization process. After a brief reminder of the relevant notions, we describe Lagrangian strategies, their limitations, and we present the level set method and its applications in shape and topology optimization as a potential remedy. In \cref{sec.bf}, we describe the level set based mesh evolution method at stake in this article, and we describe the main numerical operations involved. In \cref{sec.num}, we then present several numerical results obtained with this strategy. This article ends with a fairly long appendix where a proposed open-source implementation of this framework is presented: each step of the method is carefully described, with the hope that it is easy for the reader to get into the code and elaborate upon it for its own purposes. \section{Presentation of the shape optimization framework}\label{sec.th} In this section, we introduce the theoretical framework of this article, and we discuss a few important notions. To set ideas, our discussion unfolds in the relatively simple physical setting of linearly elastic structures, which is presented in \cref{sec.sopb}; some basic facts about Hadamard's boundary variation method and shape derivatives are then recalled in \cref{sec.hadamard}. Finally, we sketch a generic shape optimization algorithm in \cref{sec.tonum}, as a support for the subsequent practical considerations. \subsection{A model shape optimization problem in the context of structural mechanics}\label{sec.sopb} Throughout this presentation, a shape is a bounded, Lipschitz domain $\Omega \subset \R^d$ ($d=2,3$), whose boundary is the reunion of three disjoint, open regions $\Gamma_D$, $\Gamma_N$, $\Gamma$: \[ \partial \Omega = \overline{\Gamma_D} \cup \overline{\Gamma_N} \cup \overline{\Gamma}. \] In this decomposition, \begin{itemize} \item $\Gamma_D$ is a given piece of hypersurface where $\Omega$ is clamped; \item $\Gamma_N$ is another given piece of hypersurface where loads $g \in L^2(\Gamma_N)^d$ are applied; \item The remaining part $\Gamma$ of $\partial\Omega$ is traction free; it is the only part of $\partial\Omega$ which is subject to optimization; \end{itemize} see \cref{fig.setso} for an illustration. \begin{figure}[!ht] \centering \includegraphics[width=0.6\textwidth]{./figures/setso} \caption{Physical setting of mechanical structures considered in \cref{sec.sopb}.}\label{fig.setso} \end{figure} Omitting body forces for simplicity, the displacement of $\Omega$ is the unique solution $u_\Omega$ in the~space \begin{equation*}\label{eq.defH1GD} H^1_{\Gamma_D}(\Omega)^d := \left\{ u \in H^1(\Omega)^d, \:\: u =0 \text{ on } \Gamma_D \right\} \end{equation*} to the linearized elasticity system: \begin{equation}\label{eq.elas} \left\{ \begin{aligned} -\dv(Ae(u_\Omega)) &= 0 && \text{in } \Omega, \\ u_\Omega &= 0 && \text{on } \Gamma_D,\\ Ae(u_\Omega) n &= g && \text{on }\Gamma_N, \\ Ae(u_\Omega) n &= 0 && \text{on }\Gamma. \end{aligned}\right. \end{equation} In this formulation, $H^1(\Omega)$ is the usual Sobolev space of functions in $L^2(\Omega)$ whose first-order derivatives are also in $L^2(\Omega)$, see~\cite{adams2003sobolev}. We have denoted by $n : \partial \Omega \to \R^d$ the unit normal vector to $\partial \Omega$, pointing outward $\Omega$, and the symmetric $d \times d$ matrix $e(u) := \frac12(\nabla u + \nabla u^T)$ is the strain tensor associated to a displacement field $u : \Omega \to \R^d$. The Hooke's law $A$ is a fourth-order tensor characterizing the physical properties of the constituent material of $\Omega$: it relates the state of stress $\sigma(u)$ inside the structure with the strain $e(u)$ via the relation \[ \sigma(u) = Ae(u), \text{ where, for any symmetric }d\times d\text{ matrix } e, \quad A e = 2\mu e + \lambda \tr(e) \I, \] and $\lambda$, $\mu$ are the Lam\'e coefficients of the material. In this context, a generic shape optimization problem reads \begin{equation}\label{eq.sopb} \min_\Omega J(\Omega) \text{ s.t. } \begin{cases} G(\Omega) = 0, \\ H(\Omega) \leq 0. \end{cases} \end{equation} Several important examples about the objective functional $J(\Omega)$ are: \begin{itemize} \item The compliance of $\Omega$ \begin{equation}\label{eq.compliance} C(\Omega) = \int_\Omega Ae(u_\Omega):e(u_\Omega) \:\dx = \int_\Gamma g \cdot u_\Omega \:\ds, \end{equation} where $:$ is the Frobenius inner product over the set of $d \times d$ matrices. Equivalently, $C(\Omega)$ measures the total elastic energy stored inside $\Omega$ or the work done by the applied loads $g$; \item A least-square discrepancy criterion \[ D(\Omega) = \int_\Omega k(x) |u_\Omega - u_T(x) |^2 \:\dx \] between the elastic displacement of $\Omega$ and a target displacement $u_T(x)$, weighted by a function $k(x)$. \item An integral measure of the stress inside the structure \begin{equation}\label{eq.stressfunc} S(\Omega) = \int_\Omega k(x) \| \sigma(u_\Omega) \|^2 \:\dx, \end{equation} where for any $d \times d$ matrix $\sigma$, we have denoted $\| \sigma \|^2 := \sigma : \sigma$. \end{itemize} In the formulation \cref{eq.sopb}, $G(\Omega) = (G_i(\Omega))_{i=1,\ldots,p}$ and $H(\Omega) = (H_j(\Omega))_{j=1,\ldots,q}$ are collections of $p$ real-valued equality constraints and $q$ inequality constraints, respectively. These functionals may be, for instance: \begin{itemize} \item The volume $\Vol(\Omega)$ or the perimeter $\Per(\Omega)$ of the shape, \begin{equation}\label{eq.volume} \Vol(\Omega) = \int_\Omega \:\dx, \quad\text{and}\quad \Per(\Omega) = \int_{\partial\Omega} \:\ds. \end{equation} \item Other constraints related to the geometry of $\Omega$ (thickness, curvature radii, etc.), as imposed by the manufacturing process. \end{itemize} In the sequel, when it allows to simplify the exposition, we shall sometimes consider unconstrained versions of \cref{eq.sopb}, of the form \begin{equation}\label{eq.sopbuc} \min_\Omega J(\Omega), \end{equation} where the minimized function $J(\Omega)$ may represent e.g. a weighted sum of some of the above shape functionals. \subsection{The boundary variation method of Hadamard}\label{sec.hadamard} Most numerical algorithms dedicated to problems of the form \cref{eq.sopb} rely on the derivatives of the objective and constraint functions with respect to the optimized variable. In the present context, this raises the need to account for derivatives with respect to the domain. This task can be achieved in different fashions, and we rely on the Hadamard's boundary variation method, pioneered by~\cite{hadamard1908memoire,murat1976controle}, see also~\cite{allaire2007conception,delfour2011shapes,henrot2018shape,sokolowski1992introduction}. In this framework, variations of a given shape $\Omega$ are considered under the form \begin{equation}\label{eq.varhadamard} \Omega_\theta := (\Id + \theta)(\Omega), \:\: \theta \in \Winfty, \:\: \| \theta \|_{\Winfty} < 1, \end{equation} where $\Winfty$ is the Sobolev space of Lipschitz vector fields on $\R^d$, see~\cite{evans2015measure}. Intuitively, \cref{eq.varhadamard} expresses the fact that $\Omega$ is deformed according to the ``small'' vector field $\theta$, see \cref{fig.hadamard}. \begin{figure}[!ht] \centering \includegraphics[width=0.6\textwidth]{./figures/hadamard} \caption{Deformed version $\Omega_\theta$ of a shape $\Omega$ according to Hadamard's boundary variation method.}\label{fig.hadamard} \end{figure} \begin{rema} Often in practice, the vector fields $\theta$ featured in Hadamard's method are required to enjoy higher regularity, or to vanish on a region of space which is not subject to modifications, so that they are actually confined to a subset of $\Winfty$. To keep the presentation elementary, we ignore this technicality in the following, see also \cref{sec.descentth}. \end{rema} One functional of the domain $F(\Omega)$ is then said to be shape differentiable at a particular shape $\Omega$ if the underlying mapping $\theta \mapsto F(\Omega_\theta)$, from $\Winfty$ into $\R$ is Fr\'echet differentiable at $\theta=0$. The corresponding derivative $\theta \mapsto F^\prime(\Omega)(\theta)$ satisfies the following expansion: \begin{equation}\label{eq.variationtheta} F(\Omega_\theta) = F(\Omega) + F^\prime(\Omega)(\theta) + \orm(\theta), \text{ where } \frac{|\orm(\theta)|}{\|\theta\|_{\Winfty}} \xrightarrow{\theta\to 0} 0. \end{equation} The shape derivatives of the objective and constraint functions $J(\Omega)$, $G(\Omega)$ and $H(\Omega)$ of an optimization problem of the form \cref{eq.sopb} are handful for a variety of purposes. From the theoretical point of view, they are the building blocks of the necessary conditions for a shape $\Omega$ to be locally optimal with respect to \cref{eq.sopb}; from the numerical vantage, they make it possible to calculate descent directions as vector fields $\theta : \R^d \to \R^d$ encoding deformations of $\Omega$ improving its ``performance''. Like those introduced in \cref{sec.sopb}, the functionals of interest in concrete applications usually depend on the shape $\Omega$ in a quite intricate way, via a ``state'' $u_\Omega$, characterized as the solution to a ``physical'' partial differential equation posed on $\Omega$. Nevertheless, their shape derivatives can be calculated thanks to the adjoint method, pertaining to the more general field of optimal control. This stake is conceptual and by no means trivial; however, it is well-understood in the literature, and we do not expand on the subject, referring to~\cite{lions1971optimal}, or~\cite{plessix2006review} for a comprehensive introduction to this method; see also the recent review~\cite{allaire2020survey} in the more specific context of shape optimization. Let us simply recall that the shape derivatives of such functionals usually involve the function $u_\Omega$ as well as an ``adjoint state'' $p_\Omega$, satisfying a partial differential equation very similar to that for $u_\Omega$, with a different right-hand side. Here are a few examples, working under mild regularity assumptions on the shape $\Omega$ or the state $u_\Omega$, which are omitted for brevity: \begin{itemize} \item The volume $\Vol(\Omega)$ has the shape derivative \[ \Vol^\prime(\Omega)(\theta) = \int_{\partial \Omega} \theta \cdot n \:\ds. \] \item The shape derivative of the compliance $C(\Omega)$ given by \cref{eq.compliance} reads \[ C^\prime(\Omega)(\theta) = - \int_{\Gamma} Ae(u_\Omega): e(u_\Omega) \:\theta\cdot n \:\ds. \] \item The stress functional $S(\Omega)$ in \cref{eq.stressfunc} has the shape derivative \begin{equation}\label{eq.Sprime} S^\prime(\Omega)(\theta) = \int_\Gamma \Big(k(x) \| \sigma(u_\Omega) \|^2 + Ae(u_\Omega):e(p_\Omega) \Big) \:\theta\cdot n \: \ds, \end{equation} where the adjoint state $p_\Omega$ is defined as the solution in $H^1_{\Gamma_D}(\Omega)^d$ to the equation \begin{equation}\label{eq.adjstress} \left\{ \begin{aligned} -\dv (Ae(p_\Omega)) &= \dv(k(x) A \sigma(u_\Omega)) && \text{in } \Omega, \\ p_\Omega &= 0 && \text{on } \Gamma_D, \\ Ae(p_\Omega)n &= 0 && \text{on } \Gamma \cup \Gamma_N. \end{aligned} \right. \end{equation} \end{itemize} The above expressions exemplify a quite general phenomenon. Shape derivatives are naturally calculated in volume form, i.e. their expressions are made of volume integrals on $\Omega$, involving $u_\Omega$, $\theta$, $\nabla \theta$, etc.; under suitable regularity assumptions on the shape $\Omega$, the state $u_\Omega$ and the possible adjoint state $p_\Omega$, they can often be given a surface expression of the type: \begin{equation}\label{eq.surfsd} J^\prime(\Omega)(\theta) = \int_\Gamma v_\Omega \:\theta\cdot n \:\ds, \end{equation} where $v_\Omega : \Gamma \to \R$ is a scalar field whose expression depends on $J(\Omega)$, $u_\Omega$ and $p_\Omega$. In particular, such a shape derivative $J^\prime(\Omega)(\theta)$ only depends on the values of the normal component $\theta\cdot n$ of the deformation $\theta$ on the boundary $\partial \Omega$, in agreement with the so-called \emph{Structure theorem} for shape derivatives, see e.g.~\cite[Section~5.9]{henrot2018shape} or~\cite[Chapter~9, Section~3.4]{delfour2011shapes}. Surface expressions \cref{eq.surfsd} for shape derivatives are convenient for a variety of purposes; for instance, when the unconstrained minimization problem \cref{eq.sopbuc} of $J(\Omega)$ is considered, a descent direction $\theta$ is easily obtained by imposing that \begin{equation}\label{eq.thetavn} \theta = - v_\Omega n \text{ on } \Gamma, \text{ so that } J^\prime(\Omega)(\theta) = -\int_\Gamma v_\Omega ^2 \:\ds < 0. \end{equation} On the other hand, however more difficult to exploit, volume expressions may lend themselves to more accurate numerical discretization, see~\cite{hiptmair2015comparison}. \begin{rema}\label{rem.dertopo} \CD{Variations of the domain of a different nature from \cref{eq.varhadamard} are possible, leading to as many different notions of ``derivative with respect to the domain'', or, more accurately, of asymptotic expansion with respect to perturbations of the domain}. Notably, it is possible to account for variations of a shape $\Omega$ of the form $\Omega_{x_0,r} := \Omega \setminus \overline{B(x_0,r)}$, where $x_0$ is a given point inside $\Omega$ and $r >0$. This leads to expansions of the form \[ J(\Omega_{x_0,r}) = J(\Omega) + r^d dJ_T(x_0) + \orm(r^d), \] where $dJ_T(x_0)$ is called the topological derivative of $J$ at $x_0$ and measures the sensitivity of $J$ with respect to the nucleation of an infinitesimally small hole inside $\Omega$. We refer to~\cite{novotny2012topological} about this concept, and to~\cite{allaire2005structural,amstutz2006new} about different ways of using it in the context of shape and topology optimization; see also \cref{sec.2gofurther} for an implementation. Let us also evoke the existence of more exotic ``topological ligament expansions'', evaluating the sensitivity of a shape with respect to the addition of a thin ligament, see~\cite{dapogny2020topolig,kobayashi2021cellular,nazarov2004topological}. \end{rema} \subsection{An abstract shape optimization algorithm}\label{sec.tonum} The notion of shape derivative underlies a wide range of algorithms for dealing with shape optimization problems of the form \cref{eq.sopb}; in broad outline, each iteration in their implementation can be decomposed into three main stages: \begin{enumerate} \item The physical quantities attached to the current shape $\Omega^n$ (the elastic displacement $u_{\Omega^n}$ and the adjoint state $p_{\Omega^n}$ in the context of \cref{sec.sopb}) are computed by solving the corresponding partial differential equations. For instance, this can be conveniently realized thanks to a finite element solver, if a mesh of the shape $\Omega^n$ is available. \item A decent direction $\theta^n$ for the problem \cref{eq.sopb} is inferred: $\theta^n$ is a vector field such that the deformed version $(\Id + \tau^n \theta^n)(\Omega^n)$ of $\Omega^n$ along $\theta^n$ for a ``small enough'' time step $\tau^n$ has ``improved'' performance over $\Omega^n$. For instance, when the constraint-free problem~\cref{eq.sopbuc} is considered, $\theta^n$ is simply a vector field such that $J^\prime(\Omega^n)(\theta^n) < 0$. When a more general constrained optimization problem of the form \cref{eq.sopb} is considered, $\theta^n$ is calculated from the knowledge of the shape derivatives $J^\prime(\Omega)$, $G^\prime(\Omega)$, $H^\prime(\Omega)$ owing to a constrained optimization algorithm. \item The shape $\Omega^n$ is updated into $\Omega^{n+1} := (\Id + \tau^n \theta^n)(\Omega^n)$. \end{enumerate} This generic procedure is summarized in \cref{algo.sketchsg}. \begin{algorithm}[ht] \caption{Generic shape gradient algorithm.}\label{algo.sketchsg} \begin{algorithmic}[0] \STATE \textbf{Initialization:} Initial shape $\Omega^0$. \FOR{$n=0,\dots,$ until convergence} \STATE \begin{enumerate}\relabelitem \item \label{1.1} Calculate the solution $u_{\Omega^n}$ (resp. $p_{\Omega^n}$) to the state (resp. adjoint) equation posed in~$\Omega^n$. \item \label{1.2} From the theoretical formulas for the shape derivatives $J^\prime(\Omega)(\theta)$, $G^\prime(\Omega)(\theta)$ and $H^\prime(\Omega)(\theta)$, infer a descent direction $\theta^n$ from $\Omega^n$ for the shape optimization problem~\cref{eq.sopb}. \item \label{1.3} Deform $\Omega^n$ according to $\theta^n$ for a small descent step $\tau^n>0$, so that the new shape \[ \Omega^{n+1} := (\Id + \tau^n \theta^n) (\Omega^n) \] is ``better'' than the previous one in view of \cref{eq.sopb}. \end{enumerate} \ENDFOR \RETURN $\Omega^n$ \end{algorithmic} \end{algorithm} One critical issue lies in the difficulty of finding a numerical description of the shape and its evolution which is appropriate to each of these three stages. For instance, choosing to represent the shape with a computational mesh conveniently allows to carry out the mechanical analyses implied by the first stage. Unfortunately, such decision makes it notoriously difficult to realize the update of the shape in the third stage in a robust way, all the more so as when large deformations (not to say topological changes) are at stake. This dilemma is faced by all implementations of \cref{algo.sketchsg}, and we focus specifically on this point from the next \cref{sec.meshissue}. \section{Meshing aspects of shape optimization}\label{sec.meshissue} The present section is dedicated to the numerical representation of the shape in the perspective of realizing its update between two successive steps of a shape optimization procedure. At first, we briefly recall a few basic definitions and important facts about meshing in \cref{sec.mesh}, referring for instance to the books~\cite{borouchaki2017meshing,botsch2010polygon,frey2007mesh} for more in-depth presentations. We then present the early ``Lagrangian'' strategies for shape optimization, emphasizing on their inherent limitations related with meshing aspects. Finally, we describe the level set method as one attempt to circumvent them, and as a cornerstone of the numerical framework discussed in this article. \subsection{Basic notions about meshes}\label{sec.mesh} Let $\Omega$ be a polygonal domain; a simplicial mesh of $\Omega$ is a collection $\calT=\left\{ T_i \right\}_{i=1,\ldots,N}$ of open simplices (i.e. triangles in 2d, tetrahedra in 3d) which constitute a covering of $\Omega$, that is: \[ \overline \Omega = \bigcup_{i=1}^{N} \overline{T_i}. \] In addition, one usually demands that \begin{itemize} \item $\calT$ should be \emph{valid}, in the sense that the $T_i$ are mutually disjoint: $T_i \cap T_j = \emptyset$ when $i \neq j$. \item $\calT$ should be \emph{conforming}: for $i \neq j$ the intersection $\overline{T_i} \cap \overline{T_j}$ is either a vertex, an edge, or (in 3d) a face of $\calT$. \end{itemize} The volume mesh $\calT$ of $\Omega$ additionally bears the information of a \emph{surface mesh} $\calS_{\calT}$, that is, a mesh composed of edges in 2d, triangles in 3d, accounting for the boundary $\partial\Omega$ and for internal interfaces, delimiting distinct material regions within $\Omega$. Numerous methods are available when it comes to creating such a mesh $\calT$, depending on how the information related to $\Omega$ is supplied. Often, a surface mesh $\calS$ of the boundary $\partial \Omega$ is provided, which has for instance been constructed thanks to a CAD software; the interior of $\Omega$ is then filled with simplices agreeing with the surface elements of $\calS$. The most efficient procedures to fulfill this goal are the constrained Delaunay algorithm, or advancing front strategies, \cite{borouchaki2017meshing,botsch2010polygon,frey2007mesh}. Without entering into details, let us stress that, in spite of its importance and the attention that has been brought to its analysis for several decades, this task remains delicate. One crucial aspect of meshes is their \emph{quality}, a notion which actually takes on two quite different natures, both illustrated on \cref{fig.meshqual}. \begin{itemize} \item The \emph{finite element quality}. The accuracy of most numerical simulations -- conducted with e.g. the finite element method -- is strongly dependent on how close the elements of the computational mesh $\calT$ are from being regular, see for instance~\cite{ciarlet2002finite} about this classical issue. In practice, a quality factor $\calQ(T)$ is used to evaluate the aspect ratio of each simplex $T\in \calT$, with the meaning that $\calQ(T) \approx 1$ when $T$ is close to regular, and $\calQ(T) \approx 0$ when it is nearly degenerate. One popular criterion used in the literature is: \[ \calQ(T) = \alpha \frac{\Vol(T)}{\left(\sum_{i=1}^{d(d+1)/2} |e_i |^2 \right)^{\frac{d}{2}}}, \] where the $e_i$ are the edges of $T$, and $\alpha$ is a normalization factor. \item The \emph{geometric quality}. Often in practice, the considered domain $\Omega$ (or the internal regions within $\Omega$) is not polygonal; the surface triangulation $\calS_{\calT}$ is only an approximation of $\partial \Omega$, and it is crucial to ensure that this approximation is accurate enough. This may for instance be expressed by demanding that the Hausdorff distance $d^H(\calS_{\calT}, \partial \Omega)$ between $\calS_\calT$ and the ``ideal'', continuous boundary $\partial \Omega$ be smaller than a user-defined tolerance $\e$. \end{itemize} \begin{figure} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/wheelbad} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/wheelok} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/falbad} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/falk} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{center} \caption{(a) Ill-shaped mesh of a wheel; (b) high quality mesh of the same geometry; (c) mesh accounting for a poor geometric approximation of a plane; (d) fine geometric approximation of the same plane.}\label{fig.meshqual} \end{figure} \subsection{``Geometric'' shape optimization}\label{sec.geom} One early attempt to implement the abstract program of \cref{algo.sketchsg} follows an intuitive ``Lagrangian'' strategy, see e.g.~\cite{mohammadi2010applied} or~\cite{pironneau1982optimal}. At each iteration $n$, the shape $\Omega^n$ is explicitly represented by means of a computational mesh $\calT^n$. This choice is particularly convenient for the first stage of \cref{algo.sketchsg}, since the solution $u_{\Omega^n}$ to the linear elasticity system \cref{eq.elas} (and likewise, that $p_{\Omega^n}$ to the adjoint system) can be accurately computed thanks to the finite element method; in principle, any commercial solver can be used in a non intrusive way to this end. The calculation of a descent direction $\theta^n$ for the shape optimization problem \cref{eq.sopb}, which is the second stage of \cref{algo.sketchsg}, is easily conducted from these data. The major difficulty posed by this framework is related to the final stage of \cref{algo.sketchsg}, which is about passing from (the mesh $\calT^n$ of) $\Omega^n$ to (a mesh $\calT^{n+1}$ of) $\Omega^{n+1} := (\Id + \tau^n\theta^n)(\Omega^n)$. Certainly, one mesh $\calT^{n+1}$ of $\Omega^{n+1}$ can be obtained by relocating the vertices of $\calT^n$ according to the rule \begin{equation}\label{eq.relocvertex} \forall \text{vertex } x \text{ of } \calT, \quad x \longmapsto x + \tau^n \theta^n(x), \end{equation} while keeping the connectivities of the mesh unchanged, see \cref{fig.movmesh}. This simple practice unfortunately suffers from serious limitations, since some elements of the mesh are prone to become seriously ill-shaped in the process, not to say downright invalid, see \cref{fig.movmesh}$\MK$(c). \begin{figure}[!ht] \centering \resizebox{\textwidth}{!}{ \begin{tabular}{ccc} \begin{overpic}[width=0.36\textwidth]{figures/cantidisp} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} & \begin{overpic}[width=0.33\textwidth]{figures/cantiok} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} & \begin{overpic}[width=0.36\textwidth]{figures/cantinvalid} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \end{tabular}} \caption{(a) Mesh $\calT^n$ of the shape $\Omega^n$ with the descent direction $\theta^n$ discretized at its vertices; (b) updated mesh $\calT^{n+1} $ of $\Omega^{n+1} :=(\Id + \tau^n \theta^n)(\Omega^n)$ obtained by using the rule \cref{eq.relocvertex} with a ``small enough'' time step $\tau^n$; (c) invalid mesh $\calT^{n+1}$ when the chosen time step $\tau^n$ is ``too large''.}\label{fig.movmesh} \end{figure} Admittedly, such an update of the computational mesh can be conducted in a relatively efficient manner thanks to a number of heuristics. For instance, \begin{itemize} \item One may extend the displacement field $\theta^n$ from the boundary $\partial \Omega^n$ to the interior vertices by solving an elasticity system, with the hope that the extended field induces little compression of the mesh elements. Large mesh displacements have been realized based on this idea in~\cite{baker1999dynamic}. \item The displacement \cref{eq.relocvertex} of the vertices of $\calT^n$ can be realized within several sub-stages intertwined with remeshing operations, whose aim is to improve the quality of elements and thereby to postpone the onset of degenerate or invalid elements, insofar as possible, see for instance~\cite{barral2019three,hassan2007method} in this direction. \end{itemize} The so-called Deformable Simplicial Complex (DSC) method, which leverages such ideas together with further heuristics for coping with topological changes, has recently achieved impressive motions of shapes in the course of the shape and topology optimization process, see~\cite{christiansen2015combined, christiansen2014topology}. \subsection{Level set methods for shape and topology optimization}\label{sec.lsso} The level set method is a general paradigm for tracking dramatic evolutions of domains or interfaces, which may even feature topological changes. Since its inception in~\cite{osher1988fronts} in the context of curvature-driven interface motion, it has proved to be a very efficient and robust framework for fluid simulations and image processing, to name just a few applications. We outline the basic stakes of this method in \cref{sec.lsgen}, referring to~\cite{giga2006surface,osher2006level,sethian1999level} for more exhaustive presentations. The use of this method in the context of shape optimization, as an attempt to overcome the aforementioned meshing issues raised by the update of the shape, is then discussed in the next \cref{sec.lsspec}. \subsubsection{A brief reminder about the level set method}\label{sec.lsgen} Let $D \subset \R^d$ be a fixed computational domain. The level set method is a general philosophy whereby an arbitrary shape $\Omega \subset D$ can be equivalently described as the negative subdomain of a scalar ``level set'' function $\phi : D \to \R$, i.e.: \begin{equation}\label{eq.ls} \begin{cases} \phi(x) < 0 & \text{if } x \in \Omega, \\ \phi(x) = 0 & \text{if } x \in \partial \Omega, \\ \phi(x) > 0 & \text{if } x \in D \setminus \overline \Omega, \end{cases} \end{equation} see \cref{fig.lsrep} for an illustration. \begin{figure}[!ht] \centering \resizebox{\textwidth}{!}{ \begin{tabular}{cc} \begin{minipage}{0.46\textwidth} \begin{overpic}[width=1.\textwidth]{figures/cantidom} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \end{minipage}& \begin{minipage}{0.54\textwidth} \begin{overpic}[width=1.\textwidth]{figures/cantigraphlsb} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{minipage} \end{tabular}} \caption{(a) One shape $\Omega \subset \R^2$; (b) graph of an associated level set function $\phi$ defined on (a mesh of) a larger domain $D$.}\label{fig.lsrep} \end{figure} This representation conveniently allows to reformulate the motion of a shape $\Omega(t)$ over a time period $(0,T)$, with respect to a velocity field $V(t,x)$, in terms of an associated level set function $\phi(t,\cdot\,)$ (i.e. \cref{eq.ls} holds at every time $t >0$). Formally, the latter satisfies the following ``advection-like'' equation: \begin{equation}\label{eq.advls} \begin{dcases} \frac{\partial \phi}{\partial t}(t,x) + V(t,x) \cdot \nabla \phi(t,x) = 0 & \text{for } t \in (0,T), \:\: x\in \R^d, \\ \phi(0,x) = \phi_0(x) & \text{for } x \in \R^d, \end{dcases} \end{equation} where $\phi_0$ is a level set function for the initial shape $\Omega(0)$. Hence, the intricate geometric evolution problem of $\Omega(t)$ translates into the partial differential equation \cref{eq.advls} on the fixed domain $D$. From the practical vantage, the computational domain $D$ is equipped with a fixed mesh $\calT$, for instance a simplicial mesh, or a finite difference grid. The level set function $\phi(t,x)$ and the velocity field $V(t,x)$ are discretized in time and at the vertices of the mesh $\calT$. The evolution equation \cref{eq.advls} can then be solved efficiently thanks to an adapted numerical scheme, see for instance \cref{sec.advect}. \begin{rema}\label{rem.sdf} The above theoretical framework does not rely on any assumption about the nature of the level set function $\phi$ chosen to represent the shape $\Omega \subset D$. Although, in principle, any function satisfying \cref{eq.ls} could be used, it is well-known~\cite{chopp1993computing} that the numerical stability of the level set method is significantly improved when $\phi$ is the signed distance function $d_\Omega$ to $\Omega$. The latter is defined by: \begin{equation}\label{eq.sdf} d_\Omega(x) = \begin{cases} -d(x,\partial\Omega) & \text{if } x \in \Omega, \\ 0 & \text{if } x \in \partial \Omega, \\ d(x,\partial \Omega) & \text{if } x \in D \setminus \overline \Omega, \end{cases} \end{equation} where \begin{equation}\label{eq.disteucl} d(x, \partial \Omega) := \min_{p \in \partial \Omega} |x-p| \end{equation} is the usual Euclidean distance from $x$ to $\partial \Omega$. Among others, $d_\Omega$ enjoys the desirable ``Eikonal'' property, whereby its gradient has unit norm wherever it is defined: \begin{equation}\label{eq.Eikonal} | \nabla d_\Omega (x) | = 1\:\: \text{ for a.e. } x \in D, \end{equation} which expresses a regular spacing of its level sets. Aside from its relevance in the context of the level set method, the signed distance function $d_\Omega$ enjoys multiple interesting properties related to the geometry of $\Omega$. Hence, the calculation of the signed distance function to a shape is a topic of interest on its own, see e.g.~\cite{allaire2016thickness} about the modeling of thickness constraints in shape optimization using distance functions. \end{rema} \subsubsection{Application of the level set method to shape and topology optimization}\label{sec.lsspec} The application of the level set method in the context of shape and topology optimization was originally proposed in~\cite{allaire2004structural, osher2001level, sethian2000structural, wang2003level}. It brings into play a computational domain $D$ equipped with a fixed mesh $\calT$ (e.g. simplicial or Cartesian). At any time $t$, the shape $\Omega(t) \subset D$ is represented by a level set function $\phi(t,\cdot\,) : D \to \R$, discretized, e.g., at the vertices of $\calT$. The update of the shape between any two iterations $n$, $(n+1)$ of the optimization process is efficiently achieved by solving the equation \cref{eq.advls} for the evolution of $\phi(t,x)$. In the present context, the velocity field $V(t,x)$ is the descent direction $\theta^n(x)$ for the shape optimization problem \cref{eq.sopb}, which depends on the elastic displacement $u_{\Omega^n}$ and the adjoint state $p_{\Omega^n}$. The bottleneck of this approach lies in the calculation of $u_{\Omega^n}$ and $p_{\Omega^n}$. Indeed, at a given iteration $n$ of the process (whose reference in notation is omitted for brevity), no mesh of $\Omega$ is available, as $\Omega$ is solely known via the function $\phi$, defined on (the mesh of) the larger domain $D$. One idea to conduct this calculation leverages a so-called ``fictitious domain approach'': the displacement $u_\Omega$ is approximated by the solution $u_\e$ to an equation posed on the total domain $D$. In the present context of linear elasticity, where the traction-free part of shapes is optimized, the void region $D \setminus \overline \Omega$ is filled with a very soft ``ersatz'' material, with Hooke's law $\e A$, $\e \ll 1$, so that an approximate counterpart to \cref{eq.elas} is given by the following equation posed on the fixed domain~$D$: \begin{equation}\label{eq.ersatz} \begin{cases} -\dv (A_\e e(u_\e)) = 0 & \text{in } D, \\ u_\e = 0 & \text{on } \Gamma_D, \\ A_\e e(u_\e) n = g & \text{on } \Gamma_N, \\ A_\e e(u_\e) n = 0 & \text{on } \Gamma, \\ \end{cases}\quad \text{where} \quad A_\e(x) = \begin{cases} A & \text{if } x \in \Omega, \\ \e A & \text{if } x \in D \setminus \Omega. \end{cases} \end{equation} see for instance~\cite{dambrine2010ersatz} for a justification of this procedure. In practice, the tensor $A_\e$ is easily calculated from the knowledge of the level set function $\phi$ for $\Omega$. In a different spirit, advanced finite element techniques, featuring enriched basis functions, have been employed to ensure a more accurate approximation of the displacement $u_\Omega$ in such a context where only a fixed mesh of a large computational domain is available; see~\cite{duysinx2006generalized} about the use of the eXtended Finite Element Method (XFEM) and~\cite{burman2018shape} about that of the cut Finite Element Method (cutFEM). Let us point out that both approaches to cope with the absence of a body-fitted mesh for $\Omega$ are intrusive, insofar as they do not lend themselves to a black-box use of a commercial finite element solver. A typical implementation of the level set method for the shape optimization problem \cref{eq.sopb} is sketched in \cref{algo.lsso}. Note that we deliberately omit the details of Step (2), about the practical calculation of a descent direction from the derivatives of the shape functionals at stake, as it will be the focus of the later \cref{sec.nullspacealgo}. \begin{algorithm}[ht] \caption{The level set method for shape and topology optimization.}\label{algo.lsso} \begin{algorithmic}[0] \STATE \textbf{Initialization:} \noindent \begin{itemize} \item Mesh $\calT$ of the computational domain $D$; \item Level set function $\phi^0: D \to \R$ representing the initial shape $\Omega^0$. \end{itemize} \FOR{$n=0,\dots,$ until convergence} \STATE \begin{enumerate}\relabelitem \item \label{2.1} Calculate an approximate version of the elastic displacement $u_{\Omega^n}$ and of the adjoint state $p_{\Omega^n}$ on $\calT$ by solving the ersatz material problem \cref{eq.ersatz}. \item \label{2.2} From the theoretical formulas for the shape derivatives $J^\prime(\Omega)(\theta)$, $G^\prime(\Omega)(\theta)$ and $H^\prime(\Omega)(\theta)$, infer a descent direction $\theta^n: D \to \R^d$ for \cref{eq.sopb} according to the selected constrained optimization algorithm. \item \label{2.3} Solve the level set advection equation \cref{eq.advls} on $\calT$ with (time-independent) velocity $V(t,x) = \theta^n(x)$, final time $T = \tau^n$ and initial datum $\phi_0 = \phi^n$; a level set function $\phi^{n+1}$ for the new shape $\Omega^{n+1}$ is obtained. \end{enumerate} \ENDFOR \RETURN Level set function $\phi^n$ for the optimized design $\Omega^n$. \end{algorithmic} \end{algorithm} Summarizing, the level set method conveniently alleviates meshing issues and allows to describe dramatic evolutions of the optimized shape $\Omega$, including changes in its topology. Unfortunately, it raises the issue of solving the mechanical equations posed on $\Omega$: although the aforementioned fictitious domain approaches (notably, the ersatz material method) and extended finite element methods work reasonably well in the context of linear elasticity, corresponding strategies are not readily available in more challenging physical contexts, such as that of fluid-structure interactions. The level set based mesh evolution strategy, presented in the next section, is an increment over this level set method for shape and topology optimization: it enjoys all the assets of the latter, and additionally features an exact mesh of the shape at each stage of the process. \section{The level-set based mesh evolution method for shape and topology optimization}\label{sec.bf} We now turn to the numerical framework introduced in~\cite{allaire2011topology,allaire2013mesh,allaire2014shape} as an attempt to overcome the individual difficulties posed by Lagrangian and level set methods when tracking the evolution of the shape during the optimization process. To set ideas, the presentation unfolds in the structural optimization context of \cref{sec.sopb}, where one aims to solve a problem of the form \cref{eq.sopb}. After sketching the main stages of the method in \cref{sec.bfpres}, we overview its main ingredients in the subsequent sections. \subsection{General description of the method}\label{sec.bfpres} As the name suggests, the level set based mesh evolution method at stake in this article combines the meshed and level set representations of shapes discussed in \cref{sec.geom,sec.lsso}. Efficient algorithms make it possible to switch from one to the other, so that the most convenient of them with respect to the ongoing operation can be used. Let $D$ be a large ``hold-all'' domain containing all the considered shapes $\Omega$. At each iteration $n$ of the process, $D$ is endowed with a valid and conforming mesh $\calT^n$, which is modified from one iteration to the other so that the current shape $\Omega^n$ is explicitly discretized. More precisely, $\calT^n$ is consistently made of two complementary submeshes $\caliT^n$, $\caleT^n$ such that $\caliT^n$ is a valid, conforming mesh of $\Omega^n$ and $\caleT^n$ is a valid, conforming mesh of the exterior part $D \setminus \overline{\Omega^n}$. Thus, two descriptions of $\Omega^n$ are available: \begin{itemize} \item \emph{A meshed representation.} $\Omega^n$ is explicitly represented by the submesh $\caliT^n$ of $\calT^n$, see \cref{fig.bfmethod}$\MK$(a). \item \emph{A level set representation.} $\Omega^n$ is implicitly defined by a level set function $\phi^n$, defined on the whole mesh $\calT^n$ of $D$, see \cref{fig.bfmethod}$\MK$(b). \end{itemize} On the one hand, the meshed representation is handful when it comes to solving partial differential equations on $\Omega^n$, such as that \cref{eq.elas} for the elastic displacement $u_{\Omega^n}$; on the other hand, the level set description is adequate for conducting the update of the shape $\Omega^n$ into the next iterate $\Omega^{n+1}$ in a robust way. As we have mentioned, this strategy crucially hinges on efficient numerical methods for passing from one of these descriptions to the other, and notably: \begin{itemize} \item An algorithm for generating one level set function $\phi : D \to \R$ for a shape $\Omega \subset D$ which is explicitly discretized inside a mesh $\calT$ of the computational domain $D$; \item An algorithm which assumes the datum of a level set function $\phi : D \to \R$ defined on a mesh $\calT$ of $D$, associated to a shape $\Omega \subset D$, and creates a new mesh $\widetilde{\calT}$ of $D$ in which $\Omega$ is explicitly discretized. \end{itemize} This method is outlined in \cref{algo.bf} in the context of the model shape optimization problem of \cref{sec.th}. Its main steps are illustrated on \cref{fig.bfmethod}; as they constitute operations of general interest, we present them in a self-contained way in the next sections. The computation of the signed distance function to a shape is discussed in \cref{sec.sdf}, and the resolution of the level set advection equation is presented in \cref{sec.advect}. The aspects of the method concerned with remeshing are addressed in \cref{sec.remesh}; we then broach the versatile Hilbertian extension and regularization procedure in \cref{sec.descentth}, before finally describing the numerical constrained optimization algorithm in \cref{sec.nullspacealgo}. \begin{algorithm}[ht] \caption{The level set based mesh evolution method.}\label{algo.bf} \begin{algorithmic}[0] \STATE \textbf{Initialization:} Mesh $\calT^0$ of the computational domain $D$, enclosing a mesh $\caliT^0$ of $\Omega^0$ as a submesh. \FOR{$n=0,\dots,$ until convergence} \STATE \begin{enumerate}\relabelitem \item \label{3.1} Calculate the signed distance function $\phi^n = d_{\Omega^n}$ to $\Omega^n$ at the vertices of $\calT^n$. \item \label{3.2} Calculate the displacement $u_{\Omega^n}$ and the adjoint state $p_{\Omega^n}$ by solving the partial differential equation \cref{eq.elas} on the interior part $\caliT^n$ of $\calT^n$. \item \label{3.3} Use the Hilbertian procedure to calculate gradients $\theta_J^n$, $\theta_G^n$ and $\theta_H^n$ for $J(\Omega)$, $G(\Omega)$ and $H(\Omega)$ at $\Omega = \Omega^n$ on the whole mesh $\calT^n$. \item \label{3.4} Infer a descent direction $\theta^n$ on $\calT^n$ for \cref{eq.sopb} thanks to a constrained optimization algorithm. \item \label{3.5} Choose a small enough time step $\tau^n$ and calculate a level set function $\widetilde{\phi}^{n+1}$ for the new shape $\Omega^{n+1} := (\Id + \tau^n \theta^n)(\Omega^n)$ on the mesh $\calT^n$ by solving the level set evolution equation \cref{eq.advls} over $(0,\tau^n)$ with velocity field $V(t,x) = \theta^n(x)$ and initial condition $\phi_0 =\phi^n$. \item \label{3.6} Create a new mesh $\calT^{n+1}$ of $D$ where $\Omega^{n+1}$ is explicitly discretized, from the datum of the level set function $\widetilde\phi^{n+1}$ on $\calT^n$. \end{enumerate} \ENDFOR \RETURN Mesh $\calT^n$ of $D$ where $\Omega^n$ is discretized as a submesh $\caliT^n$. \end{algorithmic} \end{algorithm} \begin{figure} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/ini} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.47\textwidth]{figures/phib} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/solelasb} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/gradb} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/tempnew} \put(0,3){\fcolorbox{black}{white}{$e$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/temprem} \put(0,3){\fcolorbox{black}{white}{$f$}} \end{overpic} \end{center} \caption{Illustration of the main stages of the mesh evolution method sketched in \cref{sec.bfpres}. (a) Mesh $\calT^n$ of the computational domain $D$; the submesh $\caliT^n$ of $\Omega^n$ consists of the black elements, and that $\caleT^n$ of $D \setminus \overline{\Omega^n}$ is made of the white elements; (b) isovalues of the signed distance function $\phi^n$ to $\Omega^n$, calculated on $\calT^n$; (c) solution $u_{\Omega^n}$ to the elasticity system on the mesh $\caliT^n$ of $\Omega^n$; (d) descent direction $\theta^n$ on the whole mesh $\calT^n$; (e) level set function $\widetilde \phi^{n+1}$ on the mesh $\calT^n$; the $0$ level set of $\widetilde \phi^{n+1}$ is depicted in red; (f) new mesh $\calT^{n+1}$ with the new shape $\Omega^{n+1}$ enclosed as the submesh $\caliT^{n+1}$ (made of the black elements).}\label{fig.bfmethod} \end{figure} \begin{rema} Appealing features of this approach are its modularity and non intrusiveness: the aforementioned steps can be carried out with independent numerical codes. In particular, the resolution of the state (and adjoint) equation for the elastic displacement does not demand any exotic treatment; it can be realized with any solver, without entering into its implementation details. \end{rema} \subsection{Calculation of the signed distance function}\label{sec.sdf} This section deals with the construction of one particular level set function $\phi$ for a shape $\Omega \subset D$, namely the signed distance function $d_\Omega$ in \cref{eq.sdf}, out of a meshed representation of the latter (Step~\eqref{3.1} in \cref{algo.bf}), see the discussion in \cref{rem.sdf}. Let the computational domain $D$ be endowed with a simplicial mesh $\calT$ and let $\Omega \subset D$ be a shape; we wish to calculate the values $\phi(x) = d_\Omega(x)$ at all vertices $x\in\calT$. Note that in the shape optimization workflow of \cref{sec.bfpres}, $\Omega$ is supplied as an explicit submesh $\caliT$ of $\calT$, but for the purpose of this section, it could actually be given under a different format, for instance, via a proper mesh $\calT_\Omega$ which is not necessarily a part of $\calT$. The numerical calculation of $d_\Omega$ can be conducted in a variety of manners. ``Geometric'' algorithms involve an exhaustive calculation of the distance $d(x,\partial \Omega)$ from $x$ to $\partial \Omega$ at every vertex $x$ of $\calT$, defined as the minimum value featured in \cref{eq.disteucl}. Although this operation can be made efficient owing to a number of heuristics -- see e.g.~\cite{mauch2000fast,tsai2002rapid} and the references therein -- we rather rely on ``propagation algorithms'', which comprise two steps: \begin{enumerate} \item The function $\phi$ is initialized with (an approximation of) the exact value of $d_\Omega$ at the vertices of $\calT$ which are ``close'' to $\partial \Omega$ (for instance, at the vertices of the simplices $T \in \calT$ intersecting $\partial \Omega$), and with large, positive or negative values elsewhere. This stage is elementary; however, depending on the input format and of the complexity of the geometry of $\Omega$, it may prove tedious from the implementation viewpoint, and time-consuming in practice, see e.g.~\cite{dapogny2012computation}. \item The calculation of the signed distance $\phi(x) = d_\Omega(x)$ is realized from the vertices $x \in \calT$ closest to $\partial \Omega$ to farther ones, by relying on a discretization of the Eikonal equation \cref{eq.Eikonal}.\pagebreak{} This purpose is greatly simplified when the computational mesh $\calT$ is a Cartesian grid, since high-order finite different schemes are available. \end{enumerate} The most popular algorithm in this second category is certainly the fast marching method, see~\cite{sethian1999fast} for an overview and~\cite{kimmel1998computing} for an adaptation to the case where the computational mesh is simplicial; let us also mention the fast sweeping method~\cite{zhao2005fast}. In the present work, we rely on an algorithm based on the properties of the so-called redistancing equation, which is described in~\cite{dapogny2012computation} and is available in the open-source library \texttt{mshdist}\footnote{\texttt{\url{https://github.com/ISCDtoolbox/Mshdist}}}. \subsection{Advection of the level set function}\label{sec.advect} In this section, we turn to the numerical resolution of the equation \cref{eq.advls} accounting for the evolution of the shape. The generic situation, occurring at each iteration $n$ of the process described in \cref{sec.bfpres}, is the following: a level set function $\phi = \phi^n : D \to \R$ accounting for the current shape $\Omega = \Omega^n$ is supplied via its values at the vertices of a simplicial mesh $\calT = \calT^n$ of the computational domain $D$. Analogously, the velocity field $V : D \to \R^d$ driving the evolution of the shape, which is the descent direction $\theta^n$ at the current iteration, is supplied at the vertices of $\calT$. We aim to compute the solution $\psi(t,x)$ to the following advection equation: \begin{equation}\label{eq.advection} \begin{dcases} \frac{\partial \psi}{\partial t}(t,x) + V(x) \cdot \nabla \psi(t,x) = 0 & \text{for } t \in (0,T), \:\: x \in D, \\ \psi(0,x) = \phi(x) & \text{for } x \in D, \end{dcases} \end{equation} and notably its values $\psi(T,x)$ at the final time $t = T$, which stands for the time step $\tau^n$. Numerical methods for the resolution of \cref{eq.advection} are manifold when the computational support is a Cartesian grid, which allows for the use of high-order finite difference methods. This issue is however a little less classical and deserves a few comments in our context where it is a simplicial mesh. We rely on the method of characteristics proposed in~\cite{pironneau1989finite}, which is close in spirit to the semi-Lagrangian scheme developed in~\cite{strain1999semi}; see alternatively~\cite{ern2006discontinuous} about the use of discontinuous Galerkin methods, and~\cite{abgrall1996numerical,barth1998numerical} for the construction of numerical schemes for more general Hamilton--Jacobi equations on simplicial meshes. The method of characteristics is based on the analytical formula for the solution to \cref{eq.advection}. The latter is expressed in terms of the characteristic curves $t \mapsto X(t,t_0,x)$ of the velocity field $V(x)$, emerging from an arbitrary point $x \in D$ at a time $t_0$, defined as the solution to the ordinary differential equation: \begin{equation}\label{eq.characcurv} \begin{dcases} \frac{\dd}{\dt} X(t,t_0,x) = V(X(t,t_0,x)) & \text{for } t \in \R, \\ X(t_0,t_0,x) = x.& \end{dcases} \end{equation} Intuitively, $t\mapsto X(t,t_0,x)$ is the trajectory of a particle located in $x$ at time $t=t_0$, which is transported according to the velocity field $V(x)$. The exact solution to \cref{eq.advection} is then given by: \begin{equation}\label{eq.solcharac} \forall t \in (0,T), \:\: x \in D, \quad \psi(t,x) = \phi(X(0,t,x)), \end{equation} which expresses the natural fact that the value of $\psi$ at time $t$ and point $x$ is the value of the initial datum $\phi$ at the initial position $X(0,t,x)$ of the particle lying in $x$ at time $t$. One simple means to exploit the closed-form formula \cref{eq.solcharac} is to directly discretize it: for each vertex $x \in \calT$, the ordinary differential equation \cref{eq.characcurv} is solved for the ``backward'' characteristic curve $t \mapsto X(t,T,x)$, for instance by a Runge--Kutta 4 scheme; $\phi$ is then evaluated at the ``foot'' $X(0,T,x)$ of this characteristic line. \begin{rema} \CD{In practice, the origin $X(0,T,x)$ of the characteristic curve passing through $x$ at $t=T$ may lie outside $D$. This notably happens when the velocity field $V(x)$ is pointing inward $D$ on at least one portion of the boundary $\partial D$. In such a case, the information supplied in \cref{eq.advection} is not sufficient to guarantee the well-posedness of this equation (as one would have to add a boundary condition on the ``entrant'' part of $\partial D$); in numerical practice, it is customary to endow $\psi(T,x)$ with a consistent value by extrapolating the value of $\phi$ (and the characteristic curve $t\mapsto X(t,T,x)$) outside $D$.} \end{rema} An open-source implementation of this algorithm is proposed in the \texttt{advection}\footnote{\texttt{\url{https://github.com/ISCDtoolbox/Advection}}} code. \subsection{Using remeshing to pass from a level set to a meshed representation of a shape}\label{sec.remesh} This section is devoted to the delicate operation of passing from a level set to a meshed representation of a shape $\Omega$, which is Step~\eqref{3.6} in \cref{algo.bf}. Let us consider the following generic situation: the computational domain $D$ is equipped with a simplicial mesh $\calT$, and a level set function $\phi : D \to \R$ for a shape $\Omega \subset D$ is provided via its values at the vertices of $\calT$; we aim to construct a new mesh $\widetilde \calT$ of $D$ in which $\Omega$ appears as a submesh, see \cref{sec.bfpres}. This task is achieved in two steps: \begin{enumerate} \item The $0$ level set of $\phi$ is explicitly discretized into $\calT$. This crude procedure yields a valid, conforming mesh $\calTt$ of $D$, in which $\Omega$ exists as a submesh; unfortunately, $\calTt$ has poor finite element quality, and the interior submesh is a poor geometric approximation of the shape $\Omega$; see \cref{fig.mmgls}$\MK$(b). \item Local remeshing operations are applied to modify $\calTt$ into a new, high-quality mesh $\widetilde{\calT}$ of $D$ in which $\Omega$ is still explicitly discretized, see \cref{fig.mmgls}$\MK$(c). \end{enumerate} These two steps are presented with a little more details in the next \cref{sec.isodisc,sec.remeshgen}, respectively. They are implemented in the general purpose open-source library \texttt{mmg}\footnote{\texttt{\url{https://github.com/MmgTools/mmg}}} devoted to simplicial remeshing; we refer to~\cite{dapogny2014three} for a more exhaustive presentation of the latter, and to~\cite{balarac2021tetrahedral} for recent developments on the subject. Throughout the rest of this section, for notational simplicity, the mesh of $D$ will be systematically denoted by $\calT$, although it is constantly subject to modifications. \begin{figure} \begin{center} \begin{overpic}[width=0.33\textwidth]{figures/mmgiso} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.311\textwidth]{figures/mmgbad} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \begin{overpic}[width=0.315\textwidth]{figures/mmgok} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \end{center} \caption{(a) Isovalues of a level set function discretized at the vertices of a mesh $\calT$ of a computational domain $D$; (b) Ill-shaped mesh $\calTt$ resulting from the rough discretization of the $0$ level set of $\phi$ into $\calT$; (c) High-quality mesh $\widetilde{\calT}$ obtained from $\calTt$ by remeshing.}\label{fig.mmgls} \end{figure} \subsubsection{Explicit discretization of the \texorpdfstring{$0$}{0} level set of \texorpdfstring{$\phi$}{phi} in the mesh \texorpdfstring{$\calT$}{T}}\label{sec.isodisc} From a simplicial mesh $\calT$ of $D$ and a level set function $\phi :D \to \R$ for a subdomain $\Omega \subset D$, supplied at the vertices of ${\mathcal T}$, we aim to construct a new, valid but possibly ill-shaped mesh of $D$ in which $\Omega$ is explicitly discretized. To achieve this, we rely on a simplicial variant of the well-known marching cubes algorithm~\cite{lorensen1987marching}, namely the marching triangles (in 2d) or tetrahedra (in 3d) algorithm~\cite{doi1991efficient}. Briefly, we first identify all the elements $T \in \calT$ which intersect the $0$ level set of $\phi$, by looking at the signs of $\phi$ at their vertices. After linear interpolation of $\phi$ inside any of these elements $T$, the intersection $\partial \Omega \cap T$ is a portion of line (in 2d) or plane (in 3d) which is completely determined by the values of $\phi$ at the vertices of $T$. Then, using a pre-defined pattern, $T$ is split into several sub-triangles (in 2d) or sub-tetrahedra (in 3d), so that $\partial \Omega \cap T$ explicitly appears in the resulting mesh, see \cref{fig.pattern}. \begin{figure} \begin{center} \begin{overpic}[width=0.75\textwidth]{figures/cuttria} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.75\textwidth]{figures/cuttetra} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \caption{(a) One pattern for splitting a triangle, based on the values of $\phi$ at its vertices ($\phi$ is positive at the red vertex, negative at the blue ones); (b) one possible pattern for splitting a tetrahedron.}\label{fig.pattern} \end{figure} The mesh $\calT$ resulting from this rough operation is valid, conforming, and it encloses two submeshes $\caliT$ and $\caleT$ of $\Omega$ and $D \setminus \overline \Omega$, respectively. Unfortunately, it generally contains elements with very bad quality, and it accounts for a poor geometric approximation of $\partial \Omega$. \subsubsection{Improvement of the quality of \texorpdfstring{$\calT$}{T} by local remeshing operations}\label{sec.remeshgen} The algorithm of the previous section has produced a valid and conforming mesh $\calT$ of the domain $D$, in which the considered shape $\Omega \subset D$ explicitly appears as a submesh. Unfortunately, $\calT$ has bad quality, in both senses introduced in \cref{sec.mesh}: it contains nearly degenerate elements, and it accounts for a poor geometric approximation of the boundary of $D$ and of its internal interfaces, in particular of the boundary $\partial \Omega$ of the shape, see again \cref{fig.mmgls}$\MK$(b). We wish to modify $\calT$ by repeatedly applying local remeshing operations, so as to improve both features.\pagebreak{} The presentation of this section focuses on the case of three space dimensions, since it contains numerous specificities, and it is on any aspect more involved than its 2d counterpart. One immediate issue posed by the desire to better approximate the domains $D$ and $\Omega$ is that no continuous description is provided: the only available information about $D$ or $\Omega$ is the supplied discrete mesh $\calT$ itself. To cope with this difficulty, we first reconstruct local geometric information about $D$ and $\Omega$ from the discrete datum of $\calT$; this information will be updated throughout the remeshing process. More precisely, normal vectors $n(x)$ are calculated at the vertices $x$ of the surface mesh $\calS_\calT$, for instance as weighted sums of the normal vectors to the nearby surface triangles. Besides, remarkable geometric entities such as sharp edges and corner points are identified. This information allows to infer a piece of the ``ideal'' continuous surface $\partial D$ or $\partial \Omega$ associated to a given surface triangle $T \in \calS_\calT$: a cubic B\'ezier parametrization $\sigma : \widehat T \to \R^3$ (defined on the reference simplex $\widehat T \subset \R^2$) is calculated, which passes through the three vertices $a_0$, $a_1$, $a_2$ of $T$ and fits the attached normal vectors $n_0$, $n_1$, $n_2$, see \cref{fig.bezier}$\MK$(a). This local continuous model is helpful for a variety of purposes, in particular when it comes to measuring ``how far'' $\calT$ stands from the continuous domains $D$ and $\Omega$, as depicted in \cref{fig.bezier}$\MK$(b). \begin{figure} \resizebox{\textwidth}{!}{ \quad \begin{tabular}{cc} %\centering \begin{minipage}{0.45\textwidth} \begin{overpic}[width=1.\textwidth]{figures/bezier} \put(-15,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \end{minipage} & \begin{minipage}{0.35\textwidth} \begin{overpic}[width=1.02\textwidth]{figures/beziercontrol} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{minipage} \end{tabular} } \caption{(a) Creation of a cubic B\'ezier patch for the region of $\partial \Omega$ associated to a surface triangle $T \in {\mathcal S}_{\mathcal T}$; (b) measurement of the gap between the continuous and discrete domains from this local patch.}\label{fig.bezier} \end{figure} Another key component of the remeshing strategy is the construction and the use of a size map $h : D \to \R$, which encodes at each vertex $x \in \calT$ the desired size for edges surrounding $x$. This value is calculated from user-defined requirements, such as lower and upper bounds for the size of edges in $\calT$, or a more detailed local size prescription stemming from a priori or a posteriori error estimates related to the finite element method. Additionally, when $x$ is a surface vertex, the calculation of $h(x)$ has to take into account the needed local size to comply with the geometric approximation requirements. The remeshing procedure then starts, properly speaking. The elements of $\calT$ are travelled repeatedly; the action of four local remeshing operators is simulated, and their outcome is effectively retained if it shows an improved mesh quality. These operators are quite commonly used in remeshing practice; their action is exemplified on \cref{fig.meshops} in the case of two space dimensions, and it can be summarized as follows: \begin{itemize} \item \emph{Edge split.} When an edge $pq$ in the mesh is ``too long'' (with respect to the size map~$h$), it is split into two edges $pm$ and $qm$ after introduction of a new point $m$ in $\calT$; the tetrahedra that were formerly in the shell of $pq$ --~i.e. that were sharing $pq$ as an edge~-- are appropriately reconnected. \item \emph{Edge collapse.} When an edge $pq$ is ``too short'', $q$ is merged with $p$; the tetrahedra in the shell of $pq$ are suppressed, and the other elements formerly connected to $q$ are suitably updated. \item \emph{Edge swap.} One edge $pq$ in the mesh is suppressed and the elements in the shell of $pq$ are reconnected so that the mesh remains valid and conforming. \item \emph{Vertex relocation.} One vertex $p$ of the mesh is assigned to a new position, while all the mesh connectivities are unaltered. \end{itemize} \begin{figure} \centering \begin{overpic}[width=0.8\textwidth]{figures/split} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.8\textwidth]{figures/collapse} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \begin{overpic}[width=0.8\textwidth]{figures/swap} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic}\\[0.25em] \begin{overpic}[width=0.8\textwidth]{figures/relocate} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \caption{Two-dimensional illustration of the four remeshing operations described in \cref{sec.remeshgen}. (a) Split of the edge $pq$: the midpoint $m$ is introduced, and the two red triangles are replaced by the four blue triangles; (b) Collapse of vertex $q$ onto $p$; the ``shell'' of $pq$, made of the two red triangles, is removed; (c) Swap of the edge $pq$; the two red triangles are replaced by the two blue triangles; (d) Relocation of the vertex $p$, while maintaining the connectivities of the mesh.}\label{fig.meshops} \end{figure} \smallskip Each operator exists under two versions, dedicated to surface and internal configurations, respectively. For instance, when a boundary edge $pq$ is split, the new vertex $m$ ought to be inserted directly on the curved segment on $\partial D$ or $\partial \Omega$ which is associated to $pq$ via the local surface model $\sigma$ discussed above; on the contrary, when $pq$ is an internal edge of $\calT$, $m$ is simply introduced at the center of $pq$. \subsection{Extension - regularization of shape derivatives via the Hilbertian method}\label{sec.descentth} The shape derivatives $J^\prime(\Omega)(\theta)$, $G^\prime(\Omega)(\theta)$ and $H^\prime(\Omega)(\theta)$ of the objective and constraint functions $J(\Omega)$, $G(\Omega)$ and $H(\Omega)$ are continuous linear forms on $\Winfty$, which is a Banach space. The nature and structure of these derivatives do not lend themselves to a simple numerical treatment. This fact is conveniently illustrated by the unconstrained shape optimization problem \cref{eq.sopbuc}. As we have mentioned, when the shape derivative of $J(\Omega)$ is known under the structure \cref{eq.surfsd}, a descent direction is immediately revealed as $\theta = -v_\Omega n$, see \cref{eq.thetavn}. However appealing, this choice is awkward, since this formula only characterizes the values of $\theta$ on the boundary $\partial \Omega$, while we have seen that for a variey of purposes, including the practice of Lagrangian mesh update strategies as in \cref{sec.geom}, or that of the level set method described in \cref{sec.lsso}, $\theta$ ought to be ``adequately'' defined on the total computational domain $D$. Moreover, in practice, the choice \cref{eq.thetavn} results in a descent direction which lacks smoothness; this often causes the shape to develop numerical artifacts in the course of the optimization, see~\cite{mohammadi2010applied} for a discussion about this feature. Returning to the treatment of a general problem of the form \cref{eq.sopb}, many efficient infinite-dimensional optimization algorithms leverage a Hilbertian structure when it comes to calculating descent directions, see for instance the algorithm described in the next \cref{sec.nullspacealgo}. One handful practice to circumvent these difficulties, with countless additional benefits, is the so-called Hilbertian extension and regularization method, which consists in endowing the space of deformations $\theta$ with the structure of a Hilbert space, composed of vector fields obeying certain user-specified properties, such as their smoothness, their domain of definition, or imposed values on particular regions of space (e.g. $\theta=0$ where shapes are not subject to optimization). We refer to~\cite{azegami1996domain,burger2003framework,de2006velocity} about this idea, see also~\cite{allaire2020survey} for a summary. The key idea is to introduce a Hilbert space $V$, with inner product $a(\,\cdot\,,\cdot\,)$, which is continuously embedded in $\Winfty$. Thus, the derivative $J^\prime(\Omega)(\theta)$ of a function of the domain $J(\Omega)$ induces a continuous linear form on $V$; according to the Riesz representation theorem, we may represent the latter by the element $\theta_J \in V$ defined by: \begin{equation}\label{eq.Rieszgradient} \forall w \in V, \quad a(\theta_J,w) = J^\prime(\Omega)(w). \end{equation} This element $\theta_J \in V$ is the gradient associated to the derivative $\theta \mapsto J^\prime(\Omega)(\theta)$ via the inner product $a(\,\cdot\,,\cdot\,)$. As expected, $-\theta_J$ is a descent direction for $J(\Omega)$ since \[ J^\prime(\Omega)(-\theta_J) = -a(\theta_J,\theta_J) \leq 0, \] where the latter quantity vanishes if and only if $\Omega$ is already a critical shape for $J(\Omega)$. Multiple possibilities are available as regards the choice of the space $V$ and its inner product $a(\,\cdot\,,\cdot\,)$, see again~\cite{burger2003framework}. For instance, one natural choice is: \[ V = H^m(D)^d, \quad\text{with the usual inner product } a(u,v) := \sum_{|\alpha| \leq m} \int_D \partial^\alpha u \cdot \partial^\alpha v \:\dx, \] where the index $m$ is chosen large enough so that the continuous inclusion $V \subset \Winfty$ holds by virtue of the Sobolev embedding theorem~\cite{adams2003sobolev}. The gradient $\theta_J \in V$ produced by the identification problem \cref{eq.Rieszgradient} is thus a vector field on the whole computational domain $D$ (and not only on the boundary $\partial \Omega$ of the shape), and it enjoys the higher regularity of a vector field in $H^m(D)$. Another popular practice, albeit formal, is that retained in this article: only the normal component of the velocity field $\theta$ is extended and regularized. More precisely, we consider the Hilbert space \begin{equation}\label{eq.velext1} V = H^1(D), \quad\text{with inner product } a(u,v) = \alpha^2 \int_D \nabla u \cdot \nabla v \:\dx+ \int_D uv \:\dx, \end{equation} where $\alpha >0$ is a user-defined parameter, and we solve the following identification problem: \begin{equation}\label{eq.velext2} \text{Search for } v \in V \text{ s.t. } \forall w \in V, \quad a(v,w) = J^\prime(\Omega)(wn), \end{equation} where $n$ stands for an extension of the unit normal vector to $\partial \Omega$. We then take $\theta_J = vn$ as the ``gradient'' of $J(\Omega)$, which is defined on the total computational domain $D$. Although not rigorous, since $H^1(D)^d$ is not embedded into $\Winfty$, this practice yields good results; intuitively, the identification problem \cref{eq.velext1,eq.velext2} amounts to smoothing of the ``$L^2(\partial \Omega)$ gradient'' $v_\Omega$ of $J(\Omega)$ featured in \cref{eq.surfsd} over a thickness $\alpha$ around $\partial \Omega$. \subsection{A null space algorithm for constrained optimization}\label{sec.nullspacealgo} Only relatively few of the numerous constrained optimization algorithms (see e.g.~\cite{nocedal2006numerical} for an overview) lend themselves to an efficient treatment of infinite-dimensional problems, and in particular shape optimization problems. Let us nevertheless mention the article~\cite{svanberg1987method} introducing the popular Method of Moving Asymptotes (MMA) in the context of density-based topology optimization, or that~\cite{dunning2015introducing} about the adaptation of the Sequential Linear Programming (SLP) method to the shape optimization context. In this section, we present a numerical algorithm for dealing with constrained optimization which is particularly well-suited for our applications; it is presented in detail in~\cite{feppon2020null}\footnote{An open source implementation is available at the address: \[ \texttt{\url{https://gitlab.com/florian.feppon/null-space-optimizer}} \]}, see also the article~\cite{barbarosie2018algorithm} where a similar idea was developed independently. For simplicity, we restrict ourselves to a heuristic presentation in the case of an optimization problem featuring only equality constraints, although the method allows to treat an arbitrary number of equality and inequality constraints. We consider a shape optimization problem of the form \begin{equation}\label{eq.optpbdisc} \min_\Omega J(\Omega) \: \text{ s.t. } G(\Omega) = 0, \end{equation} where $G(\Omega) = (G_i(\Omega))_{i=1,\ldots,p} \in \R^p$ is a collection of $p$ real-valued equality constraint functionals. Let us place ourselves at a given iteration $n$ of the execution of \cref{algo.bf}, dropping all superscripts referring to the latter for notational simplicity. The current shape is denoted by $\Omega$, and we assume that a Hilbert space $V$ and an inner product $a(\,\cdot\,,\cdot\,)$ have been chosen for the considered deformation fields $\theta$ according to the Hilbertian method of \cref{sec.descentth}. As discussed in there, the shape derivatives $J^\prime(\Omega)$ and $G_i^\prime(\Omega)$ are accounted for by gradients $\theta_J$, $\theta_{G_i} \in V$, $i=1,\ldots,p$, that is: \[ \forall \xi \in V, \quad J^\prime(\Omega)(\xi) = a(\theta_J,\xi),\text{ and } G_i^\prime(\Omega)(\xi) = a(\theta_{G_i}, \xi). \] The next iterate $(\Id + \tau \theta)(\Omega)$ in the solution of the shape optimization problem \cref{eq.optpbdisc} is obtained from a deformation $\theta$ which is sought as a linear combination of $\theta_J$ and the $\theta_{G_i}$; we may write the latter under the form \begin{equation}\label{eq.thetanoptim} \theta = - \left(\alpha_J \xi_J + \alpha_G \xi_G \right),\quad \text{ where } \begin{cases} \xi_J = \theta_J - \sum_{i=1}^p \lambda_ i \theta_{G_i},\\ \xi_G = \sum_{i=1}^p \beta_i \theta_{G_i}, \end{cases} \end{equation} and the coefficients $\lambda_i, \beta_i$ are characterized by the following requirements: \begin{itemize} \item The vector field $\xi_J \in V$ is a descent direction for $J(\Omega)$ lying in the null-space of the constraint functional $G(\Omega)$, i.e. \[ \forall i = 1,\ldots,p, \quad G^\prime_i(\Omega)(\xi_J) = a(\theta_{G_i}, \xi_J) = 0. \] This property rewrites under the form of a $p \times p$ matrix system for the coefficients $\lambda = (\lambda_i)_{i=1,\ldots,p}$: \[ S \lambda = b, \quad \text{where } S_{ij} := a(\theta_{G_i},\theta_{G_j}), \text{ and } b_i := a(\theta_{G_i}, \theta_J). \] \item The contribution $\xi_G \in V$, which is then the only one able to modify the value of the constraint functional, should ensure a reduction in the absolute value of $G(\Omega)$ by a unit: \[ G((\Id + \tau \theta)(\Omega)) = (1- \alpha_G \tau) G(\Omega). \] This requires that $|1-\alpha_G \tau |< 1$, that is $\alpha_G < \frac{2}{\tau}$. Using the asymptotic expansion \cref{eq.variationtheta}, one obtains, at first order: \[ \forall i =1,\ldots,p, \quad \sum_{j=1}^p \beta_j G_i^\prime(\Omega)(\theta_{G_j}) = - G_i(\Omega), \] which leads to the $p\times p$ matrix system: \[ S \beta = - c, \text{ where } c_i := -G_i(\Omega). \] \end{itemize} Intuitively, the weights $\alpha_J$ and $\alpha_G$ in \cref{eq.thetanoptim} encode the relative rates at which the algorithm imposes the reduction in the value of the objective function and the fulfillment of the constraints. The practical implementation of this strategy relies on the following adaptations of these considerations: \begin{itemize} \item The coefficients $\alpha_J$ and $\alpha_G$ in \cref{eq.thetanoptim} are actually adapted from one iteration to the other, in order to control the maximum amplitude $\| \theta \|_{L^\infty(D)^d}$ of the descent direction $\theta$; we set: \begin{align*} \alpha_J& = \begin{dcases} \frac{A_J h}{\| \xi_J \|_{L^\infty(D)^d}} & \text{if the iteration number }n \text{ is } \leq n_0,\\ \frac{A_J h}{\max(\| \xi_J \|_{L^\infty(D)^d},\Xi)} & \text{otherwise},\\ \end{dcases} \\ \text{and}\quad \alpha_G &= \min\left(\frac{A_G h}{\| \xi_G \|_{L^\infty(D)^d}},2 \right), \end{align*}%\ici where $h$ is the average size of an edge in the mesh, and $A_J$, $A_G$ are fixed ratios accounting for the largest displacement (relatively to the mesh size) of the boundary $\partial \Omega$ induced by the requirements to reduce the value of the objective function and to satisfy the constraints, respectively. In the above formula, a given number of iterations $n_0$ is introduced, and the normalization of $\xi_J$ expressed in the definition of $\alpha_J$ is only applied during the first $n_0$ iterations of the optimization process; when $n > n_0$, the normalization factor is replaced by the maximum between the norm $ \| \xi_J \|_{L^\infty(D)^d}$ of the current direction $\xi_J$, and its value $\Xi$ at iteration $n_0$. This leaves the room for the contribution $\alpha_J \xi_J$, which is in charge of making the value of $J(\Omega)$ decrease, to tend to $0$ as $n$ grows, and thus to avoid oscillations of the algorithm in the final iterations of the method. \item A merit function is used to appraise the size of the step $\Delta t$ during which the direction $\theta$ in \cref{eq.thetanoptim} allows to improve $\Omega$ with respect to the problem \cref{eq.optpbdisc}. More precisely, we introduce the augmented Lagrangian-like shape functional \begin{equation}\label{eq.merit} M(\Omega) := \alpha_J \left(J(\Omega) - \sum_{i=1}^p \lambda_i G_i(\Omega) \right) + \frac{\alpha_G}{2} S^{-1}G(\Omega) \cdot G(\Omega), \end{equation} whose shape gradient (with respect to the inner product $a(\,\cdot\,,\cdot\,)$) is precisely the deformation field $\theta$ in \cref{eq.thetanoptim}. The time step $\Delta t$ must then be chosen in such a way that $M((\Id + \Delta t \theta)(\Omega)) < M(\Omega)$, at least up to a certain tolerance. \end{itemize} \begin{rema} One interesting feature of this algorithm is that it does not involve many user-defined parameters which may prove difficult to tune: only $A_J$ and $A_G$ have to be specified, whose physical meaning is quite clear. The ratio $A_J / A_G$ allows to control the pace at which constraints become satisfied. \end{rema} \section{Numerical illustrations}\label{sec.num} In this section, we present several numerical examples illustrating the main features of our shape optimization framework, in physical situations where a body-fitted approach is particularly desirable. The first \cref{sec.crane} arises in the exact physical setting of 2d linearly elastic structures introduced in \cref{sec.sopb}. We then turn to more complex physical situations; in \cref{sec.thermoelastic}, we consider the optimization of a two-dimensional thermoelastic device; in \cref{sec.liftdrag}, we deal with the three-dimensional optimization of an aerodynamic profile, and finally, in \cref{sec.heatTransfer}, we consider the optimization of a three-dimensional thermal-fluid heat exchanger. \subsection{Optimization of the shape of a two-dimensional crane}\label{sec.crane} Our first example deals with the optimal design of a 2d crane, as depicted on \cref{fig.cranesetting}. The considered shapes $\Omega$ are enclosed in a computational box $D$ with size $5 \times 5$; they are clamped on two small regions $\Gamma_D$ near the corners of their lower part, and uniform traction loads $g=(0,-1)$ are applied on the two regions $\Gamma_N$ located at the left and right ends of their upper part, accounting for the counterweight of the crane and the weight of the lifted object, respectively. The shapes are made of a linearly elastic material with normalized Lam\'e parameters $\lambda = 0.5769$ and $\mu=0.3846$. As presented in \cref{sec.sopb}, the physical behavior of the shape in this situation is described by the linearized elasticity system \cref{eq.elas}, and we solve the problem \[ \min_\Omega C(\Omega) \text{ s.t. } \Vol(\Omega) = V_T, \] where $C(\Omega)$ stands for the compliance \cref{eq.compliance} of $\Omega$, and where the target value $V_T$ for the volume $\Vol(\Omega)$ is set to $V_T = 2.5$. \begin{figure} \includegraphics[width=0.5\textwidth]{./figures/setcrane} \caption{Setting of the 2d crane optimization example of \cref{sec.crane}.}\label{fig.cranesetting} \end{figure} Several intermediate shapes are represented on \cref{fig.craneit}, exemplifying in particular the multiple topological changes undergone by the shape in the course of its evolution. The histories of the values of $C(\Omega$ and $\Vol(\Omega)$ are reported in \cref{fig.histocrane}. \begin{figure} \centering \begin{overpic}[width=0.45\textwidth]{figures/crane0} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/crane16} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \medskip \begin{overpic}[width=0.45\textwidth]{figures/crane25} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/crane50} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \medskip \begin{overpic}[width=0.45\textwidth]{figures/crane100} \put(0,3){\fcolorbox{black}{white}{$e$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/crane200} \put(0,3){\fcolorbox{black}{white}{$f$}} \end{overpic} \caption{(From (a) to (f)) Intermediate shapes obtained at iterations $0$, $16$, $25$, $50$, $100$ and $200$ of the crane optimization test case of \cref{sec.crane}}\label{fig.craneit} \end{figure} \begin{figure}[!ht] \includegraphics[width=0.8\textwidth]{./figures/histocrane} \caption{Evolution of the compliance and of the volume in the crane optimization example of \cref{sec.crane}.}\label{fig.histocrane} \end{figure} \subsection{Minimum compliance problem in thermoelasticity}\label{sec.thermoelastic} In this section, we illustrate the application of the level set based mesh evolution methodology to the minimization of the compliance of a thermoelastic structure. This test case is partially reproduced from the articles~\cite{feppon2020null, xia2008topology}. The considered shapes $\Omega$ are contained in the fixed computational domain $D=(0,2)\times (0,1)$; they are clamped on the left and right sides of $D$, whose reunion is denoted by $\Gamma_D$, and they are subjected to a traction load applied on a small region $\Gamma_N$ with size $0.0125$ centered at the middle of the bottom boundary, see \cref{fig.setthermoelas}. \begin{figure} \begin{center} \includegraphics[width=0.6\textwidth]{./figures/setthermoelas} \end{center} \caption{Physical setting of the compliance minimization problem of a thermoelastic structure considered in \cref{sec.thermoelastic}, issued from~\cite{xia2008topology}.}\label{fig.setthermoelas} \end{figure} The shape $\Omega$ is filled by a thermoelastic material, characterized by the Lam\'e parameters $\lambda=11510$, $\mu=7673$, and the thermal expansion coefficient $\alpha=0.77$ with reference temperature $T_{\mathrm{ref}}=0$. The quantities $\alpha$ and $T_{\mathrm{ref}}$ express that the structure experiences an additional stress induced by thermal dilation when the temperature $T$ inside the material is larger than $T_{\mathrm{ref}}$. In this example, the constant temperature field $T=5$ is imposed to the structure, causing expansion of the structure. The displacement $u_\Omega$ of the shape in this setting is the solution to the following linear thermoelasticity system: \begin{equation}\label{eqn:linearElasticity} \left\{ \begin{aligned} -\dv(\sigma(u_\Omega,T))&=0&& \In\Omega\\ u_\Omega&=0&& \On\Gamma_D\\ \sigma(u_\Omega,T) n&=g&& \On\Gamma_N\\ \sigma(u_\Omega,T) n&=0 &&\On\Gamma, \end{aligned} \right. \end{equation} where the stress tensor $\sigma(u,T)$ reads \begin{equation} \sigma(u,T)=Ae(u)-\alpha(T-T_{\textrm{ref}})\I \quad\text{with}\quad Ae =2\mu e +\lambda\tr(e)\I. \end{equation} We solve the following minimum compliance minimization problem: \[ \min_{\Omega \subset D} \quad J(\Omega):=\int_{\Omega} A e(u_\Omega):e(u_\Omega)\: \dx \quad \text{ s.t.} \quad \Vol(\Omega) \leq V_{T}, \] where the target volume is $V_{T}=0.7$. The shape derivative of the compliance in the present context has been calculated in~\cite{feppon2020null, xia2008topology}, and the formulas are omitted for brevity. The optimization results are shown on \cref{fig:history_thermo,fig:optimized_thermo}, which include respectively plots of a few intermediate shapes, and the mesh of the optimized design obtained at convergence. \begin{figure} \centering \begin{overpic}[width=0.31\textwidth]{figures/2Dit_0} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/2Dit_5} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/2Dit_10} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/2Dit_30} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/2Dit_50} \put(0,3){\fcolorbox{black}{white}{$e$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/2Dit_100} \put(0,3){\fcolorbox{black}{white}{$f$}} \end{overpic} \caption{(From (a) to (f)) Intermediate shapes obtained at iterations $0$, $5$, $10$, $30$, $50$ and $100$ of the thermoelastic test case of \cref{sec.thermoelastic}}\label{fig:history_thermo} \end{figure} \begin{figure} \centering \begin{overpic}[width=0.4\textwidth]{figures/2Dit_200} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.4\textwidth]{figures/2Dit_200_mesh} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \caption{(a) Optimized design (iteration $200$) in the thermoelastic test case of \cref{sec.thermoelastic}, (b) associated mesh ${\mathcal T}$ of the computational domain $D$.}\label{fig:optimized_thermo} \end{figure} \subsection{Lift--drag topology optimization for three-dimensional aerodynamic design}\label{sec.liftdrag} This section presents an application of the level set based mesh evolution method in the context of 3d aerodynamic design. We aim to optimize the shape of a flying solid obstacle immersed into an incoming fluid flow, so that it generate the largest possible lift force while keeping friction forces small. This classical industrial problem, which is commonly referred to as ``lift-drag'' optimization in the literature, has been mostly addressed in the context where the optimized design is described by a set of CAD parameters, see for instance~\cite{epstein2009comparative, gray2019open, Jameson, jameson1993computational, lambe2016matrix, mohammadi2010applied, pironneau1973optimum}. The test case presented in this section is a variation from that discussed in~\cite[Section~5.3]{feppon2020topology}. The physical setting is that depicted on \cref{fig:liftSetting}: the fixed computational domain is the cube $D=(0,1)^{3}$. It is filled by a fluid, entering from the left-hand boundary $\Gamma_{\In}$ with a velocity $v_0=(1,0,0)$, and exiting with a stress-free boundary condition at the opposite face $\Gamma_{\text{out}}$. The domain is divided into a fluid part $\Omega$, denoted by $\Omega_f$, and a solid obstacle $\Omega_s := D \setminus \overline{\Omega_f}$ which is such that $\Omega_s \Subset D$. The boundary of the solid obstacle is denoted by $\Gamma=\partial\Omega_s$. The physical state of the fluid is determined by its velocity field $v_\Omega$ and pressure field $p_\Omega$ which are characterized as the solutions to the system of incompressible steady-state Navier--Stokes equations: \begin{equation}\label{eqn:NavierStokes} \left\{ \begin{aligned} -\dv(\sigma_f(v_\Omega,p_\Omega))+\rho\nabla v_\Omega\,v_\Omega &=0 && \In\Omega_f\\ \dv(v_\Omega) &=0&&\In\Omega_f\\ v_\Omega &= v_0&& \On \Gamma_{\In}\\ \sigma_f(v_\Omega,p_\Omega) n &=0&&\On \Gamma_{\text{out}}\\ v_\Omega &= 0&&\On\Gamma \\ v_\Omega \cdot n &= 0 && \On \partial D\backslash (\overline{\Gamma_{\In}} \cup \overline{\Gamma_{\text{out}}}), \end{aligned} \right. \end{equation} where $\nu$ is the fluid viscosity, $\rho$ its density and where $\sigma_f(v,p)$ denotes the fluid stress tensor which is related to the rate of strain tensor $ e(v):= \frac{1}{2}(\nabla v+\nabla v^T)$ via Newton's law: \begin{equation*} \sigma (v,p):=2\nu e(v)-p\I, \end{equation*} and where $n$ denotes the unit normal vector to $\partial \Omega_f$, pointing outward $\Omega_f$. \begin{figure} \centering \begin{tikzpicture} \draw node at (-2.5,0.3) {{\color{blue}$\Gamma_{\text{in}}$}}; \draw node at (3,0) {{\color{red}$\Gamma_{\text{out}}$}}; \draw node at (2.5,-1.7) {{$D$}}; \draw node at (0.1,0.8) {{$\Omega_s$}}; \draw node at (0.3,-1.2) {{$\Omega_f$}}; \draw node at (0,0) {\includegraphics[height=4.5cm]{figures/dessin_lift}}; \end{tikzpicture} \caption{Physical setting of the lift-drag optimal design problem considered in \cref{sec.liftdrag}. }\label{fig:liftSetting} \end{figure} The considered problem is to maximize the lift generated by the obstacle $\Omega_s$ (which is the force allowing it to fly), under an upper bound constraint on the drag (i.e. the reaction force of the fluid medium on $\Omega_s$), so that the obstacle remains aerodynamic. In addition, the volume $\Vol(\Omega_s)$ and the center of mass $X(\Omega_s) \in \R^3$ of the obstacle are prescribed. The corresponding optimization program reads: \begin{equation} \begin{aligned} \min_\Omega & \qquad -\lift(\Omega)\\ \text{s.t.} &\begin{dcases} \drag(\Omega) \leqslant \drag_0\\ \Vol(\Omega_s) = V_T\\ X(\Omega_s) := \frac{1}{|\Omega_s|}\int_{\Omega_s} x \D x = x_0, \end{dcases} \end{aligned} \end{equation} where the lift functional $\lift(\Omega)$ accounts for the vertical force exerted by the fluid on the solid obstacle: \[ \lift(\Omega):=-\int_{\Gamma}\e_y\cdot \sigma_f(v_\Omega,p_\Omega) n\D s, \] where $\e_y=(0,1,0)$. The Drag functional $\drag(\Omega)$ is defined by \[ \drag(\Omega):=\int_{\Omega_f} \sigma_f(v_\Omega,p_\Omega):\nabla v_\Omega \D x= \int_{\Omega_f}2\nu e(v_\Omega):e(v_\Omega)\D x. \] The computation of the shape derivatives of these functionals is fairly classical and can be found e.g. in~\cite{dapogny2017geometrical, feppon2020, henrot2010optimal}. The parameters considered for the present test case are set to $\nu=1/200$, $\rho=1$, corresponding to a Reynolds number of the order of $\Reop \simeq 200$, $V_T:=0.03$ and $x_0=0$. The maximum value allowed for the drag is set to $\drag_0:=0.078405$, corresponding to $1.5$ times the minimum possible drag value for the prescribed solid volume (obtained from another minimization). This computationally intense problem is solved with parallel computing using 20 CPUs and the methodology described in~\cite{feppon2020topology}. The typical edge length imposed on the fluid-solid interface equals $0.01$. The final mesh features 154,982 vertices and 846,138 tetrahedra, including 134,637 vertices and 700,520 tetrahedra in the fluid domain. A few intermediate shapes obtained during the optimization procedure are shown on \cref{fig:intermediate_lift}. The optimized design, its mesh as well as its surrounding velocity field can be appraised on \cref{fig:cutview_lift,fig:lift_views}. Despite the low Reynolds number which makes this study quite unrealistic with respect to concrete aeronautic applications, is it quite remarkable that it is able to predict a design whose longitudinal profile recalls that of an airfoil, and whose ``V-shape'' at the back reminds that of airplanes. \begin{figure} \centering \begin{overpic}[width=0.31\textwidth]{figures/it_0l.png} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/it_1.png} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/it_2.png} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \medskip \begin{overpic}[width=0.31\textwidth]{figures/it_3.png} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/it_4.png} \put(0,3){\fcolorbox{black}{white}{$e$}} \end{overpic} \begin{overpic}[width=0.31\textwidth]{figures/it_6.png} \put(0,3){\fcolorbox{black}{white}{$f$}} \end{overpic} \caption{(From (a) to (f)) Intermediate designs obtained at iterations $0$, $5$, $10$, $50$, $100$ and $200$ of the optimization process for the lift-drag test case of \cref{sec.liftdrag}.}\label{fig:intermediate_lift} \end{figure} \begin{figure} \centering \includegraphics[height=6cm]{figures/it_200_cut.png} \caption{Optimized design in the lift-drag test case of \cref{sec.liftdrag}, with a cut view of the magnitude of the surrounding velocity field $v_\Omega$.}\label{fig:cutview_lift} \end{figure} \begin{figure} \centering \resizebox{\textwidth}{!}{\begin{tabular}{ccc} \begin{minipage}{0.37\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/it_200_z.png} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \end{minipage}& \begin{minipage}{0.30\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/it_200_top.png} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{minipage}& \begin{minipage}{0.30\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/it_200_below.png} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \end{minipage} \\[5em] \begin{minipage}{0.35\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/it_200_back.png} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{minipage}& \begin{minipage}{0.31\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/it_200_x.png} \put(0,3){\fcolorbox{black}{white}{$e$}} \end{overpic} \end{minipage}& \begin{minipage}{0.31\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/it_200_mesh.png} \put(0,3){\fcolorbox{black}{white}{$f$}} \end{overpic} \end{minipage} \end{tabular}} \caption{Different views of the optimized design obtained in the left-drag test case of \cref{sec.liftdrag}; (a) side view, (b) bottom view, (c) top view, (d) back view, (e) front view and (f) perspective view.}\label{fig:lift_views} \end{figure} \subsection{Optimal design of a three-dimensional fluid heating device} \label{sec.heatTransfer} The numerical example of this section is reproduced from~\cite[Section~5.5]{feppon2020topology}; it deals with the optimal design of a heating device whose purpose is to warm up a cold fluid under the coupled effects of convection and conduction. The situation is that depicted on \cref{fig:heatSetting}. The computational domain $D \subset \R^3$ is the box $[0,2]\times [0,1]\times [0,1] $. It is divided into a fluid phase $\Omega \subset D$, denoted by $\Omega_f$, and a solid phase $\Omega_s = D \setminus \overline{\Omega}$, with respective thermal conductivities $k_f$ and $k_s$; these are separated by the interface $\Gamma$. A ``cold'' fluid with density $\rho$ and viscosity $\nu$ enters from the left side of $D$ through a disk-shaped inlet boundary $\partial\Omega_f^{\text{in}}$ of radius $0.1$ with a unit parabolic velocity profile $v=v_0$ and a prescribed temperature $T=0$. The fluid absorbs the heat diffused from a ``hot'' stripe $\partial \Omega_T^D$ at the middle of $\partial D$ which is maintained at the constant temperature $T=10$, before exiting through a similar disk-shaped outlet boundary $\partial\Omega_f^{\text{out}}$ at the right side of $\partial D$ with a stress-free boundary condition. No heat is lost through the remaining parts of $\partial D$ which are considered adiabatic. \begin{figure} \centering \begin{tikzpicture} \draw node at (-4,0) {$T=0$}; \draw node at (-2.5,0.5) {$\partial\Omega_f^{\text{in}}$}; \draw node at (2.5,1.3) {$\partial\Omega_f^{\text{out}}$}; \draw node at (2.4,-1.7) {$T=10$}; \draw node at (0,0) {\includegraphics[height=4.5cm]{figures/drawing}}; \end{tikzpicture} \caption{Setting of the optimal heat transfer problem of \cref{sec.heatTransfer}. A ``cold'' fluid flows from the left to the right blue disk-shaped regions and captures the heat communicated by the red stripe $\partial \Omega^D_T$ where the ``hot'' temperature $T=10$ is imposed.}\label{fig:heatSetting} \end{figure} In this context, we aim to optimize the shape $\Omega_f = \Omega$ of the fluid phase, and thereby the repartition of $\Omega_f$ and $\Omega_s$ within $D$, in order to maximize the total heat carried by the fluid, under constraints on the static pressure loss through the device and on the volume of $\Omega_f$. Mathematically, the following optimization program is considered, see~\cite{feppon2020,feppon2020topology}: \begin{equation}\label{eqn:mt25v} \begin{aligned} \min_{\Omega} & \quad J(\Omega):=-\int_{\Omega_f} \rho c_p v_\Omega \cdot\nabla T_\Omega \D x,\\ \text{s.t.} & \quad \left\{ \begin{aligned} \text{DP}(\Omega):=\int_{\partial\Omega_{f}^{\text{in}}} p_\Omega \D s -\int_{\partial\Omega_{f}^{\text{out}}} p_\Omega \D s & \leqslant \text{DP}_{T},\\ \Vol(\Omega_f) &= V_{T}. \end{aligned}\right. \end{aligned} \end{equation} Here, the velocity $v_\Omega$ and the pressure $p_\Omega$ of the fluid are governed by the static, incompressible Navier--Stokes equations: \begin{equation}\label{eqn:NavierStokes2} \left\{ \begin{aligned} -\dv(\sigma_f(v_\Omega,p_\Omega))+\rho\nabla v_\Omega \,v_\Omega & =0 && \In\Omega_f, \\ \dv(v_\Omega)& =0&&\In\Omega_f,\\ v_\Omega& = v_0&& \On\partial\Omega_f^{\text{in}},\\ \sigma_f(v_\Omega,p_\Omega) n &=0&&\On\partial\Omega_f^{\text{out}},\\ v_\Omega &= 0&&\On\Gamma. \end{aligned}\right. \end{equation} where $\sigma_f(v,p)$ is the stress tensor within a fluid with velocity $v$ and pressure $p$, which is related to the rate of strain tensor $ e(v):= \frac{1}{2}(\nabla v+\nabla v^T)$ via the Newton's law: \begin{equation*} \sigma_f(v,p):=2\nu e(v)-p\I. \end{equation*} The temperature field $T_\Omega$ inside $D$ depends on the velocity $v_\Omega$ of the fluid via the static convection-diffusion equation: \begin{equation}\label{eqn:heatTransfer} \left\{ \begin{aligned} -\dv(k_f\nabla T_\Omega)+\rho c_p v_\Omega \cdot\nabla T_\Omega& =0&& \In\Omega_f,\\ -\dv(k_s\nabla T_\Omega)& = 0&&\In\Omega_s,\\ T_\Omega&=0&& \On\partial\Omega_{f}^{\text{in}},\\ T_\Omega &=10&& \On\partial\Omega_{T}^{D},\\ -\partialn{T_\Omega}&=0&&\On \text{ other boundaries},\\ T_\Omega &\text{ is continuous}&&\On\Gamma,\\ \text{The flux } -k\partialn{T_\Omega } &\text{ is continuous} && \On\Gamma, \end{aligned}\right. \end{equation} where $c_p$ is the thermal capacity of the fluid. The values of the physical parameters entering the formulation \cref{eqn:mt25v,eqn:NavierStokes,eqn:heatTransfer} are chosen as follows. The density and the viscosity of the fluid are respectively $\rho=10$ and $\nu=0.167$, so that the Reynolds number of the flow equals $\Reop\simeq 60$. Its thermal conductivity and capacity coefficients are respectively set to $k_f =1$ and $c_p=500$; hence, the P\'eclet number $\mathrm{Pe}$ equals $5,000$. The thermal conductivity of the solid phase is $k_s=10$. The imposed upper bound on the pressure drop in \cref{eqn:mt25v} is set to $\text{DP}_{T}=0.85$, while the target volume for the fluid phase equals 20\% of the total volume: $V_{T}=0.2\Vol(D)$. The numerical realization of this program relies on the application of the level set based mesh evolution \cref{algo.bf}, where the formulas for the derivatives of the shape functionals at play have been established in~\cite{feppon2018shape}. The initial repartition of the phases $\Omega_f$, $\Omega_s$ features islands of solid spherical inclusions, see \cref{fig:convheatIters}$\MK$(a). A few intermediate shapes arising in the course of the optimization process can be seen on \cref{fig:convheatIters}, illustrating the dramatic topological changes undergone by the system. The resulting optimized design is displayed on \cref{fig:convheatOpt}. Interestingly, the latter features thin inclusions of fluid attached to the main pipes so as to take advantage of the insulating effect induced by the low thermal conductivity of the fluid with respect to that of the solid. \begin{figure} \centering \begin{overpic}[width=0.45\textwidth]{figures/it_0.png} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/it_5.png} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/it_15.png} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/it_25.png} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/it_50.png} \put(0,3){\fcolorbox{black}{white}{$e$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/it_152.png} \put(0,3){\fcolorbox{black}{white}{$f$}} \end{overpic} \caption{(From (a) to (f)) Iterations 0, 5, 15, 25, 50 and 152 of the optimization process of a heating device considered in \cref{sec.heatTransfer}.}\label{fig:convheatIters} \end{figure} \begin{figure} \centering \begin{overpic}[width=0.55\textwidth]{figures/heatex.png} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.44\textwidth]{figures/fluid.png} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/solid.png} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/fluid_cut.png} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \caption{(a) Optimized heating device obtained in the example of \cref{sec.heatTransfer}; (b) mesh of the fluid component of the final design; (c) sectional view of the optimized solid domain (d) sectional view of the optimized fluid domain. In all the pictures, the colors correspond to the temperature profile.}\label{fig:convheatOpt} \end{figure} The sizes of the meshes of $D$ considered for this test-case range from approximately 22,000 to 278,100 vertices. The numerical resolutions of the Navier--Stokes and the advection-diffusion equations \cref{eqn:NavierStokes,eqn:heatTransfer} are efficiently achieved by combining the finite element method with domain decomposition techniques on 30 processes. The total computation takes approximately 2 days. \section{Conclusion and perspectives} In this article, we have presented a numerical framework for shape and topology optimization which combines the level set method with remeshing techniques to reconcile the antagonistic needs for accurate mechanical computations on the shape and for a robust description of its evolution. This numerical strategy has been successfully applied in a wide variety of physical contexts, including those of linearly elastic structures, fluid-structure interacting devices, heat exchangers, etc. It could also be used fruitfully beyond the realm of shape optimization, and notably in the simulation of multi-phase problems (for instance, in the study of multi-phase flows), or in the solution of inverse problems implying the detection or reconstruction of objects. Beyond ever more complicated physical situations, we are currently working on extending this method to deal with the optimization of the shape and topology of regions embedded in a fixed surface. This would allow, for instance, to optimize the regions of the boundary of a domain bearing various types of boundary conditions, in the spirit of our previous work~\cite{dapogny2020optimization}. We also wish to extend this methodology to handle the evolution of open curves in 2d, and open surfaces in 3d, which would make it possible to optimize shell structures. \subsection*{Acknowledgements} %The work of C.D. is partially supported by the project ANR-18-CE40-0013 SHAPO, financed by the French Agence Nationale de la Recherche (ANR). The body of work presented in this article owes much to the diverse expertise of many colleagues: in particular, we are very grateful to M.~Albertelli, F.~Bordeu, J.~Cortial, C.~Dobrzynski, A.~Froehly, G.~Michailidis. The authors are especially indebted to G.~Allaire and P.~Frey, who were at the inception of this project, and made it thrive with their key insight and relentless suggestions. \appendix \section{Manual of the companion code}\label{app.man} In this appendix, we describe the \texttt{python}-based implementation of the level set mesh evolution method accompanying this article. Our aim is to provide a code which meets a compromise between simplicity, re-usability and computational efficiency: we hope that it is simple enough so that a user can easily get acquainted with its main features, and at the same time general and modular enough so that it can be modified with a minimum amount of effort to address a whole gamut of new and more challenging applications. The source code can be downloaded from the \texttt{Github} repository \[ \texttt{\url{https://github.com/dapogny/sotuto}} \] which is organized in several folders: the one named \texttt{./base} contains the material devoted to the benchmark 2d cantilever test case; it serves as support for this presentation. The other folders are as many variations around this main theme, which are briefly presented in \cref{sec.2gofurther}. \CD{This manual is organized as follows. After a brief description of a guiding test case in \cref{sec.descexman}, we provide in \cref{app.maninstall} the information for installing the needed external libraries in our framework. We then present the overall architecture of the program in \cref{sec.archi}. In \cref{sec.manmain}, we outline the practical shape optimization strategy realizing the abstract \cref{algo.bf}, and we introduce the functions corresponding to its main stages.} These functions are detailed in the subsequent sections, with a particular emphasis on the parts that the user may have to modify in order to implement his own examples. Note that, for clarity, we have allowed some redundancy between the various sections of this manual, as it is not necessarily meant to be read in linear fashion. \subsection{Description of the 2d cantilever test-case}\label{sec.descexman} Our model example arises in the linearized elasticity setting described in \cref{sec.sopb}. The considered shape $\Omega$ represents a cantilever beam; $\Omega$ is comprised in a box $D$ with size $2 \times 1$, it is clamped on the left-hand side $\Gamma_D$ of $\partial D$, and traction loads $g= (0,-1)$ are applied on a given region $\Gamma_N$ near the center of its right-hand side, see \cref{fig.settingmanual}. \begin{figure} \begin{center} \includegraphics[width=0.75\textwidth]{./figures/setcanti} \end{center} \caption{Setting of the 2d cantilever test-case considered in \cref{app.man}.}\label{fig.settingmanual} \end{figure} In this context, we aim to solve the problem \begin{equation}\label{eq.sopbman} \min_\Omega \: C(\Omega), \quad\text{ s.t.}\quad \Vol(\Omega) = V_T, \end{equation} where the compliance $C(\Omega)$ and the volume $\Vol(\Omega)$ of a shape $\Omega$ have been defined in \cref{sec.sopb} as: \[ C(\Omega) = \int_\Omega A e(u_\Omega) : e(u_\Omega) \:\dx \quad\text{ and}\quad \Vol(\Omega) = \int_\Omega \:\dx, \] and where $V_T$ is a volume target. The elastic displacement $u_\Omega$ involved in the definition of $C(\Omega)$ is the unique solution in $H^1_{\Gamma_D}(\Omega)^d$ to the linear elasticity system \cref{eq.elas}, whose variational formulation reads: \begin{equation}\label{eq.varfelas} \forall v \in H^1_{\Gamma_D}(\Omega)^d, \quad \int_\Omega Ae(u_\Omega) : e(v) \:\dx = \int_{\Gamma_N} g \cdot v \:\ds. \end{equation} For the convenience of the reader, we recall from \cref{sec.hadamard} the expressions of the shape derivatives of $C(\Omega)$ and $\Vol(\Omega)$: \begin{equation}\label{eq.derCpman} C^\prime(\Omega)(\theta) = -\int_\Gamma Ae(u_\Omega) : e(u_\Omega) \:\theta\cdot n \:\ds, \quad\text{ and}\quad \Vol^\prime(\Omega)(\theta) = \int_\Gamma \theta\cdot n \:\ds. \end{equation} \subsection{Getting started: download and installation of the files}\label{app.maninstall} As a preliminary stage, the use of our program requires the installation of four open-source libraries, written in \texttt{C} or \texttt{C++}: \begin{itemize} \item The library \texttt{mshdist}~\cite{dapogny2012computation} performs the computation of the signed distance function to an input domain at the vertices of a simplicial computational mesh, as detailed in \cref{sec.sdf}. It can be freely downloaded from the repository \begin{center} \texttt{\url{https://github.com/ISCDtoolbox/Mshdist}} \end{center} and compiled by following the instructions supplied at this address. \item The library \texttt{advection}~\cite{bui2012accurate} solves the advection equation \cref{eq.advls} governing the evolution of a level set function on a simplicial computational mesh in 2d and 3d, see \cref{sec.advect}. It is available at the following link: \begin{center} \texttt{\url{https://github.com/ISCDtoolbox/Advection}} \end{center} where instructions are given to compile the library. \item The library \texttt{mmg}~\cite{dapogny2014three} is in charge of the remeshing operations of the strategy, as presented in \cref{sec.remesh}. It can be downloaded through the webpage of the project \begin{center} \texttt{\url{https://www.mmgtools.org/}} \end{center} where a series of tutorials is available to exemplify the main features of this program, see also the \texttt{github} link \begin{center} \texttt{\url{https://github.com/MmgTools/mmg}} \end{center} \item The software \texttt{FreeFem}~\cite{hecht2012new} is a numerical environment where partial differential equations can be solved in a few lines of commands through the input of their variational formulation. More generally, it conveniently and efficiently allows to handle finite element spaces and functions thanks to an intuitive \texttt{C++} based pseudo-language. The latest version of \texttt{Freefem} can be downloaded from the webpage\footnote{A \texttt{python} package interfacing \texttt{FreeFem} is available at the following link: \texttt{\url{https://gitlab.com/florian.feppon/pyfreefem}} }: \begin{center} \texttt{\url{https://freefem.org/}} \end{center} which also contains an exhaustive documentation, and multiple worked examples. \end{itemize} In addition to these four required codes, we recommend the installation of the open-source vizualization software \texttt{medit}, provided at the following address: \begin{center} \texttt{\url{https://github.com/ISCDtoolbox/Medit}} \end{center} where a short manual is also available\footnote{The \texttt{python} package \texttt{PyMedit} can be used to conveniently handle and vizualize the mesh structures processed by \texttt{medit}: \texttt{\url{https://gitlab.com/florian.feppon/pymedit}}}. Our numerical implementation of the level set based mesh evolution method relies on a main frame written in \texttt{python}, supplemented by \texttt{FreeFem} scripts dedicated to the operations involving finite element computations. The calls to the four aforementioned external libraries are encapsulated in \texttt{python} functions which merely execute \texttt{unix} command lines via subprocesses\footnote{These command lines are tailored to \texttt{Linux} and \texttt{MacOs} environments, and some adaptations may be necessary for systems relying on \texttt{Windows}.}. These interface functions are not supposed to be modified by the user, and we limit ourselves with presenting their syntax in \cref{sec.manmmgmshdist} below. Their proper functioning requires that, after installation of the libraries \texttt{mshdist}, \texttt{advection}, \texttt{mmg} and \texttt{FreeFem}, the user set the paths of the corresponding executables in the {\texttt{\modu{path.py}}} file, as described in \cref{sec.pathpy}. The communications between the main \texttt{python} program and these codes (such as the exchange of input parameters or output results) are realized via an exchange file referred to as \texttt{path.EXCHFILE} in the \texttt{python} program; this mechanism should be invisible to the user. \CD{\begin{rema} In principle, there is no obstruction in replacing \texttt{FreeFem++} with another, e.g. commercial finite element solver in the proposed implementation. Essentially, this task would require minor adaptations to convert the meshes processed by \texttt{mmg} (and the other aforementioned libraries) into a suitable format for the latter. \end{rema}} \subsection{Global architecture of the program}\label{sec.archi} The source code in the \texttt{./base} folder is made of files of three different types: \begin{itemize} \item The \texttt{.py} files are \texttt{python} modules; \item The \texttt{.edp} files are \texttt{FreeFem} scripts, which are executed by calling \texttt{FreeFem} from \texttt{python} functions; \item The \texttt{.idp} files are also processed by \texttt{FreeFem}; they contain general information and routines shared by most of the \texttt{.edp} scripts. \end{itemize} Two of the above files are dedicated to the input of user-defined parameters: \begin{itemize} \item The module {\texttt{\modu{path.py}}} contains global variables and the paths to the directories where input and output objects are stored, to the executables, etc. For instance, \texttt{path.EXCHFILE} is the address of the aforementioned exchange file between \texttt{python} and \texttt{FreeFem} modules, the file \texttt{path.HISTO} contains the data (values of the compliance and volume) of the successive shapes produced during the optimization process, \texttt{path.FREEFEM} is the command line for calling \texttt{FreeFem}, etc. This file also contains the definition of shorthands used throughout the program: \texttt{path.step(n,"mesh")} is the string of characters containing the name of the mesh of $D$ at the $n^{\text{th}}$ iteration, \texttt{path.step(n,"u.sol")} refers to the file containing the elastic displacement, etc. The contents of this file are more thoroughly described in \cref{sec.pathpy}. \item The file \texttt{macros.idp} contains global variables and macros used in \texttt{FreeFem} scripts; it is described in \cref{sec.macrosidp}. \end{itemize} The main program is executed by calling the file {\texttt{\modu{main.py}}}, which orchestrates the functions from the other \texttt{python} modules. These modules are: \begin{itemize} \item The module {\texttt{\modu{inout.py}}} contains the functions for reading and printing data from and to external files; \item The module {\texttt{\modu{inigeom.py}}} is devoted to the creation of the initial shape, see \cref{sec.inigeom} for its description; \item The module {\texttt{\modu{mechtools.py}}} gathers the functions in charge of the mechanical calculations of the framework, and notably the resolution of the linear elasticity system \cref{eq.varfelas}, the calculation of the compliance of shapes, that of a descent direction for the problem \cref{eq.sopbman}, etc. These are described in \cref{sec.tutoelas,man.descent}; \item The module {\texttt{\modu{lstools.py}}} contains the functions dedicated to the treatment of level set functions, and notably the interface functions with the \texttt{C} libraries \texttt{mshdist} and \texttt{advection}. The practical use of these functions is described in \cref{sec.manmmgmshdist}. \item The module {\texttt{\modu{mshtools.py}}} contains the interface functions with the library \texttt{mmg}, see \cref{sec.manmmgmshdist}; \end{itemize} The execution of our program mainly results in the production of data files of two types, written in the common formats used by \texttt{mshdist}, \texttt{advect}, \texttt{FreeFem}, \texttt{mmg} and \texttt{medit}: \begin{itemize} \item \texttt{.mesh} files contain meshes, \item \texttt{.sol} files contain (scalar- or vector-valued) $\mathbb{P}_1$ finite element functions given by their values at the vertices of a companion mesh. \end{itemize} For brevity, and since the structure of such files is not needed for our purpose, we do not describe them in the present manual, referring for this matter to the documentation of any of the codes mentioned above. \subsection{Overview of the global strategy: description of the file \ {\tt main.py}}\label{sec.manmain} The code contained in the folder \texttt{./base} is executed by entering the following command line from the root folder. \begin{lstlisting}[style=shell,escapechar=|] python./base/main.py \end{lstlisting} In the present section, we describe in some detail how the file \texttt{main.py} implements the strategy of \cref{algo.bf}; the outline of this practical realization is provided in \cref{algo.bfprac}. \begin{algorithm}[ht] \caption{Practical realization of the level set based mesh evolution method.}\label{algo.bfprac} \begin{algorithmic}[0] \STATE \textbf{Initialization:} \noindent \begin{itemize} \item Function {\texttt{\modu{inout}}}\texttt{.iniWF}: Initialization of folders and exchange files. \item (Optional) Function {\texttt{\modu{inout}}}\texttt{.testLib}: Test of the links with \texttt{python} functions and external \texttt{C} libraries. \item Function {\texttt{\modu{inigeom}}}\texttt{.iniGeom}: Creation of the initial mesh $\calT^0$ of $D$. \item Function {\texttt{\modu{mechtools}}}\texttt{.elasticity}: Calculation of $u_{\Omega^0}$. \item Functions {\texttt{\modu{mechtools}}}\texttt{.compliance} and {\texttt{\modu{mechtools}}}\texttt{.volume}: Calculation of $C(\Omega^0)$ and $\Vol(\Omega^0)$. \end{itemize} \FOR{$n=0,\dots,\texttt{path.MAXIT}-1$} \STATE \begin{enumerate}\relabelitem \item Function {\texttt{\modu{lstools}}}\texttt{.mshdist}: Calculation of the signed distance function $\phi^n$ on $\calT^n$. \item Functions {\texttt{\modu{mechtools}}}\texttt{.gradCp} and {\texttt{\modu{mechtools}}}\texttt{.gradV}: Calculation of the gradients of $C(\Omega)$ and $\Vol(\Omega)$ on $\calT^n$. \item Functions {\texttt{\modu{mechtools}}}\texttt{.descent}: Calculation of a descent direction for \cref{eq.sopbman}. \item Functions {\texttt{\modu{mechtools}}}\texttt{.merit}: Calculation of the merit of $\Omega^n$. \item \textbf{Line search:} \textbf{for} $k=0, \ldots, \texttt{path.MAXITLS}-1$ \textbf{do} \begin{enumerate}\relabeliitem \item Function {\texttt{\modu{lstools}}}\texttt{.advect}: Resolution of the level set advection equation, for the new function $\phi^{n+1}$ on the old mesh $\calT^n$ of $D$. \item Function {\texttt{\modu{mshtools}}}\texttt{.mmg}: Creation of the new mesh $\calT^{n+1}$ of $D$, where the $0$ level set of $\phi^{n+1}$ is explicitly discretized. \item Function {\texttt{\modu{mechtools}}}\texttt{.elasticity}: Calculation of $u_{\Omega^{n+1}}$ on $\calT^{n+1}$. \item Functions {\texttt{\modu{mechtools}}}\texttt{.compliance}, {\texttt{\modu{mechtools}}}\texttt{.volume} and {\texttt{\modu{mechtools}}}\texttt{.merit}: Calculation of $C(\Omega^{n+1})$, $\Vol(\Omega^{n+1})$ and $M(\Omega^{n+1})$. \item Decision: Acceptation / Rejection of the new iterate. \end{enumerate} \item[] \textbf{end for} \end{enumerate} \ENDFOR \RETURN Mesh $\calT^n$ of $D$ where $\Omega^n$ is discretized as a submesh $\caliT^n$. \end{algorithmic} \end{algorithm} \subsubsection{Initialization} The program starts with the creation of the folder \texttt{./res/}, where all the mesh and solution files produced during the execution will be saved. Meanwhile, the files \texttt{path.EXCHFILE}, \texttt{path.HISTO} are created, as well as other exchanges files whose descriptions are not needed by the user. \begin{lstlisting}[style=python,escapechar=|] # Initialize folders and exchange files inout.iniWF() \end{lstlisting} Then, the calls to the four external libraries \texttt{FreeFem}, \texttt{mshdist}, \texttt{advection} and \texttt{mmg} are tested. \begin{lstlisting}[style=python,escapechar=|] # Test the links with external C libraries inout.testLib() \end{lstlisting} This should result with the following terminal output. \begin{lstlisting}[style=shell,escapechar=|] FreeFem installation working. Mshdist installation working. Advect installation working. Mmg installation working. All external libraries working. \end{lstlisting} Of course, once the user is confident with his settings, this testing stage can be deactivated by simply removing this instruction. The next call to the function \texttt{iniGeom} from the module {\texttt{\modu{inigeom.py}}}, which is described in \cref{sec.inigeom}, is meant to create the initial mesh $\calT^0$ of $D$, equipped with the boundary conditions of the test-case, and featuring an explicit discretization of the initial shape $\Omega^0$ as a submesh $\caliT^0$. \begin{lstlisting}[style=python,escapechar=|] # Creation of the initial mesh inigeom.iniGeom(path.step(0,"mesh"),path.step(0,"phi.sol") \end{lstlisting} The resulting mesh file \texttt{step.0.mesh} is stored in the folder \texttt{path.RES}; it is depicted in \cref{fig.cantiman}$\MK$(a), where the interior submesh $\caliT^0$ is that made of the black triangles (identified by the reference \texttt{path.REFINT}), while the exterior submesh $\caleT^0$ is the collection of white triangles (with reference \texttt{path.REFEXT}). The linear elasticity system \cref{eq.varfelas} for the elastic displacement $u_{\Omega^0}$ of $\Omega^0$ is then solved on the interior submesh $\caliT^0$ thanks to the function \texttt{elasticity} from the {\texttt{\modu{mechtools.py}}} module. \begin{lstlisting}[style=python,escapechar=|] # Resolution of the state equation mechtools.elasticity(path.step(0,"mesh"),path.step(0,"u.sol")) \end{lstlisting} This function is described in \cref{sec.tutoelas} below. Finally, the values $C(\Omega^0)$ and $\Vol(\Omega^0)$ of the compliance and the volume of $\Omega^0$ are calculated owing to the following commands: \begin{lstlisting}[style=python,escapechar=|] # Calculation of the compliance and the volume of the shape newCp = mechtools.compliance(path.step(0,"mesh"), path.step(0,"u.sol")) newvol = mechtools.volume(path.step(0,"mesh")) \end{lstlisting} The functions \texttt{compliance} and \texttt{volume} from the module {\texttt{\modu{mechtools.py}}} are detailed in \cref{sec.tutoelas}. \subsubsection{Main optimization loop} At the beginning of each iteration $n$, ranging from $0$ to \texttt{path.MAXIT} $-1$, a few shorthands are defined, such as for instance: \begin{itemize} \item \texttt{curmesh} is the address \texttt{path.step(n,"mesh")} (corresponding, by default, to the file \texttt{step.n.mesh} in the folder \texttt{path.RES}) for the mesh $\calT^n$ of $D$; \item \texttt{curphi} is the address \texttt{path.step(n,"phi.sol")} (corresponding, by default, to the file \texttt{step.n.phi.sol} in \texttt{path.RES}) for the level set function $\phi^n : D \to \R$. \end{itemize} At this moment, the following objects are available: \begin{itemize} \item The mesh $\calT^n$ of the domain $D$, called \texttt{curmesh}; the shape $\Omega^n$ (resp. the exterior part $D \setminus \overline{\Omega^n}$) is represented by the submesh $\caliT^n$ (resp. $\caleT^n$), whose elements are the triangles with label \texttt{path.REFINT} (resp. \texttt{path.REFEXT}). \item The displacement $u_{\Omega^n}$, solution to the linear elasticity equation \cref{eq.varfelas}) on $\Omega^n$, is called \texttt{curu}. It is encoded as a $\mathbb{P}_1$ finite element function on the whole mesh \texttt{curmesh}, with the convention that only its values in the interior part of \texttt{curmesh} are to be considered. \item The respective values \texttt{newCp}, \texttt{newvol} of the compliance and volume~$C(\Omega^n)$ and $\Vol(\Omega^n)$ of the current shape $\Omega^n$. \end{itemize} The $n^{\text{th}}$ iteration starts with the calculation of the signed distance function $\phi^n = d_{\Omega^n}$ to $\Omega^n$ at the vertices of the mesh $\calT^n$. This relies on the function \texttt{mshdist} from the module {\texttt{\modu{lstools.py}}}: \begin{lstlisting}[style=python,escapechar=|] # Generation of a level set function for $\Omega^n$ on $D$ lstools.mshdist(curmesh,curphi) \end{lstlisting} The use of this function is described in \cref{sec.manmmgmshdist}. The shape gradients of the objective and constraint functionals $C(\Omega)$ and $\Vol(\Omega)$ are then calculated. This is realized thanks to the functions \texttt{gradCp} and \texttt{gradV} from the module {\texttt{\modu{mechtools.py}}}, which rely in particular on the Hilbertian extension-regularization procedure presented in \cref{sec.descentth}. \begin{lstlisting}[style=python,escapechar=|] # Calculation of the gradients of compliance and volume mechtools.gradCp(curmesh,curu,curCPgrad) mechtools.gradV(curmesh,curVgrad) \end{lstlisting} These functions are described in \cref{sec.tutogcpv}. A descent direction for the optimization problem \cref{eq.sopbman} is then inferred from these gradients, following the constrained optimization strategy presented in \cref{sec.nullspacealgo}. \begin{lstlisting}[style=python,escapechar=|] # Calculation of a descent direction mechtools.descent(curmesh,curphi,curCp,curCpgrad,curvol, curVgrad,curgrad) \end{lstlisting} This function is described in \cref{sec.tutodescent}. It takes as arguments \begin{itemize} \item A string of characters \texttt{curmesh} for the name of the current mesh of $D$; \item The name \texttt{curphi} of the level set function for the shape $\Omega^n$; \item The values \texttt{curCp}, \texttt{curvol} of the objective (compliance) and constraint (volume) functionals; \item The respective names \texttt{curCpgrad}, \texttt{curVgrad} for the gradients of these functionals; \item The name \texttt{curgrad} for the file where the calculated descent direction should be saved. \end{itemize} \begin{rema}\label{rem.mannormdescent} The descent direction produced by this function is normalized, so that its $L^\infty$ norm equals the typical size {\texttt{path.MESHSIZ}} of an element in the mesh, see \cref{sec.tutodescent} about this practice. Hence, a motion of the shape $\Omega^n$ in the direction of \texttt{curgrad} for a pseudo-time step $\tau^n=1$ corresponds to a maximum amplitude \texttt{path.MESHSIZ} for a motion of a point of the shape. \end{rema} The above function \texttt{descent} prints in the file \texttt{path.EXCHFILE} the values of the coefficients $\lambda$, $S$ (which boils down to a scalar value in the present example), $\alpha_J$, $\alpha_G$ used in the optimization strategy, as they are needed in the evaluation of the merit function $M(\Omega)$, see \cref{eq.merit}. This value of $M(\Omega^n)$ is then calculated, thanks to the following function, which is described in \cref{sec.manevalmerit}. \begin{lstlisting}[style=python,escapechar=|] # Evaluation of the merit of $\Omega^n$ merit = mechtools.merit(curCp,curvol) \end{lstlisting} A line search procedure then begins, under the form of a nested loop indexed by $k$; for each substep, the following operations are carried out. \begin{enumerate} \item The motion of the current shape $\Omega^n$ given by \texttt{curmesh} and \texttt{curphi} in the direction \texttt{curgrad} is realized, for a certain pseudo-time \texttt{coef}. Because of the chosen normalization for \texttt{curgrad} (see \cref{rem.mannormdescent}), the maximum displacement of a point of the shape equals \texttt{coef} $*$ \texttt{path.MESHSIZ}. \begin{lstlisting}[style=python,escapechar=|] # Advection of the level set function print(" Level Set advection") lstools.advect(curmesh,curphi,curgrad,coef,newphi) \end{lstlisting} The use of the function {\texttt{\modu{lstools}}}.\texttt{advect} is explained in \cref{sec.manmmgmshdist}. The level set function $\phi^{n+1}$ for the new shape $\Omega^{n+1}$ resulting from this operation, stored in the file \texttt{newphi}, is still attached to the mesh \texttt{curmesh}. \item The $0$ level set $\phi^{n+1}$, named \texttt{newphi}, is explicitly discretized in the mesh \texttt{curmesh} by calling \texttt{mmg}. \begin{lstlisting}[style=python,escapechar=|] # Creation of a mesh associated to the new shape print(" Local remeshing") retmmg = mshtools.mmg2d(curmesh,1,newphi,path.HMIN, path.HMAX,path.HAUSD,path.HGRAD,1,newmesh) \end{lstlisting} The function {\texttt{\modu{mshtools}}}.\texttt{mmg2d} in charge of the interface between the library \texttt{mmg} and our \texttt{python} framework is described in \cref{sec.manmmgmshdist}. \item \CD{The elastic displacement $u_{\Omega^{n+1}}$ is calculated on the new shape defined by the mesh $\calT^{n+1}$ of $D$, named \texttt{newmesh}, or more accurately, on the submesh $\caliT^{n+1}$ of $\calT^{n+1}$ made of the triangles with reference \texttt{path.REFINT}. The corresponding values of the compliance and volume functionals, and then the new value of the merit are evaluated. These operations are executed unless \texttt{mmg} has accidentally failed, in which case the result of the current iteration of the inner, line search loop is immediately rejected (see below).} \begin{lstlisting}[style=python,escapechar=|] if ( retmmg ) : # Resolution of the state equation on the new shape print(" Resolution of linearized elasticity system") mechtools.elasticity(newmesh,newu) # Calculation of the new values of compliance and volume print(" Evaluation of the new merit") newCp = mechtools.compliance(newmesh,newu) newvol = mechtools.volume(newmesh) newmerit = mechtools.merit(newCp,newvol) \end{lstlisting} \end{enumerate} The final step of the inner line-search loop is the decision stage: the new shape \texttt{newmesh} is accepted if the corresponding value \texttt{newmerit} of the merit function is less than that of the previous shape \texttt{curmesh} (up to some tolerance), or if a maximum number of iterations has been made in the line search procedure with this descent direction. Otherwise, the new iterate is refused; the algorithm goes back to the beginning of the line search procedure, with a reduced value of the ``time-step'' \texttt{coef}. \begin{lstlisting}[style=python,escapechar=|] # Accept iteration: break and increase slightly the "time step" if ( retmmg and ( newmerit < merit + path.TOL*abs(merit) ) or ( k == 2 ) or ( coef < path.MINCOEF ) ) : coef = min(path.MAXCOEF,1.1*coef) print(" Iteration {} - subiteration {} accepted\n".format(n,k)) break # Reject iteration: go to start of line search # with a decreased "time step" else : print(" Iteration {} - subiteration {} rejected".format(n,k)) proc = subprocess.Popen(["rm {nmesh}".format(nmesh=newmesh)], shell=True) proc.wait() coef = max(path.MINCOEF,0.6*coef) \end{lstlisting} \subsubsection{Outputs}\label{sec.outputs} At each iteration $n=0,\ldots,$ \texttt{path.MAXIT}$-1$ of the main loop of our optimization procedure, the following files are saved in the \texttt{path.RES} repository: \begin{itemize} \item The current mesh $\calT^n$ of $D$, referred to via the shortcut \texttt{curmesh}, is saved as the file \texttt{path.step(n,"mesh")} (which is, unless specified otherwise by the user, the file \texttt{step.n.mesh} in the folder \texttt{path.RES}); \item The level set function $\phi^n$ for the shape $\Omega^n$, called \texttt{curphi}, is saved as the file \texttt{path.step(n,"phi.sol")}; \item The solution $u_{\Omega^n}$ to the elasticity system associated to $\Omega^n$ is saved in the file \texttt{path.step(n,"u.sol")}; let us recall that, although $u_{\Omega^n}$ is only defined on $\Omega^n$ from the mathematical point of view, it is actually stored as a $\mathbb{P}_1$ finite element function on the total mesh $\calT^n$, being understood that only the values at the vertices of the interior part $\caliT^n$ are relevant; \item The gradients \texttt{curCpgrad} and \texttt{curVgrad} of the objective and constraint functionals are saved in the files \texttt{path.step(n,"CP.grad.sol")} and \texttt{path.step(n,"V.grad.sol")}, \item The descent direction $\theta^n\!$, called \texttt{curgrad}, is saved in the file \texttt{path.step(n,"grad.sol")}. \end{itemize} Eventually, the values of the compliance and volume attached to each shape $\Omega^n$are stored in the \texttt{path.HISTO} file. A few shapes produced in the course of the optimization are depicted in \cref{fig.cantiman}. The histories describing the evolution of the compliance and volume of the shape can be displayed by calling the {\texttt{\modu{plot.py}}} module via the following command line: \begin{lstlisting}[style=shell,escapechar=|] python base/plot.py \end{lstlisting} see \cref{fig.histocantiman}. \begin{figure} \begin{center} \begin{overpic}[width=0.49\textwidth]{figures/canti0} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/canti20} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.49\textwidth]{figures/canti50} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/canti150} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{center} \caption{Iterations (a) $0$, (b) $20$, (c) $50$ and (d) $150$ (final) of the optimization process in the cantilever test-case discussed in \cref{app.man}.}\label{fig.cantiman} \end{figure} \begin{figure}[!ht] \centering \includegraphics[width=0.8\textwidth]{./figures/histocanti} \caption{Evolution of the compliance and the volume functionals in the cantilever example of \cref{app.man}.}\label{fig.histocantiman} \end{figure} \subsection{Description of the two files dedicated to the specification of the global variables and paths to executables} \subsubsection{The file {\texttt{\modu{path.py}}}}\label{sec.pathpy} The file {\texttt{\modu{path.py}}} contains the global variables used in all the \texttt{python} functions. \begin{lstlisting}[style=python,escapechar=|] # Global parameters REFDIR = 1 # Reference for Dirichlet B.C REFNEU = 2 # Reference for Neumann B.C. REFISO = 10 # Reference for the boundary edges # of the shape REFINT = 3 # Reference of the triangles # in the interior domain REFEXT = 2 # Reference of the triangles # in the exterior domain \end{lstlisting} The parameters regarding the desired size of the elements in the mesh are specified next; their meanings are described in \cref{sec.manmmgmshdist}. \begin{lstlisting}[style=python,escapechar=|] # Parameters of the mesh MESHSIZ = 0.02 HMIN = 0.001 HMAX = 0.02 HAUSD = 0.001 HGRAD = 1.3 # Other parameters EPS = 1e-10 # Precision parameter EPSP = 1e-20 # Precision parameter for packing MAXIT = 150 # Maximum number of iterations # in the shape optimization process MAXITLS = 3 # Maximum number of iterations # in the line search procedure [...] \end{lstlisting} The file {\texttt{\modu{path.py}}} also contains the command lines needed to call the external programs \texttt{mshdist}, \texttt{advect}, \texttt{mmg} and \texttt{FreeFem} from the terminal; as hinted at in \cref{app.maninstall}, these instructions should be supplied by the user after their installation. \begin{lstlisting}[style=python,escapechar=|] # Call for the executables of external codes FREEFEM = "FreeFem++ -nw" MSHDIST = "mshdist" ADVECT = "/Users/dapogny/Advection/build/advect" MMG2D = "/Users/dapogny/mmg/build/bin/mmg2d_O3" \end{lstlisting} The names of the various \texttt{Freefem} scripts of the archive used in the \texttt{python} functions are provided: \begin{lstlisting}[style=python,escapechar=|] # Path to FreeFem scripts FFTEST = SCRIPT + "testFF.edp" FFDESCENT = SCRIPT + "descent.edp" FFEVALOBJ = SCRIPT + "evalobj.edp" ... \end{lstlisting} Eventually, a shorthand is defined for the names of various output files at each iteration of the process, see \cref{sec.outputs}. \begin{lstlisting}[style=python,escapechar=|] # Shortcut for various file types def step(n,typ) : return STEP + "." + str(n) + "." + typ \end{lstlisting} \subsubsection{Global variables and macros related to the \texttt{Freefem} functions}\label{sec.macrosidp} The archive contains two files \texttt{macros.idp} and \texttt{inout.idp} which are loaded in the preamble of all the \texttt{FreeFem} scripts, via the following lines. \begin{lstlisting}[style=ff,escapechar=|] include "./sources/inout.idp" include "./sources/macros.idp" \end{lstlisting} The file \texttt{macros.idp} contains the global variables and the macros shared by the \texttt{FreeFem} scripts of the archive. It is organized as follows: \begin{lstlisting}[style=ff,escapechar=|] /* File for communication of data with python */ string EXCHFILE = "./res/exch.data"; /* Inner product for extension / regularization */ real alpha = getrParam(EXCHFILE,"Regularization"); macro psreg(u,v) ( int2d(Th)( alpha^2*(dx(u)*dx(v)+dy(u)*dy(v)) + u*v) ) // EOM /* Linear elasticity parameters */ real lm = 0.5769; real mu = 0.3846; /* Load case */ real loadx = 0.0; real loady = -1.0; \end{lstlisting} In particular, this file contains the definition of the inner product $a(\,\cdot\,,\cdot\,)$ endowing the space of deformations $\theta$ with a Hilbertian structure. In the present context, $a(\,\cdot\,,\cdot\,)$ is a bilinear form, associating to two scalar-valued functions $u,v : D \to \R$ the value \begin{equation}\label{eq.aman} a(u,v) = \alpha^2 \int_D \nabla u \cdot \nabla v \:\dx + \int_D uv \:\dx, \end{equation} where $\alpha$ can be interpreted as a regularization length, see \cref{sec.descentth}. The file \texttt{inout.idp} gathers several functions in charge of reading and printing parameters or files. These functions ensure the communication with the \texttt{python} part of the implementation; their syntaxes are meant to be explicit enough so that the user does not need to enter their implementation. \begin{lstlisting}[style=ff,escapechar=|] /* Get real parameter from file */ func real getrParam(string file,string kwd) { [...] } [...] /* Read a .sol file containing a scalar-valued solution */ func int loadsol(string sin, real[int] & u) { [...] } /* Read a .sol file containing a vector-valued solution in 2d */ func int loadvec2(string sin, real[int] & ux, real[int] & uy) { [...] } /* Save scalar function u as a .sol file */ func int printsol(string sout, real[int] & u) { [...] } /* Save vector field [ux,uy] as a .sol file */ func int printvec2(string sout, real[int] & ux, real[int] & uy) { [...] } \end{lstlisting} \subsection{Creation of the initial geometry and specification of the test case}\label{sec.inigeom} The function dedicated to the creation of the mesh $\calT^0$ of $D$, in which the initial shape $\Omega^0$ is explicitly discretized as a submesh $\caliT^0$, is contained in the module {\texttt{\modu{inigeom.py}}}. \begin{lstlisting}[style=python,escapechar=|] def iniGeom(mesh) : \end{lstlisting} This function proceeds in several stages: \begin{enumerate} \item The exchange file \texttt{path.EXCHFILE} is filled with the names of the output mesh and of the intermediate level set function used in the creation of the latter: \begin{lstlisting}[style=python,escapechar=|] # Fill in exchange file inout.setAtt(file=path.EXCHFILE,attname="MeshName", attval=mesh) inout.setAtt(file=path.EXCHFILE,attname="PhiName", attval=path.TMPSOL) \end{lstlisting} \item A mesh of the computational domain $D$ is generated via the call to the \texttt{FreeFem} script \texttt{path.FFINIMSH}, whose contents are described in the next listing. \begin{lstlisting}[style=ff,escapechar=|] /* Get mesh and sol names */ string MESH = getsParam(EXCHFILE,"MeshName"); int REFDIR = getiParam(EXCHFILE,"Dirichlet"); int REFNEU = getiParam(EXCHFILE,"Neumann"); /* Create mesh */ /* Mesh definition */ border left(t=0.0,1.0){x=0.0; y=1.0-t; label=REFDIR;}; border bot(t=0.0,2.0){x=t; y=0.0; label=0;}; border right1(t=0.0,0.45){x=2.0; y=t; label=0;}; border right2(t=0.45,0.55){x=2.0; y=t; label=REFNEU;}; border right3(t=0.55,1.0){x=2.0; y=t; label=0;}; border top(t=0.0,2.0){x=2.0-t; y=1.0; label=0;}; mesh Th = buildmesh(left(100)+bot(200)+right1(45) +right2(10)+right3(45)+top(200)); /* Save mesh */ savemesh(Th,MESH); \end{lstlisting} \item The resulting mesh is modified towards quality improvement, thanks to the external library \texttt{mmg}. \begin{lstlisting}[style=python,escapechar=|] # Call to mmg2d for remeshing the background mesh mshtools.mmg2d(mesh,0,None,path.HMIN,path.HMAX, path.HAUSD,path.HGRAD,0,mesh) \end{lstlisting} Again, the syntax of the function {\texttt{\modu{mshtools}}}.\texttt{mmg2d} in charge of making the interface with this external library and our \texttt{python} implementation is described in detail in \cref{sec.manmmgmshdist}. \item In this example, $\Omega^0$ is the total computational domain $D$, deprived from a collection of holes, see \cref{fig.cantiman}$\MK$(a). In order to create the associated mesh $\calT^0$ of $D$, one level set function $\phi^0$ for this shape is generated via the \texttt{FreeFem} script \texttt{path.FFINILS}. \begin{lstlisting}[style=ff,escapechar=|] /* Get mesh and sol names */ string MESH = getsParam(EXCHFILE,"MeshName"); string PHI = getsParam(EXCHFILE,"PhiName"); /* Read mesh */ mesh Th = readmesh(MESH); /* Finite element space and functions */ fespace Vh(Th,P1); Vh phi; /* Definition of the initial level set function */ func real iniLS(real xx,real yy) { [...] return (dd); } phi = -iniLS(x,y); /* Save LS function */ printsol(PHI,phi[]); \end{lstlisting} \item This level set function $\phi^0$ is explicitly discretized in the mesh by another call to \texttt{mmg}, resulting in the mesh displayed in \cref{fig.cantiman}$\MK$(a). \end{enumerate} \subsection{External calls to the libraries \texttt{mshdist}, \texttt{advection} and \texttt{mmg}}\label{sec.manmmgmshdist} This section is devoted to the \texttt{python} functions through which the external libraries \texttt{mshdist}, \texttt{advection} and \texttt{mmg} are executed in our implementation. Since they are not meant to be modified by the user, we solely present their syntax. We recall that, after the compilation of these three codes, the command lines for their executions have to be specified in the file {\texttt{\modu{path.py}}}, as described in \cref{app.maninstall,sec.pathpy}. Let $D$ be a domain equipped with a mesh $\calT$ in which the shape $\Omega$ of interest is discretized, as the submesh $\caliT$ of triangles with labels \texttt{path.REFINT}. The calculation of the signed distance function $\phi = d_\Omega$ to $\Omega$ is carried out by the function \begin{lstlisting}[style=python,escapechar=|] def mshdist(mesh,phi) : \end{lstlisting} from the module {\texttt{\modu{lstools.py}}}. Here, \begin{itemize} \item \texttt{mesh} is a string of characters, bearing the name of the \texttt{.mesh} file for $\calT$. \item \texttt{phi} is a string of characters, containing the name of the \texttt{.sol} file where $\phi$ should be saved. \end{itemize} When $D$ is a domain equipped with a mesh $\calT$, the resolution of the level set advection equation \cref{eq.advls}, starting from $\phi^0$, with velocity field $V : D \to \R^d$, over a time period $(0,T)$, is realized via the function \begin{lstlisting}[style=python,escapechar=|] def advect(mesh,phi,vel,step,newphi) : \end{lstlisting} from the module {\texttt{\modu{lstools.py}}}. Here, \begin{itemize} \item \texttt{mesh} is a string of characters, containing the name of the \texttt{.mesh} file for $\calT$. \item \texttt{phi} is a string of characters referring to the \texttt{.sol} file for the initial level set function $\phi^0$. \item \texttt{vel} is a string of characters, containing the \texttt{.sol} file for the velocity field $V$. \item \texttt{step} is a real value, standing for the length $T$ of the advection time period. \item \texttt{newphi} is a string of characters, containing the name of the \texttt{.sol} file where the resulting level set function should be written. \end{itemize} Last but not least, the remeshing operations of an input mesh $\calT$ of a domain $D$ involved in our framework are realized by the \texttt{mmg} library. All the calls to this library are encapsulated in the proposed \texttt{python} implementation via the function \begin{lstlisting}[style=python,escapechar=|] def mmg2d(mesh,ls,phi,hmin,hmax,hausd,hgrad,nr,out) : \end{lstlisting} defined in the module {\texttt{\modu{mshtools.py}}}. The parameters of this function are the following: \begin{itemize} \item \texttt{mesh} is a string of characters bearing the name of the input mesh $\calT$. \item \texttt{ls} is an integer, taking the values $0$ or $1$. When \texttt{ls}$=0$, the mesh $\calT$ is remeshed towards improving the quality of its elements; when \texttt{ls}$=1$, the $0$ level set of the function \texttt{phi} is explicitly discretized in $\calT$, which is then remeshed towards improving the quality of its elements. \item \texttt{phi} is a string of characters, bearing the name of the level set function in the case where the \texttt{ls} parameter is set to $1$. \item \texttt{hmin} is a real value, corresponding to the minimum authorized length for an edge in the mesh. \item \texttt{hmax} is a real value, corresponding to the maximum authorized length for an edge in the mesh. \item \texttt{hausd} is a real value, encoding the maximum gap authorized between an edge on the surface part $\calS_{\calT}$ of $\calT$ and the underlying curve drawn on the ideal surface, see \cref{sec.remeshgen} and \cref{fig.bezier}. \item \texttt{hgrad} is a gradation parameter controlling the shocks of lengths between adjacent edges in the mesh. It is a real value, specifying the maximum ratio allowed between the lengths of two adjacent edges (typical values are $1.3$ or $1.4$). \item \texttt{nr} is an integer taking two values $0$, or $1$. It equals $1$ when sharp angles are to be detected by \texttt{mmg}, then treated as corners during remeshing, $0$ otherwise. \item \texttt{out} is a string of characters containing the name of the output mesh. \end{itemize} The function \texttt{mmg2d} returns $1$ if the remeshing process is successful and $0$ if an error occurred. \subsection{Mechanical calculations}\label{sec.tutoelas} Let $D$ be a computational domain, equipped with a mesh $\calT$ in which the shape $\Omega$ of interest is explicitly discretized, as the submesh $\caliT$ composed of the triangles with label \texttt{path.REFINT}. The resolution of the linear elasticity system for the displacement $u_\Omega$ is achieved by the \texttt{elasticity} function, pertaining to the module {\texttt{\modu{mechtools.py}}}. \begin{lstlisting}[style=python,escapechar=|] def elasticity(mesh,u) : \end{lstlisting} In this command, \begin{itemize} \item \texttt{mesh} is a string of characters, bearing the name of the mesh $\calT$. \item \texttt{u} is a string of characters for the name of the output vector field $u_\Omega$. The latter is treated as a $\mathbb{P}_1$ finite element vector field, stored as a \texttt{.sol} file containing values at all the vertices of $\calT$, being understood that only those values at the vertices of $\caliT$ will be considered. \end{itemize} The function proceeds by first entering this information in the exchange file \texttt{path.EXCHFILE}, then calling \texttt{Freefem} with the script \texttt{path.FFELAS}, that we now touch on. After a few lines, where notably the mesh \texttt{Th} of $D$ is loaded, this mesh is truncated thanks to the \texttt{trunc} function. \begin{lstlisting}[style=ff,escapechar=|] /* Mesh of the interior part and corresponding FE space */ mesh Thi = trunc(Th,(reg(x,y) == REFINT),label=REFINT); fespace Vhi(Thi,P1); Vhi uix,uiy,vix,viy; \end{lstlisting} Hence, \texttt{Thi} is the desired mesh $\caliT$ of the shape $\Omega$ and \texttt{Vhi} is the corresponding finite element space of $\mathbb{P}_1$ Lagrange finite element functions. Then, the linear elasticity system is solved on this interior mesh. \begin{lstlisting}[style=ff,escapechar=|] /* Variational formulation of the problem */ problem elas([uix,uiy],[vix,viy]) = int2d(Thi)(mu*(2.0*dx(uix)*dx(vix) + (dx(uiy)+dy(uix))*(dx(viy)+dy(vix)) + 2.0*dy(uiy)*dy(viy)) + lm*(dx(uix)+dy(uiy))*(dx(vix)+dy(viy))) - int1d(Thi,REFNEU)(loadx*vix+loady*viy) + on(REFDIR,uix=0.0,uiy=0.0); /* Solve problem */ elas; \end{lstlisting} Finally, the result is transferred onto the total mesh of $D$. \begin{lstlisting}[style=ff,escapechar=|] /* Transfer the problem on the full mesh */ ux = uix; uy = uiy; \end{lstlisting} The calculation of the compliance of $\Omega$ is easily realized from the datum of the mesh \texttt{mesh} and of the resulting elastic displacement \texttt{u}. In our implementation, this task relies on the function \begin{lstlisting}[style=python,escapechar=|] def compliance(mesh,u) : \end{lstlisting} from the module {\texttt{\modu{mechtools.py}}}. This function essentially calls the \texttt{FreeFem} script \texttt{path.FFCPLY}, which is organized as follows. \begin{lstlisting}[style=ff,escapechar=|] /* Get mesh and sol names, and global parameters */ string MESH = getsParam(EXCHFILE,"MeshName"); string SOL = getsParam(EXCHFILE,"DispName"); int REFINT = getiParam(EXCHFILE,"Refint"); int REFNEU = getiParam(EXCHFILE,"Neumann"); int REFISO = getiParam(EXCHFILE,"ReferenceBnd"); /* Loading mesh */ mesh Th = readmesh(MESH); [...] /* Read solution */ loadvec2(SOL,ux[],uy[]); /* Calculate compliance */ cply = int2d(Th,REFINT)(mu*(2.0*dx(ux)*dx(ux) + (dx(uy)+dy(ux))*(dx(uy)+dy(ux)) + 2.0*dy(uy)*dy(uy)) + lm*(dx(ux)+dy(uy))*(dx(ux)+dy(uy))); /* Save result */ setrParam(EXCHFILE,"Compliance",cply); \end{lstlisting} The calculation of the volume of $\Omega$ is likewise conducted via the function \begin{lstlisting}[style=python,escapechar=|] def volume(mesh) : \end{lstlisting} from the module {\texttt{\modu{mechtools.py}}}, whose syntax is analogous. \subsection{Calculation of a descent direction}\label{man.descent} In this section, we detail the calculation of shape gradients for the compliance and volume shape functionals, and that of a descent direction $\theta$ for the optimization problem \cref{eq.sopbman}. The shape $\Omega$ of interest is supplied by a mesh $\calT$ of the computational domain $D$ in which it is explicitly discretized, as the submesh $\caliT$ of triangles with reference \texttt{path.REFINT}. Let us recall from \cref{sec.descentth,sec.nullspacealgo} that these calculations hinge on the choice of an inner product acting on the vector space of deformation fields, according to the so-called Hilbertian extension-regularization procedure. The macro for this inner product is supplied in the \texttt{macros.idp} file, as described in \cref{sec.macrosidp}. \subsubsection{Calculation of the gradient of the compliance (and volume) functional}\label{sec.tutogcpv} The calculation of the shape gradient for the compliance functional is realized by the function \begin{lstlisting}[style=python,escapechar=|] def gradCp(mesh,disp,grad) : \end{lstlisting} from the {\texttt{\modu{mechtools.py}}} module. This function essentially calls \texttt{FreeFem} with the file \texttt{path.FFGRADCP}, which is organized as follows. \begin{lstlisting}[style=ff,escapechar=|] /* Get mesh and sol names */ [...] /* Loading mesh */ mesh Th = readmesh(MESH); /* Finite element spaces and functions */ fespace Vh(Th,P1); fespace Vh0(Th,P0); Vh ux,uy,g,v; Vh0 reg = region; /* Load elastic displacement */ loadvec2(DISP,ux[],uy[]); /* Mesh of the interior part and corresponding FE spaces */ mesh Thi = trunc(Th,(reg(x,y) == REFINT),label=REFINT); fespace Vhi(Thi,P1); fespace Vh0i(Thi,P0); Vhi uix,uiy; Vh0i Aeueu; /* Integrand of the shape derivative */ uix = ux; uiy = uy; Aeueu = mu*(2.0*dx(uix)*dx(uix) + (dx(uiy)+dy(uix))*(dx(uiy)+dy(uix)) + 2.0*dy(uiy)*dy(uiy)) + lm*(dx(uix)+dy(uiy))*(dx(uix)+dy(uiy)); /* Resolution of the extension - regularization problem */ problem velext(g,v) = psreg(g,v) - int1d(Thi,REFISO)(-Aeueu*v) + on(REFNEU,g=0.0); /* Solve problem and save solution */ velext; printsol(GRADCP,g[]); \end{lstlisting} In a few words, the variational problem featured in this script solves the identification problem \cref{eq.Rieszgradient} with the inner product $a(\,\cdot\,,\cdot\,)$ --~given by \cref{eq.aman} and supplied in the file \texttt{macros.idp} described in \cref{sec.macrosidp}~-- and the right-hand side $C^\prime(\Omega)(\cdot n)$ given by \cref{eq.derCpman}. The calculation of the gradient of the volume functional follows the exact same trail, and we do not detail further the contents of the corresponding function \begin{lstlisting}[style=python,escapechar=|] def gradV(mesh,grad) : \end{lstlisting} \subsubsection{Derivation of a descent direction}\label{sec.tutodescent} The null-space optimization algorithm described in \cref{sec.nullspacealgo} is implemented in the\linebreak function \begin{lstlisting}[style=python,escapechar=|] def descent(mesh,phi,Cp,gCp,vol,gV,g) : \end{lstlisting} from the {\texttt{\modu{mechtools.py}}} module. The latter essentially calls \texttt{FreeFem} with the file \texttt{path.FFDESCENT}, which proceeds as follows. After reading the computational mesh $\calT$ of $D$, finite element spaces and parameters are declared. The finite element functions associated to the level set function $\phi$ for $\Omega$ and the gradients of the compliance and volume functionals are loaded. \begin{lstlisting}[style=ff,escapechar=|] loadsol(PHI,phi[]); loadsol(GRADCP,thJ[]); loadsol(GRADV,thG[]); \end{lstlisting} The contributions $\xi_J$ and $\xi_G$ to total descent direction $\theta$, given by \cref{eq.thetanoptim}, are then calculated. \begin{lstlisting}[style=ff,escapechar=|] /* Coefficients for the descent direction */ m = psreg(thG,thG); lambda = 1.0 / m * psreg(thJ,thG); /* Null space and range contributions to descent direction */ xiJ = thJ - lambda*thG; xiG = 1.0/m*(vol-vtarg)*thG; \end{lstlisting} The coefficients $\alpha_J$ and $\alpha_G$ are next calculated, and the scalar amplitude of the descent direction is inferred, see again \cref{eq.thetanoptim} and \cref{sec.descentth}. \begin{lstlisting}[style=ff,escapechar=|] maxxiJ = max(-xiJ[].min,xiJ[].max); mMaxxiJ = max(maxxiJ,maxNormXiJ); maxxiG = max(-xiG[].min,xiG[].max); alphaJ = AJ*meshsiz / (eps^2+mMaxxiJ); alphaG = AG*meshsiz / (eps^2+maxxiG); /* Scalar descent direction */ g = - alphaJ*xiJ - alphaG*xiG; \end{lstlisting} Then, the extended normal vector field $n$ to the boundary of the shape $\Omega$ is calculated on the whole domain $D$ via the formula $n(x) = \frac{\nabla\phi(x)}{|\nabla \phi(x)|}$, first as a $\mathbb{P}_0$ vector field on $\calT$, then as a $\mathbb{P}_1$ vector field. \begin{lstlisting}[style=ff,escapechar=|] /* Extended normal vector as a P0 function over the mesh */ norm0 = sqrt(eps+dx(phi)*dx(phi)+dy(phi)*dy(phi)); nx0 = dx(phi) / norm0; ny0 = dy(phi) / norm0; /* Extended normal vector as a P1 function over the mesh */ problem extnx(nx,v) = int2d(Th)(nx*v) - int2d(Th)(nx0*v); problem extny(ny,v) = int2d(Th)(ny*v) - int2d(Th)(ny0*v); extnx; extny; norm = sqrt(eps^2+nx*nx+ny*ny); nx = nx / norm; ny = ny / norm; \end{lstlisting} Eventually, a descent direction for the optimization problem \cref{eq.sopbman} is obtained as a vector field $D \to \R^2$ and the data are saved in the appropriate files. \begin{lstlisting}[style=ff,escapechar=|] /* Save solution */ printvec2(GRAD,gx[],gy[]); /* Save coefficients for the merit function */ setrParam(EXCHFILE,"Lagrange",lambda); setrParam(EXCHFILE,"Penalty",m); setrParam(EXCHFILE,"alphaJ",alphaJ); setrParam(EXCHFILE,"alphaG",alphaG); setrParam(EXCHFILE,"NormXiJ",maxxiJ); \end{lstlisting} \subsubsection{Evaluation of the merit function}\label{sec.manevalmerit} The evaluation of the merit $M(\Omega)$ of a shape $\Omega$, defined in \cref{eq.merit}, is simply realized via the function \texttt{merit} of the {\texttt{\modu{mechtools}}} module: \begin{lstlisting}[style=python,escapechar=|] def merit(Cp,vol) : \end{lstlisting} This simple function takes as arguments the compliance \texttt{Cp} and the volume \texttt{vol} of the shape, which have been calculated via the functions defined in \cref{sec.tutoelas}. It also relies on the coefficients $\lambda$, $S$, $\alpha_J$, $\alpha_G$ produced by the optimization algorithm of \cref{man.descent}, stored in the exchange file \texttt{path.EXCHFILE}. \begin{lstlisting}[style=python,escapechar=|] # Read parameters in the exchange file [alphaJ] = inout.getrAtt(file=path.EXCHFILE,attname="alphaJ") [alphaG] = inout.getrAtt(file=path.EXCHFILE,attname="alphaG") [ell] = inout.getrAtt(file=path.EXCHFILE,attname="Lagrange") [m] = inout.getrAtt(file=path.EXCHFILE,attname="Penalty") [vtarg] = inout.getrAtt(file=path.EXCHFILE, attname="VolumeTarget") merit = alphaJ*(Cp - ell*(vol-vtarg)) + 0.5*alphaG/m*(vol-vtarg)**2 return merit \end{lstlisting} \subsection{Elaborations upon this code}\label{sec.2gofurther} The code described in this appendix has been prepared with the concern that it should be reasonably simple to adapt to the treatment of a wide variety of shape optimization problems. In this section, we briefly describe four such settings, illustrating as many features that can advantageously enrich the base code. The corresponding source files are provided in the \texttt{github} repository of the project. \begin{figure} \resizebox{\textwidth}{!}{\begin{tabular}{cc} \begin{minipage}{0.45\textwidth} %\hspace{-0.5cm} \begin{overpic}[width=1.0\textwidth]{figures/setbrtopo} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \end{minipage}& \begin{minipage}{0.45\textwidth} \begin{overpic}[width=1.0\textwidth]{figures/setlshape} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{minipage}\\ \begin{minipage}{0.59\textwidth} %\hspace{-0.5cm} \begin{overpic}[width=1.\textwidth]{figures/setmload} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \end{minipage}& \begin{minipage}{0.45\textwidth} \vspace{0.3cm} \begin{overpic}[width=1.0\textwidth]{figures/setcanti3d} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{minipage} \end{tabular}} \caption{(a) Setting of the bridge test case considered in \cref{man.topoder}; (b) setting of the L-shaped beam example considered in \cref{man.stress}; (c) setting of the multi-load bridge considered in \cref{man.mload}; (d) setting of the 3d cantilever test case considered in \cref{man.3d}.}\label{fig.moreexamples} \end{figure} \subsubsection{A combination of shape and topological derivatives: minimization of the compliance of a \texorpdfstring{$2d$}{2d} bridge} \label{man.topoder} This section illustrates the resolution of a shape and topology optimization problem by means of a combination of the boundary variation \cref{algo.bf} with an occasional use of topological derivatives, as originally proposed in~\cite{allaire2005structural}. As we have mentioned in \cref{rem.dertopo}, the topological derivative $\dJ_T(\Omega)(x)$ of a function $J(\Omega)$ at some point $x \in \Omega$ is the first non trivial term in the asymptotic expansion \[ J(\Omega \setminus \overline{B(x,r)}) = J(\Omega) + r^d \dJ_T(\Omega)(x) + \orm(r^d); \] in particular, the value of $J(\Omega)$ is decreased if a tiny hole is nucleated around a point $x \in \Omega$ where $\dJ_T(\Omega)(x)$ is negative. Let us recall the following expressions of the topological derivatives of the compliance and volume functionals $C(\Omega)$ and $\Vol(\Omega)$: \begin{multline}\label{eq.topderCV} \dC_T(\Omega)(x) = \frac{\pi(\lambda + 2\mu)}{2\mu(\lambda+\mu)} \Big(4\mu Ae(u_\Omega) : e(u_\Omega) + (\lambda - \mu) \tr(Ae(u_\Omega)) \tr(e(u_\Omega)) \Big)(x),\\ \text{and}\quad \dVol_T(\Omega)(x) = -\pi, \end{multline} see e.g.~\cite{garreau2001topological,novotny2012topological}. The considered example deals with the optimization of a 2d bridge, as depicted on \cref{fig.moreexamples}$\MK$(a). The shapes $\Omega$ of interest are contained in a box $D$ with size $2 \times 1.5$; they are clamped on their lower-left corner, their vertical displacement is prevented on their lower-right corner, and a unit vertical load $g = (0,-1)$ is applied on a small region $\Gamma_N$ around the middle of their bottom side. In this context, we aim to minimize the compliance $C(\Omega)$ of $\Omega$ under a constraint on its volume $\Vol(\Omega)$: we solve problem \cref{eq.sopbman} with a target volume~$V_T = 0.7$. To achieve this, we rely on the boundary variation \cref{algo.bfprac} with an extra ingredient: at each iteration $n$ of the main loop, a boolean value \texttt{dotopder} is calculated: \begin{lstlisting}[style=python,escapechar=|] # Equals 1 if topological derivative operation, 0 else dotopder = 1 if ( it % 5 == 0 and it <= 50 ) else 0 \end{lstlisting} The iterations $n$ where this parameter equals $0$ unfold exactly as those of the base code: the shape gradients of $C(\Omega)$ and $\Vol(\Omega)$ are calculated, a descent direction for the problem \cref{eq.sopbman} is inferred from the null space optimization algorithm described in \cref{sec.nullspacealgo,man.descent}, then the shape $\Omega^n$ is updated into the new iterate $\Omega^{n+1}$ (a level set function $\phi^n$ for $\Omega^n$ is calculated on the mesh $\calT^n$ of $D$, the level set advection equation is solved, etc.). The iterations $n$ where \texttt{dotopder} equals $1$ are associated to a stage of topological sensitivity analysis, which proceeds along the following lines. \begin{lstlisting}[style=python,escapechar=|] # Modification of the shape using topological gradients if ( dotopder ) : # Calculation of the topological derivatives of # the compliance and volume print(" Computation of topological derivatives") mechtools.topgradCp(curmesh,curu,curtopCpgrad) mechtools.topgradV(curmesh,curtopVgrad) # Calculation of the topological derivative of merit mechtools.descentTG(curmesh,curphi,curCp,curtopCpgrad,curvol, curtopVgrad,curtopgrad) # Update of the level set function of the shape lstools.creaHoles(curmesh,curphi,curtopgrad, path.VFRACTG,newphi) # Creation of a mesh associated to the new shape print(" Local remeshing") retmmg = mshtools.mmg2d(curmesh,1,newphi,path.HMIN,path.HMAX, path.HAUSD,path.HGRAD,1,newmesh) [...] # Decision [...] \end{lstlisting} Briefly, the topological derivatives of $C(\Omega)$ and $\Vol(\Omega)$ are calculated from the formulas \cref{eq.topderCV}, thanks to the functions \texttt{topgradCp} and \texttt{topgradV} of the {\texttt{\modu{mechtools.py}}} module. These essentially call the \texttt{FreeFem} scripts \texttt{path.FFTOPGRADCP} and \texttt{path.FFTOPGRADV}, whose syntaxes are quite straightforward. Then, the function \texttt{descentTG} of the {\texttt{\modu{mechtools.py}}} module infers a descent direction for \cref{eq.sopbman} as a scalar field on the computational domain, owing to a simple adaptation of the null space algorithm of \cref{sec.nullspacealgo}. Finally, the function \begin{lstlisting}[style=python,escapechar=|] def creaHoles(mesh,phi,grad,volfrac,newphi) : \end{lstlisting} of the {\texttt{\modu{lstools.py}}} module is called. The latter expects the following arguments: \begin{itemize} \item A string of characters \texttt{mesh} for the name of the current mesh of $D$. \item The name \texttt{phi} of the file containing a level set function for the shape $\Omega^n$. \item The name \texttt{grad} of the file containing the descent direction for \cref{eq.sopbman} resulting from the topological sensitivity analysis. \item A real parameter \texttt{volfrac} (taking values between $0$ and $1$). \item The name \texttt{newphi} of the output file containing the level set function of the new shape $\Omega^{n+1}$. \end{itemize} This function essentially relies on the \texttt{FreeFem} script \texttt{path.FFUPTG}. It first identifies by dichotomy the \texttt{volfrac} percentage of elements of the input shape \texttt{phi} where the values of \texttt{grad} are the most negative, and it then returns a new level set function \texttt{newphi} featuring holes in place of these elements. A few iterations of the resolution of \cref{eq.sopbman} with this strategy in the 2d bridge example are depicted on \cref{fig.brex}. The corresponding source code is contained in the following folder: \[ \texttt{\url{https://github.com/dapogny/sotuto/topder}} \] \begin{figure} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/br0} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/br26} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.45\textwidth]{figures/br40} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.45\textwidth]{figures/br150} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{center} \caption{Iterations (a) $0$, (b) $26$, (c) $40$ and (d) $150$ of the optimization of the shape of a 2d bridge with the combined use of shape and topological derivatives proposed in \cref{man.topoder}.}\label{fig.brex} \end{figure} \subsubsection{Handling shape functionals whose derivative involve an adjoint state: minimization of the stress within an L-shaped beam} \label{man.stress} The present example illustrates how the base code described in this appendix can be adapted to handle other objective functionals than the compliance, whose shape derivatives involve an adjoint state. Let us consider the physical situation depicted on \cref{fig.lb}. Shapes $\Omega$ are beams contained in an L-shaped domain $D$ with size $1 \times 1$; they are clamped on their upper side $\Gamma_D$, and a unit vertical load $g=(0,-1)$ is applied on a small region $\Gamma_N$ around the middle of their right-hand side. In this context, we consider the shape optimization problem \begin{equation}\label{eq.minS} \min_\Omega S(\Omega) \text{ s.t. } \Vol(\Omega) = V_T, \end{equation} \CD{where $S(\Omega)$ is the integral measure of the stress inside $\Omega$ defined in \cref{eq.stressfunc} with a weight function $k(x)$ equal to $1$ everywhere except in two small regions around $\Gamma_D$ and $\Gamma_N$, where it equals $0$. The volume target is set to $V_T = 0.7$.} The shape derivative $S^\prime(\Omega)(\theta)$ of $S(\Omega)$ is provided in \cref{eq.Sprime}; its expression involves the adjoint state $p_\Omega \in H^1_{\Gamma_D}(\Omega)^d$, which is defined as the solution to the system~\cref{eq.adjstress}. The numerical resolution of this problem is almost identical to that featured in the base implementation, except when it comes to the calculation of a shape gradient for the stress-based functional $S(\Omega)$. At each iteration $n$ of the main loop, this operation proceeds along the following lines: \begin{lstlisting}[style=python,escapechar=|] # Calculation of a descent direction print(" Computation of the adjoint state ") mechtools.adjoint(curmesh,curu,curp) # Calculation of the gradients of stress and volume print(" Computation of a descent direction") mechtools.gradS(curmesh,curu,curp,curSgrad) mechtools.gradV(curmesh,curVgrad) mechtools.descent(curmesh,curphi,curS,curSgrad,curvol, curVgrad,curgrad) \end{lstlisting} In short, the adjoint state $p_{\Omega^n}$ is calculated by the finite element method, thanks to the function \texttt{adjoint} of the {\texttt{\modu{mechtools.py}}} module; the latter consists in turn in executing the \texttt{FreeFem} script \texttt{path.FFADJ}, which is very similar to the script \texttt{path.FFELAS} described in \cref{sec.tutoelas}. Then, the function \texttt{gradS} of the {\texttt{mechtools.py}} module is in charge of calculating a shape gradient for $S(\Omega)$ at $\Omega^n$, in a analogous way to the computation described in the above \cref{sec.tutogcpv}. A few iterations of the resolution of \cref{eq.minS} in the setting of the L-shaped beam example are presented on \cref{fig.lb}. The source code associated to this development is contained in the following folder: \[ \texttt{\url{https://github.com/dapogny/sotuto/stress}} \] \begin{figure} \begin{center} \begin{overpic}[width=0.49\textwidth]{figures/lb0} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/lb26} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.49\textwidth]{figures/lb35} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/lb200} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{center} \caption{Iterations (a) $0$, (b) $26$, (c) $35$ and (d) $200$ (final) of the optimization of the shape of an L-beam with respect to an integral measure of the stress, as presented in \cref{man.stress}.}\label{fig.lb} \end{figure} \subsubsection{Using multiple equality or inequality constraints: optimization of the shape of a bridge submitted to multiple loads} \label{man.mload} In this section, we illustrate how shape optimization problems featuring multiple equality or inequality constraints can be handled with an adaptation of the base code. This task relies on the open-source implementation of the null space optimization algorithm of \cref{sec.nullspacealgo}, which is provided at the address \[ \texttt{\url{https://gitlab.com/florian.feppon/null-space-optimizer}} \] where a documentation of the library is available. The physical situation of interest is that depicted on \cref{fig.moreexamples}$\MK$(c): the considered shape $\Omega$ represents a two-dimensional bridge contained in a box $D$ with size $10 \times 2$. It is clamped on its lower-left and lower-right corners, and $7$ loads $g_i = (0,-1)$ are applied on as many different regions $\Gamma_{N,i}$ ($i=1,\ldots,7$) of its upper side. In this setting, we aim to solve the optimization problem \begin{equation}\label{eq.multicst} \min_\Omega \Vol(\Omega) \text{ s.t } C_i(\Omega) \leq C_T, \end{equation} In this formulation, $C_i(\Omega)$ is the compliance of $\Omega$ when the $i^{\text{th}}$ load $g_i$ is applied on $\Gamma_{N,i}$, involving the solution $u_{\Omega,i}$ to the linear elasticity system \cref{eq.elas} in this situation. The threshold $C_T$ is calculated from the values $C_i(\Omega^0)$ of the compliances of the initial design $\Omega^0$ in each load scenario: \[ C_T = 0.9 \:\max_{i=1,\ldots,7} C_i(\Omega^0). \] The treatment of this problem with the null space optimization library essentially requires a modification of the main loop of the optimization strategy. Briefly, the main steps of the optimization strategy are re-organized as the methods of a new class \texttt{Bridge}. \begin{lstlisting}[style=python,escapechar=|] ########### Definition of Optimizable class Bridge ############ class Bridge(Optimizable) : # Initialization def x0(self) : return path.step(0,"mesh") [...] # Calculation of the objective and constraint functions def evalObjective(self,x) : [...] for i in range(0,path.NC) : refneu = path.REFNEU + i gx = path.GX[i] gy = path.GY[i] ui = nam + ".u."+str(i)+".sol" mechtools.elasticity(x,refneu,gx,gy,ui) # Calculate J and H J = mechtools.volume(x) Cptab = [] for i in range(0,path.NC) : Cp = mechtools.compliance(x,nam+".u."+str(i)+".sol") [...] # Shape derivatives, sensitivity of objective and constraint def evalSensitivities(self,x) : [...] # Calculate dJ and gradJ mechtools.gradV(x,nam+".diffV.sol",nam+".gradV.sol") # Calculate dH and gradH for i in range(0,path.NC) : mechtools.gradCp(x,nam+".u."+str(i)+".sol", nam+".diffCp."+str(i)+".sol", nam+".gradCp."+str(i)+".sol") [...] # Retraction: shape update # dx = array of np values containing values # of the scalar velocity field def retract(self, x, dx) : [...] # Generation of a level set function for $\Omega^n$ on $D$ print(" Creation of a level set function") lstools.mshdist(curmesh,curphi) # Put scalar velocity defined on D in the normal direction inout.saveSol(dx,path.TMPSOL) lstools.norvec(curmesh,curphi,path.TMPSOL,curgrad) # Advection of the level set function print(" Level Set advection") lstools.advect(curmesh,curphi,curgrad,1.0,newphi) # Creation of a mesh associated to the new shape print(" Local remeshing") retmmg = mshtools.mmg2d(curmesh,1,newphi,path.HMIN, path.HMAX, path.HAUSD,path.HGRAD,1,newmesh) return newmesh # Accept step def accept(self,results) : [...] \end{lstlisting} The optimization is then launched by applying the \texttt{nlspace\_solve} function to an object of this class. \begin{lstlisting}[style=python,escapechar=|] # Run optimization solver optSettings = {"dt" : path.MESHSIZ, "alphaJ" : 1.0, "alphaC" : 2.0, "maxit" : 300, "provide_gradient" : True, "maxtrials" : 3, "itnormalisation" : 3 } results=nlspace_solve(Bridge(), optSettings) \end{lstlisting} A few intermediate shapes arising during the resolution of \cref{eq.multicst} are represented on \cref{fig.cst}; the source code corresponding to this example is contained in the following folder: \[ \texttt{\url{https://github.com/dapogny/sotuto/constraints}} \] \begin{figure}[!ht] \centering \begin{minipage}{0.9\textwidth} \includegraphics[width=1.0\textwidth]{./figures/cst0} \end{minipage} \begin{minipage}{0.9\textwidth} \includegraphics[width=1.0\textwidth]{./figures/cst15} \end{minipage} \begin{minipage}{0.9\textwidth} \includegraphics[width=1.0\textwidth]{./figures/cst200} \end{minipage}\caption{(From top to bottom) Iterations 0, 15 and 200 of the optimization of the shape of a 2d bridge submitted to 7 constraints on its compliance in as many load situations, as described in \cref{man.mload}.}\label{fig.cst} \end{figure} \subsubsection{A 3d implementation: optimization of a 3d cantilever beam} \label{man.3d} In this section, we briefly illustrate how our base implementation can be adapted to handle three-dimensional problems. The physical situation of interest is that of a 3d cantilever beam, which is sketched on \cref{fig.moreexamples}$\MK$(d): the shapes $\Omega$ are contained in a box $D$ with size $2 \times 1 \times 1$; they are clamped on their face lying in the plane $\left\{ x = (x_1,x_2,x_3) \in \partial D, \: x_1 = 0 \right\}$ and a unit vertical load $g = (0,0,-1)$ is applied on a small disk $\Gamma_N$ around the middle of the opposite side. In this context, we aim to minimize the compliance $C(\Omega)$ of $\Omega$ under a volume constraint: problem \cref{eq.sopbman} is solved with a volume target $V_T = 0.45$. The main numerical strategy is very similar to the two-dimensional base implementation. The three-dimensional version of \texttt{FreeFem} is used for the finite element calculation involved, which demands minor adaptations of the scripts described in the main part of this tutorial; see the documentation of \texttt{FreeFem} for further details. The external libraries \texttt{mshdist} and \texttt{advect} are called in exactly the same way as in the previous two-dimensional cases, and the three-dimensional version \texttt{mmg3d} of the remeshing software \texttt{mmg} is used in place of \texttt{mmg2d}. \begin{lstlisting}[style=python,escapechar=|] # Call for the executables of external codes FREEFEM = "FreeFem++ -nw" MSHDIST = "mshdist" ADVECT = "/Users/dapogny/Advection/build/advect" MMG3D = "/Users/dapogny/mmg/build/bin/mmg3d_O3" \end{lstlisting} \goodbreak{}\noindent The code \texttt{mmg3d} is called with the help of the \texttt{python} function \begin{lstlisting}[style=python,escapechar=|] def mmg3d(mesh,ls,phi,hmin,hmax,hausd,hgrad,nr,out) : \end{lstlisting} of the {\texttt{{mshtools.py}}} module, which operates in exactly the same way as the function \texttt{mmg2d} described in \cref{sec.manmmgmshdist}. A few iterations of the resolution of problem \cref{eq.sopbman} are represented on \cref{fig.canti3dman}; the associated source code is contained in the following folder: \[ \texttt{\url{https://github.com/dapogny/sotuto/base3d}} \] \begin{figure} \begin{center} \begin{overpic}[width=0.49\textwidth]{figures/canti3d0} \put(0,3){\fcolorbox{black}{white}{$a$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/canti3d50} \put(0,3){\fcolorbox{black}{white}{$b$}} \end{overpic} \end{center} \begin{center} \begin{overpic}[width=0.49\textwidth]{figures/canti3d120} \put(0,3){\fcolorbox{black}{white}{$c$}} \end{overpic} \begin{overpic}[width=0.49\textwidth]{figures/canti3d200} \put(0,3){\fcolorbox{black}{white}{$d$}} \end{overpic} \end{center} \caption{Iterations (a) $0$, (b) $50$, (c) $120$ and (d) $200$ (final) of the three-dimensional cantilever example presented in \cref{man.3d}.}\label{fig.canti3dman} \end{figure} \bibliographystyle{crplain} \bibliography{CRMATH_Dapogny_20221179} \end{document}