Fireshape: a shape optimization toolbox for Firedrake
FFireshape: a shape optimization toolbox for Firedrake
Alberto Paganini · Florian WechsungAbstract
We introduce Fireshape, an open-source andautomated shape optimization toolbox for the finite el-ement software Firedrake. Fireshape is based on themoving mesh method and allows users with minimalshape optimization knowledge to tackle with ease chal-lenging shape optimization problems constrained to par-tial differential equations (PDEs).
One of the ultimate goals of structural optimizationis the development of fully automated software thatallows users to tackle challenging structural optimiza-tion problems in the automotive, naval, and aerospaceindustries without requiring deep knowledge of struc-tural optimization theory. The scientific community isworking actively in this direction, and recent years haveseen the publication of educational material that sim-plifies the understanding of structural optimization al-gorithms and guides the development of related opti-mization software. These resources are based on differ-ent models, such as moving mesh methods [3, 5, 16, 18],level-sets [4,41], phase fields [29], and SIMP [11,53,54],and are implemented in various software environmentssuch as Matlab [53, 54], FreeFem++ [4, 5, 18], Open- A. PaganiniSchool of Mathematics and Actuarial Science, University ofLeicester, University Road, Leicester, LE1 7RH, United King-dom,E-mail: [email protected]. WechsungCourant Institute of Mathematical Sciences, New York Uni-versity, 251 Mercer St, New York, NY 10012,E-mail: [email protected] The acronym SIMP stands for Solid Isotropic Materialwith Penalisation.
FOAM CFD [16], FEMLAB [43], and FEniCS [21, 41],to mention just a few.In this work, we introduce Fireshape: an automatedshape optimization library based on the moving meshesapproach that requires very limited input from the user.Shape optimization refers to the optimization of domainboundaries and plays an important role in structuraldesign. For instance, shape optimization plays a cru-cial role in the design of airfoils [33, 51, 55] and boathulls [1, 44, 46]. Shape optimization is also a useful re-finement step to be employed after topology optimiza-tion [11, Ch. 1.4]. Indeed, topology optimization al-lows more flexibility in geometric changes and it is apowerful tool to explore a large design space. However,topology optimization usually provides slightly blurred(grey-scale) and/or staircase designs [11]. By addinga final shape optimization step, it is possible to post-process results computed with topology optimizationand devise optimal designs with sharp boundaries andinterfaces.Fireshape is based on the moving meshes shape op-timization approach [3, Ch. 6]. In this approach, ge-ometries are parametrized with meshes that can bearbitrarily precise and possibly curvilinear. The meshnodes and faces are then optimized (or “moved”) tominimize a chosen target function. The moving meshesapproach is only one among several possible shape mod-els, such as phase-fieds [13, 23, 29], level-sets [4, 12, 41],and SIMP [10, 11, 53, 54], and each model comes withadvantages and disadvantages. Fireshape has been de-veloped on the moving meshes approach because thelatter has a very neat interpreation in terms of geomet-ric transformations and is inherently compatible withstandard finite element software [49]. The main draw-back of the moving meshes approach is that it does notallow topological changes in a straightforward and con-sistent fashion. However, Fireshape has been developed a r X i v : . [ m a t h . O C ] M a y Alberto Paganini, Florian Wechsung to facilitate shape optimization, and topology optimiza-tion is beyond its scope.Fireshape couples the finite element library Fire-drake [31, 36, 37, 45, 50] with the Rapid OptimizationLibrary (ROL) [20]. Fireshape allows decoupled dis-cretization of geometries and state constraints and itincludes all necessary routines to perform shape opti-mization (geometry updates, regularization terms, geo-metric constraints, etc.). To solve a shape optimizationproblem in Fireshape, users must describe the objectivefunction and the eventual constraints using the UnifiedForm Language (UFL) [6], a Python-based scriptinglanguage that is very similar to standard mathemati-cal notation. Once objective functions and constraintshave been implemented with UFL, users need only toprovide a mesh that describes the initial design and, fi-nally, select their favorite optimization algorithm fromthe optimization library ROL.Typically the bottleneck of PDE constrained opti-mization code lies in the solution of the state and theadjoint equation. While Fireshape and Firedrake areboth written in Python, to assemble the state and ad-joint equations the Firedrake library automatically gen-erates optimized kernels in the programming languageC. The generated systems of equations are then passedto the PETSc library (also in written in C) and canbe solved using any of the many linear solver and pre-conditioners provided by PETSc [7–9, 14, 17, 42]. Thiscombination of Python for user facing code and C forperformance critical parts is well established in scien-tific computing as it provides highly performant codethat is straightforward to develop and use. Finally, wemention that Fireshape, just as Firedrake, PETSc andROL, supports MPI parallelization and hence can beused to solve even large scale three dimensional shapeoptimization problems.Fireshape is an open-source software licensed underthe GNU Lesser General Public License v3.0. Fireshapecan be downloaded at https://github.com/Fireshape/Fireshape . Its documentation contains several tutori-als and is available at https://fireshape.readthedocs.io/en/latest/ . To illustrate Fireshape’s capabilitiesand ease of use, in Section 2 we provide a tutorial andsolve a three-dimensional shape optimization problemconstrained to a nonlinear boundary value problem.The shape optimization knowledge required to under-stand this tutorial is minimal. In Section 3, we describein detail the rich mathematical theory that underpinsFireshape. In Section 4, we describe Fireshape’s mainclasses and Fireshape’s extended functionality. Finally,in Section 5, we provide concluding remarks.
For this hands-on introduction to Fireshape, we con-sider a viscous fluid flowing through a pipe Ω (see Fig-ure 1) and aim to minimize the kinetic energy dissipa-tion into heat by optimizing the design of Ω. To beginwith, we need to describe this optimization problemwith a mathematical model.
Fig. 1
Viscous fluid flows through a pipe Ω from the left tothe right. A poor pipe design can lead to an excessive amountof kinetic energy being is dissipated into heat.
We assume that the fluid velocity u and the fluidpressure p satisfy the incompressible Naver-Stokes equa-tions [25, Eqn. 8.1], which read − ν ∇ · ε u + ( u · ∇ ) u + ∇ p = 0 in Ω , (1a)div u = 0 in Ω , (1b) u = g on ∂ Ω \ Γ , (1c)p n − νε u n = 0 on Γ , (1d)where ν denotes the fluid viscosity, ε u = ε ( u ) = ( ∇ u + ∇ u (cid:62) ), ∇ u is the derivative (Jacobian matrix) of u and ∇ u (cid:62) is the derivative transposed, Γ denotes the outlet,and the function g is a Poiseuille flow velocity [25, p.122] at the inlet and is zero on the pipe walls. In ournumerical experiment, the inlet is a disc of radius 0.5in the xy-plane, and g ( x, y, z ) = (0 , , − x + y )) (cid:62) on the inlet.To model the kinetic energy dissipation, we considerthe diffusion term in (1a) and introduce the functionJ(Ω) = (cid:90) Ω ν ε u : ε u d x , (2)where the colon symbol : denotes the Frobenius innerproduct, that is, ε u : ε u = trace(( ε u ) (cid:62) ε u ).Now, we can formulate the shape optimizations prob-lem we are considering. This reads:Find Ω ∗ such that J(Ω ∗ ) = inf J(Ω) subject to (1). (3)To make this test case more interesting, we further im-pose a volume equality constraint on Ω. Otherwise, thesolution to this problem would be a pipe with arbitrar-ily large diameter.In the next subsections, we explain step by step howto solve this shape optimization problem in Fireshape. ireshape: a shape optimization toolbox for Firedrake 3 To this end, we need to create: a mesh that approxi-mates the initial guess Ω (Section 2.1), an object thatdescribes the PDE-constraint (1) (Section 2.2), an ob-ject that describes the objective function (2) (Section2.3), and, finally, a “main file” to set up and solve theoptimization problem (3) (Section 2.4).2.1 Step 1: Provide an initial guessThe first step is to provide a mesh that describes an ini-tial guess of Ω. For this tutorial, we create this mesh us-ing the software Gmsh [30]. The initial guess employedis sketched in Figure 1. For the geometric details, werefer to the code archived on Zenodo [58].2.2 Step 2: Implement the PDE-constraintThe second step is to implement a finite element solverof the Navier-Stokes equations (1). To derive the weakformulation of (1a), we multiply Equation (1a) with atest (velocity) function v that vanishes on ∂ Ω \ Γ, andEquation (1b) with a test (pressure) function q. Then,we integrate over Ω, integrate by parts [25, Eq. 3.18],and impose the boundary condition (1d). The resultingweak formulation of equations (1) reads:Find u such that u = g on ∂ Ω \ Γ and (cid:90) Ω νε u : ε v − p div( v ) + v · ( ∇ u ) u + q div( u ) d x = 0for any pair ( v , q) . (4)To implement the finite element discretization of(4), we create a class NavierStokesSolver that in-herits from the Fireshape’s class
PdeConstraint , seeListing 1. To discretize (4), we employ P2-P1 Taylor-Hood finite elements, that is, we discretize the trial andtest velocity functions u and v with piecewise quadraticLagrangian finite elements, and the trial and test pres-sure functions p and q with piecewise affine Lagrangianfinite elements. It is well known that this is a stablediscretization of the incompressible Navier-Stokes equa-tions [25, pp 136-137].To address the nonlinearity in (4), we use PETSc’sScalable Nonlinear Equations Solvers (SNES) [7–9, 14,17, 42]. In each iteration, SNES linearizes Equation (4)and solves the resulting system with a direct solver.In general, it is possible that the SNES solver failsto converge sufficiently quickly (in our code, we allowat most 100 SENS iterations). For instance, this hap-pens if the finite element mesh self-intersects, or if theinitial guess used to solve (4) is not sufficiently good. from firedrake import * from fireshape import PdeConstraint class NavierStokesSolver(PdeConstraint): """Incompressible Navier-Stokes equations.""" def __init__(self, mesh_m, viscosity): super().__init__() self.mesh_m = mesh_m self.failed_to_solve = False self.V = VectorFunctionSpace(mesh_m, "CG", 2) \ * FunctionSpace(mesh_m, "CG", 1) self.solution = Function(self.V, name="State") self.testfunction = TestFunction(self.V) self.viscosity = viscosity u, p = split(self.solution) v, q = split(self.testfunction) nu = viscosity self.F = nu*2*inner(sym(grad(u)), sym(grad(v)))*dx \ - p*div(v)*dx + div(u)*q*dx \ + inner(dot(grad(u), u), v)*dx X = SpatialCoordinate(mesh_m) rsq = X[0]**2+X[1]**2 uin = as_vector([0, 0, 1-4*rsq]) bc1 = DirichletBC(self.V.sub(0), 0., [12, 13]) bc2 = DirichletBC(self.V.sub(0), uin, [10]) self.bcs = [bc1, bc2] self.nsp = None self.params = {"snes_max_it": 100, "mat_type": "aij", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"} def solve(self): super().solve() self.failed_to_solve = False u_old = self.solution.copy(deepcopy=True) try: solve(self.F == 0, self.solution, bcs=self.bcs, solver_parameters=self.params) except ConvergenceError: self.failed_to_solve = True self.solution.assign(u_old) Lst. 1
Incompressible Navier-Stokes’s solver based onTaylor-Hood finite elements.
Most often, these situations happen when the optimiza-tion algorithm takes an optimization step that is toolarge. In these cases, we can address SNES’ failure tosolve (4) by reducing the optimization step. In prac-tice, we deal with these situations with Python’s try:... except: ... block. If the SNES solver fails witha
ConvergenceError , we catch this error and set theboolean flag self.failed to solve to True . This flag
Alberto Paganini, Florian Wechsung is used to adjust the output of the objective function J,see Section 2.3.2.3 Step 3: Implement the objective functionThe third step is to implement a code to evaluate theobjective function J defined in Equation (2). For this,we create a class
EnergyDissipation that inherits fromFireshape’s class
ShapeObjective , see Listing 2. Ofcourse, evaluating J requires access to the fluid velocity u . This access is implemented by assigning the variable self.pde solver . This variable gives access also to thevariable NavierStokesSolver.failed to solve , whichcan be used to control the output of J when the Navier-Stokes’ solver fails to converge. Here, we decide that thevalue of J is
NaN (“not a number”) if the Navier-Stokes’solver fails to converge. from firedrake import * from fireshape import ShapeObjective from PDEconstraint import NavierStokesSolver import numpy as np class EnergyDissipation(ShapeObjective): """Kinetic energy dissipation.""" def __init__(self, pde_solver: NavierStokesSolver, *args, **kwargs): super().__init__(*args, **kwargs) self.pde_solver = pde_solver def value_form(self): """Evaluate misfit functional.""" nu = self.pde_solver.viscosity if self.pde_solver.failed_to_solve: return np.nan * dx(self.pde_solver.mesh_m) else: z = self.pde_solver.solution u, p = split(z) return nu * inner(sym(grad(u)), sym(grad(u))) * dx Lst. 2
Objective function J (2), which quantifies the kineticenergy dissipation of the fluid. Note that J returns
NaN (“nota number”) if the solver of the PDE-constraint (4) fails toconverge. from firedrake import * from fireshape import * import fireshape.zoo as fsz import ROL from PDEconstraint import NavierStokesSolver from objective import EnergyDissipation mesh = Mesh("pipe.msh") bbox = [(-0.5, +0.5), (-0.5, +5.5), (2, 12)] orders = [3, 3, 3]; levels = [2, 2, 4] Q = BsplineControlSpace(mesh, bbox, orders, levels, boundary_regularities=[0, 0, 2], fixed_dims=[0, 2]) q = ControlVector(Q, LaplaceInnerProduct(Q)) viscosity = Constant(1/10.) e = NavierStokesSolver(Q.mesh_m, viscosity) out = File("solution/u.pvd") def cb(): return out.write(e.solution.split()[0]) J_ = EnergyDissipation(e, Q, cb=cb) J = ReducedObjective(J_, e) Jq = fsz.MoYoSpectralConstraint(10, Constant(0.5), Q) J = J + Jq vol = fsz.VolumeFunctional(Q) vol0 = vol.value(q, None) C = EqualityConstraint([vol], target_value=[vol0]) M = ROL.StdVector(1) params_dict = { 'General': {'Secant': {'Type': 'Limited-Memory BFGS', 'Maximum Storage': 20}}, 'Step': {'Type': 'Augmented Lagrangian', 'Augmented Lagrangian': {'Subproblem Step Type': 'Trust Region', 'Subproblem Iteration Limit': 50, 'Penalty Parameter Growth Factor': 2}}, 'Status Test': {'Gradient Tolerance': 1e-2, 'Step Tolerance': 1e-2, 'Constraint Tolerance': 1e-2, 'Iteration Limit': 10}} params = ROL.ParameterList(params_dict, "Parameters") problem = ROL.OptimizationProblem(J, q, econ=C, emul=M) solver = ROL.OptimizationSolver(problem, params) solver.solve() Lst. 3
Minimize the kinetic energy dissipation of a fluid us-ing Fireshape. • In lines 1-6, we import all necessary Python librariesand modules. • In lines 8-13, we load the mesh of the initial designand we specify how to discretize the domain Ω (seeSections 4.1 and 4.2 for more details about the Fire-shape classes
ControlSpace and
InnerProduct ). • In lines 15-16, we initiate the Navier-Stokes’ finiteelement solver. ireshape: a shape optimization toolbox for Firedrake 5 • In lines 18-19, we tell Fireshape to store the finiteelement solution to Navier-Stokes’ equations in thefolder solution everytime the domain Ω is updated.We can visualize the evolution of u along the op-timization process by opening the file u.pvd withParaview [2]. • In lines 21-22, we initiate the objective function Jand the associated reduced functional, which is usedby Fireshape to define the appropriate Lagrangefunctionals to shape differentiate J. We refer to Sec-tion 3.2 for more details about shape differentiationusing Lagrangians. We stress that Fireshape doesnot require users to derive (and implement) shapederivatives’ formulas by hand. The whole shape dif-ferentiation process is automated using pyadjoint[22, 32]. • In lines 24-25, we add an additional regularizationterm to J to promote mesh quality of domain up-dates. This regularization term controls the point-wise maximum singular value of the geometric trans-formation used to update Ω [56]. See Section 3.2 andFigure 5 for more information about the role of ge-ometric transformations in Fireshape. • In lines 27-30, we set up an equality constraint toensure that the volume of the initial and the opti-mized domains are equal. • In lines 32-43, we select an optimization algorithmfrom the optimization library ROL. More specifi-cally, we use an augmented Lagrangian approach[48, Ch. 17.3] with limited-memory BFGS Hessianupdates [48, Ch. 6.1] to deal with the volume equal-ity constraint, and a Trust-Region algorithm [48,Ch. 4.1] to solve the intermediate models generatedby the augmented Lagrangian algorithm. • In lines 44-47, finally, we gather all information andsolve the problem.2.5 ResultsRunning the code contained in Listing 3 optimizes thepipe design when the fluid viscosity ν is 0.1, whichcorresponds to Reynold-number Re = 1 ./ν = 10. Ofcourse, the resulting shape depends on the fluid vis-cosity. A natural question is how the optimized shapedepends on this parameter. To answer this question, wecan simply run Listing 3 for different Reynold-numbers(by modifying line 13 with the desired value) and in-spect the results. Here we perform this comparison forRe ∈ { , , , , } . For this computation we extend the code in Listing 1 andinclude analytic continuation in the Reynolds-number to re-compute good initial guesses for the SNES solver when this
Fig. 2
Initial design of the pipe and velocity solutions forRe = 1 , , , , Fig. 3
Optimized shapes for Re = 1 , , , , In Figure 2, we show the initial design and plot themagnitude of the fluid velocity on a cross section ofthe pipe for different Reynold-numbers. In Figure 3,we show the resulting optimized shapes, and in Fig-ure 4 the corresponding magnitudes of the fluid velocity.Qualitatively, we observe that, as the Reynolds-numberincrease, we obtain increasingly S-shaped designs thatavoid high curvature at the two fixed ends. Finally, weremark that the objective is reduced by approximately6 . . . .
7% and 12 . In this section, we describe the theory that underpinsFireshape. We begin with a brief introduction to PDE- fails to converge. This modification is necessary to solve (4)at high Reynolds-numbers. The code used to obtain theseresults can be found at [58]. Alberto Paganini, Florian Wechsung
Fig. 4
Velocity field for optimized shapes for Re =1 , , , , constrained optimization to set the notation and men-tion the main idea behind optimization algorithms (Sub-section 3.1). Then, we continue with an introduction toshape calculus (Subsection 3.2). Finally, we concludewith a discussion about the link between shape op-timization and parametric finite elements (Subsection3.3).3.1 Optimization with PDE-constraintsThe basic ingredients to formulate a PDE-constrainedoptimization problems are: a control variable q thatlives in a control space Q , a state variable u that lives ina state space U and that solves a (possibly nonlinear)state equation A ( q , u ) = , and a real function J : Q × U → R to be minimized. Example:
Equation (3) is a PDE-constrained op-timization problem. The control variable q correspondsto the domain Ω , the state variable u represents thepair velocity-pressure ( u , p) , the nonlinear constraint A represents the Navier-Stokes equations (1) , and J denotes the objective function (2) . The state space U corresponds to the space of pairs ( a , b) with weakly dif-ferentiable velocities a that satisfy a = g on ∂ Ω \ Γ andsquare integrable pressures b [25, Ch. 8.2]. The controlspace Q is specified in Section 3.2; see Equation (5) . Most numerical methods for PDE-constrained op-timization attempt to construct a sequence of controls { q ( k ) } k ∈ N and corresponding states { u ( k ) } k ∈ N such thatlim k →∞ J( q ( k ) , u ( k ) ) = inf q , u { J( q , u ) : A ( q , u ) = } . Often, the sequence { q ( k ) } k ∈ N is constructed using deriva-tives of the function J and of the constraint A . Common An alternative approach is to employ so-called “one-shot”methods, where the optimality system of the problem issolved directly [52]. Note that, since the optimality systemis often nonlinear, one-shot methods still involve iterative al-gorithms. approaches are steepest descent algorithms and New-ton methods (in their quasi, Krylov, or semi-smoothversions) [34, 40]. These algorithms ensure that the se-quence { J( q ( k ) , u ( k ) ) } k ∈ N decreases. In special cases, itis even possible to show that the sequence { q ( k ) } k ∈ N converges [34]. Although it may be difficult to ensurethese assumptions are met in industrial applications,these optimization algorithms are still powerful toolsto improve state-of-the-art designs and are widely usedto perform shape optimization.3.2 Shape optimization and shape calculusShape optimization with PDE constraints is a particu-lar branch of PDE-constrained optimization where thecontrol space Q is a set of domain shapes. The space ofshapes is notoriously difficult to characterize uniquely.For instance, one could describe shapes through theirboundary regularity, or as level-sets, or as local epigraphs[19, ch. 2]. The choice of the shapes’ space character-ization plays an important role in the concrete imple-mentation of a shape optimization algorithm, and itcan also affect formulas that result by differentiatingJ with respect to perturbations of the shape q . To bemore precise, different methods generally lead to thesame first order shape derivative formula, but differ onshape derivatives of higher order [19, ch. 9]In this work, we model the control space Q as theimages of bi-Lipschitz geometric transformations T ap-plied to an initial set Ω ⊂ R d [19, ch. 3], that is, Q := { q = T (Ω) : T : R d → R d is bi-Lipschitz } , (5)see Figure 5. We choose this model of Q because it pro-vides an explicit description of the domain boundariesvia ∂ ( T (Ω)) = T ( ∂ Ω), and because it is compatiblewith higher-order finite elements [49], as explained indetail in the Section 3.3. Henceforth, we use the termdifferomorphism to indicate that a geometric transfor-mation T is bi-Lipschitz. Ω ...
Fig. 5
To construct the control space Q , we select an initialset Ω and a set of diffeomorphisms { T } ( black arrows ), andcollect every image T (Ω).ireshape: a shape optimization toolbox for Firedrake 7 In this setting, the shape derivative of J correspondsto the classical Gˆateaux derivative in the Sobolev space W , ∞ ( R d , R d ) [34, Def. 1.29]. To see this, let us momen-tarily remove the PDE-constraint, and only considera shape functional of the form Ω (cid:55)→ J(Ω). The shapederivative of J at Ω is the linear and continuous opera-tor dJ(Ω , · ) : W , ∞ ( R d , R d ) → R defined bydJ(Ω , V ) := lim t (cid:38) J(Ω t ) − J(Ω)t , (6)where Ω t is the set Ω t = { x + t V ( x ) : x ∈ Ω } [3, Def.6.15]. By replacing Ω with q = I (Ω), where I denotesthe identity transformation defined by I ( x ) = x for any x ∈ R d , Equation (6) can be equivalently rewritten asdJ( q , V ) := lim t (cid:38) J( q + t V ) − J( q )t . We highlight that this interpretation immediately gen-eralizes to any q in Q and can be used to define higherorder shape derivatives.The same definition (6) of shape derivative holds inthe presence of PDE-constraints: the shape derivativeat q of the function J : Q × U → R subject to A ( q , u ) = is the linear and continuous operator defined bydJ( q , u , V ) := lim t (cid:38) J( q + t V , u t ) − J( q , u )t , (7)where u t is the solution to A ( q + t V , u t ) = . Com-puting shape derivative formulas using (7) may presentthe difficulty of computing the shape derivative of u ,which intuitively arises by the “chain rule” formula.The shape derivative of u , which is sometimes called“material derivative” of u , can be eliminated with thehelp of adjoint equations [34, Sec. 1.6.2]. This processcan be automated by introducing the Lagrange func-tional [34, Sec. 1.6.3] L ( q , u , p ) := J( q , u ) + (cid:104) A ( q , u ) , p (cid:105) . (8)The term (cid:104) A ( q , u ) , p (cid:105) stems from testing the equation A ( q , u ) = with a test function p in the same way itis usually done when writing a PDE in its weak form. Example: If A denotes the Navier-Stokes equa-tions (1) , then (cid:104) A ( q , u ) , p (cid:105) corresponds to the weak for-mulation (4) . The advantage of introducing the Lagrangian (8)is that, by choosing p as the solution to the adjointequation (cid:104) ∂ u A ( q , u ) , p (cid:105) = − ∂ u J( q , u ) for all u in U , where ∂ u denotes the partial differentiation with re-spect to the variable u , the shape derivative of J canbe computed asdJ( q , u , V ) = ∂ q L ( q , u , p , V ) , where ∂ q denotes partial differentiation with respectto the variable q . This is advantageous because partialdifferentiation does not require computing the shapederivative of u .The Lagrangian approach to compute derivativesof PDE-constrained functionals is well established [39],and its steps can be replicated by symbolic compu-tation software. Probably, the biggest success in thisdirection is the dolfin-adjoint project [26, 47], whichderives “the adjoint and tangent-linear equations andsolves them using the existing solver methods in FEn-iCS/Firedrake” [47]. Thanks to the shape differentia-tion capabilities of UFL introduced in [32], dolfin-adjointis also capable of shape differentiating PDE-constrainedfunctionals [22]. Remark 1
Optimization algorithms are usually basedon steepest descent directions to update the controlvariable q . A steepest descent direction is a direction V ∗ that minimizes the derivative dJ. Since dJ is linear,it is necessary to introduce a norm (cid:107) · (cid:107) on the spaceof possible perturbations { V } (the tangent space of Q at q ), and to restrict the search of a steepest descentdirection to directions V of length (cid:107) V (cid:107) = 1. A naturalchoice would be to select the W , ∞ -norm , with respectto which a minimizer has been shown to exists for mostfunctionals [49, Prop. 3.1]. However, in practice it ismore convenient to employ a norm that is induced byan inner product, so that the steepest descent direc-tion corresponds to the Riesz representative of dJ withrespect to the inner product [34, p. 98].3.3 Geometric transformations, moving meshes, andparametric finite elementsTo solve a PDE-constrained optimization problem iter-atively, it is necessary to employ a numerical methodthat is capable of solving the constraint A ( q , u ) = forany feasible control q . In shape optimization, this trans-lates into the requirement of a numerical scheme thatcan solve a PDE on a domain that changes at each iter-ation of the optimization algorithm. There are severalcompeting approaches to construct a numerical schemewith this feature, such as the level-set method [4] or thephase field approach [13], among others.In Fireshape, we employ the approach sometimesknow as “moving mesh” method [3, 49]. In its simplestversion (see Figure 6), this method replaces (or approx-imates) the initial domain Ω with a polygonal mesh Ω h .On this mesh, the state and adjoint equations are solved The norm (cid:107) V (cid:107) , ∞ is the maximum between the essentialsupremum of V and of its derivative DV . Alberto Paganini, Florian Wechsung with a standard finite element method, whose construc-tion on polygonal meshes is immediate. For instance,depending on the nature of the state constraint, onemay solve the state and adjoint equations using linearLagrangian or P2-P1 Taylor-Hood finite elements (seeSection 2.2). With the state and adjoint solutions athand, one employs shape derivatives (see Remark 1 andSection 4.2) to update the coordinates of mesh nodeswhile retaining the mesh connectivity of Ω h . This leadsto a new the mesh that represents an improved design.This update process is repeated until some prescribedconvergence criteria are met. Fig. 6
Simplest approach to PDE-constrained shape opti-mization. An initial mesh ( left ) is used to compute stateand adjoint variables. This information, together with shapederivatives, is used to devise an update of the nodes’ coordi-nates ( center ). A new and improved initial guess is obtainedby updating the nodes’ coordinates and retaining the initialmesh connectivity ( right ). This process is repeated until con-vergence.
The moving mesh method we just described is a sim-ple and yet powerful method. However, in its currentformulation, it requires polyhedral meshes, which lim-its the search space Q to polyhedra. In the remainingpart of this section, we describe an equivalent interpre-tation of the moving mesh method that generalizes tocurved domains. Additionally, this alternative interpre-tation allows approximating state and adjoint variableswith arbitrarily high-order finite elements without suf-fering from reduction of convergence order due to poorapproximation of domain boundaries [15, Ch. 4.4].We begin by recalling the standard construction ofparametric finite elements. For more details, we referto [15, Sect. 2.3 and 4.3]. To construct a parametricfinite element space, one begins by partitioning the do-main Ω into simpler geometric elements { K i } (usuallytriangles or tetrahedra, as depicted in the first row ofFigure 7). Then, one introduces a reference elementˆ K and a collection of diffeomorphisms F i that mapthe reference element ˆ K to the various K i s, that is, F i ( ˆ K ) = K i for every value of the index i (as depictedin the second row of Figure 7). Finally, one considers aset of local reference basis functions { ˆ b j } defined on thereference element ˆ K , and defines local basis functions { b ij } on each K i via the pullback b ij ( x ) := ˆ b j ( F − i ( x ))for every x in K i . These local basis functions are used to construct global basis functions that span the finiteelement space. An important property of this construc-tion is that the diffeomorphisms F i can be expressed interms of local linear Lagrangian basis functions { ˆ β m } ,that is, there are some coefficient vectors { µ im } suchthat F i ( x ) = (cid:88) m µ im ˆ β m ( x ) for every x in ˆ K . (9)Finite elements constructed following this procedureare usually called parametric, because they rely on theparametrization { F i } . Note that the most common fi-nite elements families, such as Lagrangian, Raviart–Thomas, or Nedelec finite elements, are indeed para-metric. Ω K K ...K F i ˆ K K i Fig. 7
To construct finite element basis functions, one usu-ally begins by triangulating a domain Omega ( first row ) andthen introducing diffeomorphisms F i that map the referenceelement ˆ K to the various K i in the triangulation ( secondrow ). Keeping this knowledge about parametric finite ele-ments in mind, we can revisit the moving mesh method(see Figure 6). There, the main idea was to updateonly nodes’ coordinates and keep the mesh connectivityunchanged, so that constructing finite elements on thenew mesh is straightforward. In [49], it has been shownthat the new finite element space can also be obtainedby modifying the parametric construction of finite ele-ments on the original domain Ω. In the next paragraphs,we give an extended explanation (with adapted nota-tion) of the demonstration given in [49].Let T denote the transformation employed to mod-ify the mesh Ω h on the left in Figure 6 into the newand perturbed mesh T (Ω h ) on the right in Figure 6.Additionally, let { T ( K i ) } denote the simple geomet-ric elements that constitute the latter. Using the para-metric approach, we can construct finite elements onthe new mesh by introducing a collection of diffeomor-phisms ˜ F i that map the reference element ˆ K to the For linear Lagrangian finite elements, the two set of ref-erence local basis functions { ˆ b j } and { ˆ β m } coincide.ireshape: a shape optimization toolbox for Firedrake 9 various T ( K i )s, that is, ˜ F i ( ˆ K ) = T ( K i ) for every valueof the index i ; see Figure 8. ˜ F i F i T ˆ K K i T ( K i ) Fig. 8
The standard construction of finite elements on theperturbed mesh { T ( K i ) } employs diffeomorphisms ˜ F i thatmap the reference element ˆ K to the various perturbed ele-ments T ( K i ). The behavior of the transformation T in Figure 6 isprescribed only the mesh nodes. Since its behavior onthe interior of the mesh triangles can be chosen arbi-trarily, we can decide that T is piecewise affine on eachtriangle. This convenient choice implies that T can bewritten as a linear combination of piecewise affine La-grangian finite elements defined on the first mesh, thatis, T ( x ) = (cid:80) (cid:96) ν (cid:96) B (cid:96) ( x ) for every x in Ω, where { ν (cid:96) } are some coefficient vectors and { B (cid:96) } are global basisfunctions of the space of piecewise affine Lagrangianfinite elements defined on the partition { K i } . Since La-grangian finite elements are constructed via pullbacksto the reference element, for every element K i there arecoefficient vectors { ν im } so that the restriction T | K i of T on K i can be rewritten as T | K i ( x ) = (cid:88) (cid:96) ν (cid:96) B (cid:96) | K i ( x ) = (cid:88) m ˜ ν im ˆ β m ( F − i ( x )) . Therefore, the composition T | K i ◦ F i is of the form T | K i ◦ F i ( x ) = (cid:88) m ˜ ν im ˆ β m ( x ) for every x ∈ ˆ K , (10)that is, of the same form of Equation (9). This impliesthat, to construct finite elements on the perturbed ge-ometry T (Ω h ), we only need to replace the original coef-ficients { µ im } in Equation (9) with the new coefficients { ˜ ν im } from Equation (10).This alternative and equivalent viewpoint on themoving mesh method generalizes naturally to higher-order finite element approximations. Indeed, one of thekey steps to ensure that higher-order finite elementsachieve higher-order convergence on curved domains isto employ sufficiently accurate polynomial interpola-tion of domain boundaries [15, Ch. 4.4]. This boundaryinterpolation can be encoded in the diffeomorphisms { F i } by using higher-order Lagrangian local basis func-tions. Therefore, simply employing higher-order Lagrangianfinite element transformations T leads to a natural ex-tension of moving mesh method to higher-order finiteelements.This alternative and equivalent viewpoint general-izes further to allow the use of any arbitrary discretiza-tion of the transformation T (for instance, using B-splines [35], harmonic polynomials, or radial basis func-tions [57]). The only requirement to be fulfilled to en-sure the desired order of convergence p is that the maps T ◦ F i : x (cid:55)→ T ( F i ( x )) satisfy the asymptotic algebraicestimates (cid:107) D α ( T ◦ F i ) (cid:107) = O ( h α ) for 0 ≤ α ≤ p , where D α ( T ◦ F i ) denotes the α th derivative of T ◦ F i .Using a different discretization of the transformation T can give several advantages, like increasing the smooth-ness T (because finite elements are generally only Lip-schitz continuous) or varying how shape updates arecomputed during the optimization process [24]. Remark 2
An issue that can arise with the moving meshmethod is that it can lead to poor quality (or even tan-gled) meshes. In terms of geometric transformations, amesh with poor quality corresponds to a transforma-tion T for which the value max α (cid:107) D α T (cid:107) is large (and atangled mesh to a transformation T that is not a diffeo-morphism). To a certain extent, it is possible to enforcemoderate derivatives by employing suitable metrics toextract descent directions from shape derivatives (forinstance, by using linear elasticity based inner productswith a Cauchy-Riemann augmentation [38]) and/or byadding penalty terms to the functional J (as in Section2.4). In this section, we give more details about Fireshape’simplementation and features. Fireshape is organized ina few core Python classes (and associated subclasses)that implement the control space Q , the metric to beemployed by the optimization algorithm, and the (pos-sibly PDE-constrained) objective function J. The fol-lowing subsections describe these classes.4.1 The class ControlSpace
Fireshape models the control space Q of admissible do-mains using geometric transformations T as in Equa-tion (5) (see also Figure 5). From a theoretical perspec-tive, the transformations T can be discretized in nu-merous different ways as long as certain minimal reg-ularity requirements are met. In Fireshape, the class ControlSpace allows the following options: (i) Lagrangianfinite elements defined on the same mesh employed tosolve the state equation, (ii) Lagrangian finite elementsdefined on a mesh coarser than the mesh employed tosolve the state equation, and (iii) tensorized B-splinesdefined on a separate Cartesian grid (not to be confusedwith a spline or B´ezier parametrization of the domainboundary, see Figure 9). Ω Fig. 9
Fireshape allows discretizing transformations T usingtensorized B-splines. On the left, we display a Cartesian gridthat that covers the computational domain Ω. On the right,we display a quadratic tensorized B-spline. If discretization (i) can be considered to be the de-fault option, discretization (ii) allows introducing a reg-ularization by discretizing the geometry more coarsely(so-called “regularization by discretization”), whereasB-splines allow constructing transformation with higherregularity (Lagrangian finite elements are only Lips-chitz continuous, whereas B-splines can be continuouslydifferentiable and more).The class
ControlSpace can be easily extended toinclude additional discretization options, such as radialbasis functions and harmonic polynomials.4.2 The class
InnerProduct
To formulate a continuous optimization algorithm, weneed to specify how to compute lengths of vectors (and,possibly, how to compute the angle between two vec-tors). The class
InnerProduct addresses this require-ment and allows selecting an inner product ( · , · ) H toendow the control space Q with .The choice of the inner product affects how steepest-descent directions are computed. Indeed, a steepest-descent direction is a direction V of length (cid:107) V (cid:107) H = 1such that dJ( q , V ) is minimal. Let α := (cid:107) dJ( q , · ) (cid:107) ∗ denote the length of the operator (cid:107) dJ( q , · ) (cid:107) measuredwith respect to the dual norm. Then [34, p. 103], thesteepest descent direction V satisfies the equation α ( V , W ) H = − dJ( q , W ) for all W in H , Note that specifying a norm would suffice to formulate anoptimization algorithm. However, as mention in Remark 1, itis computationally more convenient to restrict computationsto inner product spaces. which clearly depends on the inner product ( · , · ) H .The control space Q can be endowed with differentinner products. In Fireshape, the class InnerProduct allows the following options: (i) an H (Ω) inner productbased on standard Galerkin stiffness and mass matri-ces, (ii) a Laplace inner product based on the Galerkinstiffness matrix, and (iii) an elasticity inner productbased on the linear elasticity mechanical model. Theseoptions can be complemented with additional homoge-neous Dirichlet boundary conditions to specify parts ofthe boundary ∂ Ω that are not to be modified duringthe shape optimization procedure.Although all three options are equivalent from a the-oretical perspective, in practice it has been observedthat option (iii) generally leads to geometry updatesthat result in meshes of higher equality compared tooptions (i) and (ii). A thorough comparison is avail-able in [38], where the authors also suggest to considercomplementing these inner products with terms stem-ming from Cauchy-Riemann equations to further in-crease mesh quality. This additional option is readilyavailable in Fireshape.4.3 The classes
Objective and
PdeConstraint
In the vast majority of cases, users who aim to solve aPDE-constrained shape optimization problem are onlyrequired to instantiate the two classes
Objective and
PdeConstraint , where they can specify the formula ofthe function J to be minimized and the weak formula-tion of its PDE-constraint A (Ω , u ) = (see Sections2.2 and 2.3, for instance). Since Fireshape is built ontop of the finite element library Firedrake, these formu-las must be written using the Unified Form Language(UFL). We refer to the tutorials on the website [27] formore details about Firedrake and UFL.4.4 Supplementary classesFireshape also includes a few extra classes to specify ad-ditional constraints, such as volume or perimeter con-straints on the domain Ω, or spectral constraints tocontrol the singular values of the transformation T . Formore details about these extra options, we refer to Fire-shape’s documentation and tutorials [28]. We have introduced Fireshape: an open-source and au-tomated shape optimization toolbox for the finite el-ement software Firedrake. Fireshape is based on the ireshape: a shape optimization toolbox for Firedrake 11 moving mesh method and allows user with minimalshape optimization knowledge to tackle challenging shapeoptimization problems constrained to PDEs. In partic-ular, Fireshape computes adjoint equations and shapederivatives in an automated fashion, allows decoupleddiscretizations of control and state variables, and givesfull access to Firedrake’s and PETSc’s discretizationand solver capabilities as well as to the extensive opti-mization library ROL.
Acknowledgements
The work of Florian Wechsung waspartly supported by the EPSRC Centre For Doctoral Train-ing in Industrially Focused Mathematical Modelling [grantnumber EP/L015803/1].
Contributions
Alberto Paganini and Florian Wechsung havecontributed equally to the development of the software Fire-shape and to this manuscript.
Conflict of interest statement
On behalf of all authors,the corresponding author states that there is no conflict ofinterest.
Replication of results
For reproducibility, we cite archivesof the exact software versions used to produce the resultsin this paper. All major Firedrake components have beenarchived on Zenodo [58]. An installation of Firedrake withcomponents matching those used to produce the results inthis paper can by obtained following the instructions at . The exact ver-sion of the Fireshape library used for these results has alsobeen archived at [59]. The latest version of the Fireshape li-brary can be found at https://github.com/Fireshape/Fireshape . References
1. Abramowski, T., Sugalski, K.: Energy saving proceduresfor fishing vessels by means of numerical optimization ofhull resistance. Zeszyty Naukowe Akademii Morskiej wSzczecinie (2017)2. Ahrens, J., Geveci, B., Law, C.: Paraview: An end-usertool for large data visualization. The visualization hand-book (2005)3. Allaire, G.: Conception optimale de structures,
Mathe-matics & Applications , vol. 58. Springer-Verlag, Berlin(2007)4. Allaire, G., Jouve, F., Toader, A.M.: A level-set methodfor shape optimization. C. R. Math. Acad. Sci. Paris (12), 1125–1130 (2002)5. Allaire, G., Pantz, O.: Structural optimization withFreeFem++. Struct. Multidiscip. Optim. (3), 173–181(2006)6. Alnæs, M.S., Logg, A., Ølgaard, K.B., Rognes, M.E.,Wells, G.N.: Unified Form Language: A Domain-specificLanguage for Weak Formulations of Partial DifferentialEquations. ACM Transactions on Mathematical Software (2), 9:1–9:37 (2014). DOI 10.1145/25666307. Amestoy, P.R., Duff, I.S., LExcellent, J.Y., Koster, J.:MUMPS: a general purpose distributed memory sparsesolver. In: International Workshop on Applied Parallel Computing, pp. 121–130. Springer (2000). DOI 10.1007/3-540-70734-4 168. Balay, S., Abhyankar, S., Adams, M.F., Brown, J., Brune,P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp,W.D., Karpeyev, D., Kaushik, D., Knepley, M.G., May,D.A., McInnes, L.C., Mills, R.T., Munson, T., Rupp, K.,Sanan, P., Smith, B.F., Zampini, S., Zhang, H., Zhang,H.: PETSc users manual. Tech. Rep. ANL-95/11 - Revi-sion 3.11, Argonne National Laboratory (2019)9. Balay, S., Gropp, W.D., McInnes, L.C., Smith, B.F.: Ef-ficient management of parallelism in object oriented nu-merical software libraries. In: E. Arge, A.M. Bruaset,H.P. Langtangen (eds.) Modern Software Tools in Scien-tific Computing, pp. 163–202. Birkh¨auser Press (1997)10. Bendsøe, M.P.: Optimal shape design as a material dis-tribution problem. Structural optimization (4), 193–202(1989)11. Bendsoe, M.P., Sigmund, O.: Topology Optimization:Theory, Methods and Applications. Springer (2004)12. Bernland, A., Wadbro, E., Berggren, M.: Shape opti-mization of a compression driver phase plug. SIAMJ. Sci. Comput. (1), B181–B204 (2019). DOI10.1137/18M1175768. URL https://doi.org/10.1137/18M1175768
13. Blank, L., Garcke, H., Farshbaf-Shaker, M.H., Styles, V.:Relating phase field and sharp interface approaches tostructural topology optimization. ESAIM Control Op-tim. Calc. Var. (4), 1025–1058 (2014)14. Chevalier, C., Pellegrini, F.: PT-SCOTCH: a tool for effi-cient parallel graph ordering. Parallel Computing (6),318–331 (2008)15. Ciarlet, P.G.: The finite element method for elliptic prob-lems. Society for Industrial and Applied Mathematics(SIAM), Philadelphia, PA (2002)16. Courtais, A., Lesage, F., Privat, Y., Frey, P., razak Lat-ifi, A.: Adjoint system method in shape optimization ofsome typical fluid flow patterns. In: A.A. Kiss, E. Zonder-van, R. Lakerveld, L. zkan (eds.) 29th European Sympo-sium on Computer Aided Process Engineering, ComputerAided Chemical Engineering , vol. 46, pp. 871 – 876. El-sevier (2019)17. Dalcin, L.D., Paz, R.R., Kler, P.A., Cosimo, A.: Paralleldistributed computing using Python. Advances in WaterResources (9), 1124–1139 (2011). New ComputationalMethods and Software Tools18. Dapogny, C., Frey, P., Omn`es, F., Privat, Y.: Geometricalshape optimization in fluid mechanics using FreeFem++.Struct. Multidiscip. Optim. (6), 2761–2788 (2018)19. Delfour, M.C., Zol´esio, J.P.: Shapes and geometries. Met-rics, analysis, differential calculus, and optimization, Ad-vances in Design and Control , vol. 22, second edn. So-ciety for Industrial and Applied Mathematics (SIAM),Philadelphia, PA (2011)20. Denis Rldzal, D.K.: Rapid optimization library, ver-sion 00 (2014). URL
21. Dokken, J.S., Funke, S.W., Johansson, A., Schmidt, S.:Shape optimization using the finite element method onmultiple meshes with nitsche coupling. SIAM Journal onScientific Computing (3), A1923–A1948 (2019)22. Dokken, J.S., Mitusch, S.K., Funke, S.W.: Automaticshape derivatives for transient PDEs in FEniCS and Fire-drake (2020)23. Dondl, P., Poh, P.S.P., Rumpf, M., Simon, S.: Simul-taneous elastic shape optimization for a domain split-ting in bone tissue engineering. Proc. A. (2227),2 Alberto Paganini, Florian Wechsung20180718, 17 (2019). DOI 10.1098/rspa.2018.0718. URL https://doi.org/10.1098/rspa.2018.0718
24. Eigel, M., Sturm, K.: Reproducing kernel hilbert spacesand variable metric algorithms in pde-constrained shapeoptimization. Optimization Methods and Software (2),268–296 (2018)25. Elman, H.C., Silvester, D.J., Wathen, A.J.: Finite ele-ments and fast iterative solvers: with applications in in-compressible fluid dynamics, second edn. Oxford Univer-sity Press, Oxford (2014)26. Farrell, P.E., Ham, D.A., Funke, S.W., Rognes, M.E.: Au-tomated derivation of the adjoint of high-level transientfinite element programs. SIAM J. Sci. Comput. (4),C369–C393 (2013)27. Firedrake, project website.
28. Fireshape, documentation. https://fireshape.readthedocs.io/en/latest/index.html
29. Garcke, H., Hecht, C., Hinze, M., Kahle, C.: Numericalapproximation of phase field based shape and topologyoptimization for fluids. SIAM Journal on Scientific Com-puting (4), A1846–A1871 (2015)30. Geuzaine, C., Remacle, J.F.: Gmsh: A 3-d finite elementmesh generator with built-in pre-and post-processing fa-cilities. International journal for numerical methods inengineering (11), 1309–1331 (2009)31. Gibson, T.H., Mitchell, L., Ham, D.A., Cotter, C.J.: Adomain-specific language for the hybridization and staticcondensation of finite element methods (2018). URL https://arxiv.org/abs/1802.00303
32. Ham, D.A., Mitchell, L., Paganini, A., Wechsung, F.: Au-tomated shape differentiation in the Unified Form Lan-guage. Struct. Multidiscip. Optim. (5), 1813–1820(2019)33. He, X., Li, J., Mader, C.A., Yildirim, A., Martins, J.R.:Robust aerodynamic shape optimizationfrom a circle toan airfoil. Aerospace Science and Technology , 48–61(2019)34. Hinze, M., Pinnau, R., Ulbrich, M., Ulbrich, S.: Opti-mization with PDE constraints, Mathematical Modelling:Theory and Applications , vol. 23. Springer, New York(2009)35. H¨ollig, K.: Finite element methods with B-splines. So-ciety for Industrial and Applied Mathematics (SIAM),Philadelphia, PA (2003)36. Homolya, M., Kirby, R.C., Ham, D.A.: Exposing and ex-ploiting structure: optimal code generation for high-orderfinite element methods (2017). URL http://arxiv.org/abs/1711.02473
37. Homolya, M., Mitchell, L., Luporini, F., Ham, D.A.:TSFC: a structure-preserving form compiler (2017). URL http://arxiv.org/abs/1705.003667
38. Iglesias, J.A., Sturm, K., Wechsung, F.: Two-dimensionalshape optimization with nearly conformal transforma-tions. SIAM Journal on Scientific Computing (6),A3807–A3830 (2018)39. Ito, K., Kunisch, K.: Lagrange multiplier approach tovariational problems and applications, Advances in De-sign and Control , vol. 15. Society for Industrial and Ap-plied Mathematics (SIAM), Philadelphia, PA (2008)40. Kelley, C.T.: Solving nonlinear equations with Newton’smethod,
Fundamentals of Algorithms , vol. 1. Society forIndustrial and Applied Mathematics (SIAM), Philadel-phia, PA (2003)41. Laurain, A.: A level set-based structural optimizationcode using FEniCS. Struct. Multidiscip. Optim. (3),1311–1334 (2018) 42. Lawrence Livermore National Laboratory: hypre : High Performance Preconditioners. https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods
43. Liu, Z., Korvink, J., Huang, R.: Structure topology op-timization: fully coupled level set method via femlab.Structural and Multidisciplinary Optimization (6),407–417 (2005)44. Lombardi, M., Parolini, N., Quarteroni, A., Rozza, G.:Numerical simulation of sailing boats: dynamics, fsi,and shape optimization. In: Variational Analysis andAerospace Engineering: Mathematical Challenges forAerospace Design, pp. 339–377. Springer (2012)45. Luporini, F., Ham, D.A., Kelly, P.H.J.: An algorithmfor the optimization of finite element integration loops.ACM Transactions on Mathematical Software , 3:1–3:26 (2017)46. Marini´c-Kragi´c, I., Vuˇcina, D., ´Curkovi´c, M.: Efficientshape parameterization method for multidisciplinaryglobal optimization and application to integrated shiphull shape optimization workflow. Computer-Aided De-sign , 61–75 (2016)47. Mitusch, S.K., Funke, S.W., Dokken, J.S.: dolfin-adjoint2018.1: automated adjoints for FEniCS and Firedrake.The Journal of Open Source Software (2019)48. Nocedal, J., Wright, S.J.: Numerical optimization, secondedn. Springer, New York (2006)49. Paganini, A., Wechsung, F., Farrell, P.E.: Higher-ordermoving mesh methods for PDE-constrained shape opti-mization. SIAM J. Sci. Comput. (4), A2356–A2382(2018)50. Rathgeber, F., Ham, D.A., Mitchell, L., Lange, M., Lu-porini, F., McRae, A.T.T., Bercea, G.T., Markall, G.R.,Kelly, P.H.J.: Firedrake: automating the finite elementmethod by composing abstractions. ACM Trans. Math.Softw. (3), 24:1–24:27 (2016)51. Schmidt, S., Ilic, C., Schulz, V., Gauger, N.R.: Three-dimensional large-scale aerodynamic shape optimizationbased on shape calculus. AIAA journal (11), 2615–2627 (2013)52. Schulz, V., Gherman, I.: One-shot methods for aero-dynamic shape optimization. In: MEGADESIGN andMegaOpt-German Initiatives for Aerodynamic Simula-tion and Optimization in Aircraft Design, pp. 207–220.Springer (2009)53. Sigmund, O.: A 99 line topology optimization code writ-ten in matlab. Struct. Multidiscip. Optim. (2), 120127(2001)54. Talischi, C., Paulino, G.H., Pereira, A., Menezes, I.F.:Polytop: a matlab implementation of a general topologyoptimization framework using unstructured polygonal fi-nite element meshes. Structural and MultidisciplinaryOptimization (3), 329–357 (2012)55. Wang, H., Doherty, J., Jin, Y.: Hierarchical surrogate-assisted evolutionary multi-scenario airfoil shape opti-mization. In: 2018 IEEE Congress on Evolutionary Com-putation (CEC), pp. 1–8. IEEE (2018)56. Wechsung, F.: Shape optimisation and robust solvers forincompressible flow. Ph.D. thesis, University of Oxford(2019)57. Wendland, H.: Scattered data approximation. CambridgeUniversity Press, Cambridge (2005)58. Software used in ’Fireshape: a shape optimization tool-box for Firedrake’ (2020). DOI 10.5281/zenodo.3823872.URL https://doi.org/10.5281/zenodo.3823872
59. Fireshape: a shape optimization toolbox for Firedrake(2020). DOI 10.5281/zenodo.3824907. URL59. Fireshape: a shape optimization toolbox for Firedrake(2020). DOI 10.5281/zenodo.3824907. URL