Fully and Semi-Automated Shape Differentiation in NGSolve
Peter Gangl, Kevin Sturm, Michael Neunteufel, Joachim Schöberl
FFully and Semi-Automated Shape Differentiation in
NGSolve
Peter Gangl ∗ , Kevin Sturm † , Michael Neunteufel ‡ and Joachim Schöberl § April 16, 2020
Abstract
In this paper we present a framework for automated shape differentiation in the finiteelement software
NGSolve . Our approach combines the mathematical Lagrangian approachfor differentiating PDE constrained shape functions with the automated differentiation ca-pabilities of
NGSolve . The user can decide which degree of automatisation is required andthus allows for either a more custom-like or black-box-like behaviour of the software.We discuss the automatic generation of first and second order shape derivatives for un-constrained model problems as well as for more realistic problems that are constrained bydifferent types of partial differential equations. We consider linear as well as nonlinear prob-lems and also problems which are posed on surfaces. In numerical experiments we verify theaccuracy of the computed derivatives via a Taylor test. Finally we present first and secondorder shape optimisation algorithms and illustrate them for several numerical optimisationexamples ranging from nonlinear elasticity to Maxwell’s equations.
Keywords: shape optimisation, automated differentiation, shape Newton method
Numerical simulation and shape optimisation tools to solve the problems have become an in-tegral part in the design process of many products. Starting out from an initial design, non-parametric shape optimisation techniques based on first and second order shape derivatives canassist in finding shapes of a product which are optimal with respect to a given objective func-tion. Examples include the optimal design of aircrafts [
28, 29 ] , optimal inductor design [ ] ,optimisation of microlenses [ ] , the optimal design of electric motors [ ] , applications frommechanical engineering [
2, 19 ] , multiphysics problems [ ] or electrical impedance tomography(EIT) in medical sciences to name only a few [ ] .Shape optimisation algorithms are based on the concept of shape derivatives. Let (cid:65) ⊂ (cid:80) ( R d ) a set of admissible shapes and (cid:74) : (cid:65) → R a shape function. Given an admissible shape Ω ∈ (cid:65) ∗ TU Graz, Steyrergasse 30, 8010 Graz, Austria, E-Mail: gangl(at)math.tugraz.at † TU Wien, Wiedner Hauptstr. 8-10, 1040 Vienna, Austria, E-Mail: kevin.sturm(at)tuwien.ac.at ‡ TU Wien, Wiedner Hauptstr. 8-10, 1040 Vienna, Austria, E-Mail: michael.neunteufel(at)tuwien.ac.at § TU Wien, Wiedner Hauptstr. 8-10, 1040 Vienna, Austria, E-Mail: joachim.schoeberl(at)tuwien.ac.at a r X i v : . [ m a t h . O C ] A p r P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel and a sufficiently smooth vector field V , we define the perturbed domain Ω t : = ( Id + t V )( Ω ) fora small perturbation parameter t >
0. The shape derivative is defined as D (cid:74) ( Ω )( V ) : = (cid:129) dd t (cid:74) ( Ω t ) (cid:139)(cid:12)(cid:12)(cid:12)(cid:12) t = = lim t (cid:38) (cid:74) ( Ω t ) − (cid:74) ( Ω ) t . (1.1)In most practically relevant applications, the objective functional depends on the shape ofa (sub-)domain via the solution of a partial differential equation (PDE). Thus, one is facing aproblem of PDE-constrained shape optimisation of the formmin ( Ω , u ) ∈(cid:65) × X J ( Ω , u ) s.t. ( Ω , u ) ∈ (cid:65) × X : e ( Ω ; u , v ) = v ∈ X . (1.2)Here, the second line represents the constraining boundary value problem posed on a Hilbertspace X , which we assume to be uniquely solvable for all admissible Ω ∈ (cid:65) . Denoting theunique solution for a given Ω ∈ (cid:65) by u Ω , we introduce the notation for the reduced functional (cid:74) ( Ω ) : = J ( Ω , u Ω ) .In order to be able to apply a shape optimisation algorithm to a given problem of this kind,the shape derivative (1.1) has to be computed, see the standard literature [
5, 32 ] or [ ] foran overview of different approaches. In the following we focus on computing the so-called vol-ume form of the shape derivative which in a finite element context is known to give a betterapproximation compared to the boundary form; see [
4, 15 ] .The convergence of shape optimisation algorithms can be speeded up by using second ordershape derivatives. Given two sufficiently smooth vector fields V , W and an admissible shape Ω ∈ (cid:65) , let Ω s , t : = ( Id + sV + tW )( Ω ) the perturbed domain. Then, the second order shapederivative is defined as D (cid:74) ( Ω )( V )( W ) : = (cid:129) d dsd t (cid:74) ( Ω s , t ) (cid:139)(cid:12)(cid:12)(cid:12)(cid:12) s , t = . (1.3)Second order information in Newton-type algorithms has been explored in the articles [
1, 8,22, 24, 31 ] . Since the computation of second order shape derivatives is more involved and er-ror prone, several authors have employed automatic differentiation (AD) tools, see e.g. [ ] and [ ] for two approaches based on the Unified Form Language (UFL) [ ] . In [ ] , the au-thors present a fully automated shape differentiation software which uses the transformationproperties on the finite element level. In [ ] (see also the earlier work [ ] ) the automatedderivatives are computed using UFL. The strategy of [ ] and [ ] differ in that, for the latter,the software computes an unsymmetric shape Hessian since it involves the term D (cid:74) ( Ω )( ∂ V W ) .Optionally the software allows to make the shape Hessian symmetric by requiring ∂ V W = [ ] where automated shape derivatives for transient PDEs in FEniCS andFiredrake are presented.In this paper we present an alternative framework for AD of PDE constrained problems oftype (1.2). There exist several approaches for the rigorous derivation of the shape derivative ofPDE-constrained shape functionals, see [ ] for an overview. The main idea, however, is always utomated Shape Optimisation in NGSolve t which now enters via the corresponding transformation and its gradient. It is shownin [ ] that the shape derivative for a nonlinear PDE-constrained shape optimisation problemcan be computed as the derivative of the Lagrangian with respect to the perturbation parame-ter. We will illustrate this systematic procedure for a number of different applications and utilisesymbolic differentiation provided by the finite element software package NGSolve [ ] to obtainthe shape derivative for different classes of PDE-constrained optimisation problems. NGSolve allows for a fast and efficient numerical solution of a large number of different boundary valueproblems. The aim of this paper is to extend
NGSolve by the possibility of semi-automatic andfully automatic shape differentiation and optimisation.Distinctly from previous approaches we cover the following two points: • a fully automated setting requiring as input the weak formulation of the constraint andthe cost function, • a semi-automated setting which offers a highly customizable user interface, but requiresmathematical background knowledge. Structure of the paper.
In Section 2 we give a brief introduction on how to solve a PDE in
NGSolve and present its built-in auto-differentiation capabilities. The introduced syntax willalso lay the foundation for the following sections. In Section 3 we present a first unconstrainedshape optimisation problem and show how to solve it in
NGSolve . For this purpose we showhow to compute the first and second order shape derivative in a semi-automated way. Section 4extends the preceding section by incorporating a PDE constraint. The strategy is illustratedby means of a simple Poisson equation. We also show how to treat the computation of shapederivatives when the PDE is defined on surfaces. While the semi-automated shape differentiationpresented in Sections 3 and 4 requires mathematical background knowledge, in Section 5 weshow how the shape derivatives can be computed in a fully automated fashion. In the last sectionof the paper we verify the computed formulas by a Taylor test, discuss optimisation algorithmsand present several numerical optimisation examples including nonlinear elasticity, Maxwell’sequations and Helmholtz’s equation.
NGSolve
In this section, we give a brief overview over the main concepts of the finite element software
NGSolve [ ] . We first describe the main principles for numerically solving boundary valueproblems in NGSolve before focusing on its built-in automatic differentiation capabilities. In thesubsequent sections of this paper, these ingredients will be combined to implement the shapederivative of unconstrained and PDE-constrained shape optimisation problems in an automatedway.
P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
NGSolve
In this section, we illustrate the syntax of
NGSolve using the python programming languagefor the Poisson equation with homogeneous Dirichlet conditions as a model problem. We referthe reader to the online documentationhttps: // ngsolve.org / docu / latest / for a more detailed description of the many features.The weak form of the model problem on a domain Ω ⊂ R d readsFind u ∈ H ( Ω ) : ˆ Ω ∇ u · ∇ w d x = ˆ Ω f w d x for all w ∈ H ( Ω ) . (2.1)We consider a ball of radius in two space dimensions centered at the point ( ) (cid:62) , i.e. Ω = B (( ) (cid:62) , 0.5 ) , and the right hand side is defined by f ( x , x ) = x ( − x )+ x ( − x ) .We will go through the steps for numerically solving this problem by the finite element method.We begin by importing the necessary functionalities and setting up a finite element mesh. from n g s o l v e import ∗ from netgen . geom2d import SplineGeometry geo = SplineGeometry ( ) geo . A d d C i r c l e ( ( 0 . 5 , 0 . 5 ) , 0 . 5 , bc = " c i r c l e " ) mesh = Mesh ( geo . GenerateMesh (maxh = mesh . Curve (3) The first line imports all modules from the package
NGSolve . The second line includes theSplineGeometry function which enables us to define a mesh via a geometric description, in ourcase a circle centered at ( ) (cid:62) of radius 0.5. Finally the mesh is generated in line 7 andin line 8 we specify that we want to use a curved finite element mesh for a more accurateapproximation of the geometry.Next in line 9 we define an H conforming finite element space of polynomial degree 3 andinclude Dirichlet boundary conditions on the boundary of the domain ∂ Ω (referenced by thestring “circle” that we assigned in line 5). On this space we define a trial function u in line11 and a test function w in line 12. These are purely symbolic objects which are used to defineboundary value problems in weak form. f e s = H1( mesh , o r d e r =
3, d i r i c h l e t = " c i r c l e " ) u = f e s . T r i a l F u n c t i o n ( ) w = f e s . T e s t F u n c t i o n ( ) For a more compact presentation later on, we define a coefficient function X which combinesthe three spatial components: X = C o e f f i c i e n t F u n c t i o n ( ( x , y , z ) )
Now, the left and right hand side of problem (2.1) can be conveniently defined as a bilinear orlinear form, respectively, on the finite elements space fes by the following lines. utomated Shape Optimisation in NGSolve L = LinearForm ( f e s ) f 1 = (2 ∗ X [ ] ∗ (1 − X [ ] ) + ∗ X [ ] ∗ (1 − X [ ] ) ) L += f 1 ∗ w ∗ dx a = B i l i n e a r F o r m ( f e s , symmetric = True ) a += grad ( u ) ∗ grad (w) ∗ dx We assemble the system matrix coming from the bilinear form a and the load vector comingfrom L and solve the corresponding system of linear equations. a . Assemble ( ) L . Assemble ( ) gfu = G r i d F u n c t i o n ( f e s ) gfu . vec . data = a . mat . I n v e r s e ( f e s . F r e e D o f s ( ) , i n v e r s e = " s p a r s e c h o l e s k y " ) ∗ L . vec
Draw ( gfu , mesh , " s t a t e " )
Here, gfu is defined as a
GridFunction over the finite element space fes . The Dirichlet con-ditions are incorporated into the direct solution of the linear system and the numerical solutionis drawn in the graphical user interface. The numerical solution is depicted in Figure 1.
P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
NGSolve In NGSolve , symbolic expressions are stored in expression trees, see Figure 2 for an example.It is possible to differentiate an expression expr with respect to a variable var appearing in expr into a direction dir by the command expr.Diff(var, dir) .Mathematically this line corresponds to the directional derivative of g: = expr at x : = var indirection v : = d i r , that is, D g ( x )( v ) . (2.2)When calling the Diff command for expr , the expression tree of expr is gone through node bynode, and for each node the corresponding differentiation rules such as product rule or chainrule are applied. When a node represents the variable with respect to which the differentiationis carried out, it is replaced by the direction dir of differentiation.Figure 2 shows the differentiation of the expression expr = + with respect to x intothe direction given by v : v = Parameter (1) expr = ∗ x ∗ x + ∗ y dexpr = expr . D i f f ( x , v ) p r i n t ( expr ) p r i n t ( dexpr ) The output of print(expr) reads coef binary operation ’+’, realcoef binary operation ’*’, realcoef scale 2, realcoef coordinate x, realcoef coordinate x, realcoef scale 3, realcoef coordinate y, real which translates to 2 x ∗ x + y and corresponds to the expression tree depicted in Figure 2(a).The output of print(dexpr) reads coef binary operation ’+’, realcoef binary operation ’+’, realcoef binary operation ’*’, realcoef scale 2, realcoef N5ngfem28ParameterCoefficientFunctionE, realcoef coordinate x, realcoef binary operation ’*’, realcoef scale 2, realcoef coordinate x, realcoef N5ngfem28ParameterCoefficientFunctionE, realcoef scale 3, realcoef 0, real utomated Shape Optimisation in NGSolve + *
2x x + + *
2v x *
2x v (a) (b)Figure 2: Illustration of
Diff command for example expr = + . (a) Expression tree for expr . (b) Expression tree for expression obtained by call of expr.Diff(x, dir) .which translates to ( v ∗ x + x ∗ v ) + ∗ GridFunction s as well as trial and test functions do not depend on thespatial variables x , y , z . While trial and test functions are purely symbolic objects, GridFunction srepresent functions in the finite element space. The code segments u = f e s . T r i a l F u n c t i o n ( ) w = f e s . T e s t F u n c t i o n ( ) g f = G r i d F u n c t i o n ( f e s ) g f . S e t ( x ∗ x ∗ y ) p r i n t ( " D i f f u wrt x " , u . D i f f ( x ) ) p r i n t ( " D i f f w wrt x " , w. D i f f ( x ) ) p r i n t ( " D i f f g f wrt x " , g f . D i f f ( x ) ) will give the following output: Diff u wrt x: ConstantCF, val = 0Diff w wrt x: ConstantCF, val = 0Diff gf wrt x: ConstantCF, val = 0
We will illustrate the steps to be taken in order to obtain the shape derivative of a shape functionin a semi-automatic way for a simple shape optimisation problem. For Ω ⊂ R d bounded and openand a continuously differentiable function f ∈ C ( R d ) , we consider the shape differentiation of P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel the shape function (cid:74) ( Ω ) = ˆ Ω f ( x ) d x . (3.1)Clearly the minimiser of (cid:74) over all measurable sets in R d is given by Ω ∗ = { x ∈ R d : f ( x ) < } .We also refer to [ ] for the computations of first and second order variations of type (3.1)where Ω is a submanifold of R d . Given a vector field V ∈ C ( R d ) d , we define the transformation T t ( x ) : = ( Id + t V )( x ) , x ∈ R d , t ≥ Definition 3.1.
The first order shape derivative of a shape function (cid:74) at Ω in direction V ∈ C ( R d ) d is defined by D (cid:74) ( Ω )( V ) = lim t (cid:38) (cid:74) ( T t ( Ω )) − (cid:74) ( Ω ) t . (3.2) Using the transformation y = T t ( x ) and the notation F t : = ∂ T t = I + t ∂ V for the Jacobian ofthe transformation T t , we get for (cid:74) as in (3.1), (cid:74) ( Ω t ) = ˆ Ω t f ( x (cid:48) ) d x (cid:48) = ˆ Ω ( f ◦ T t )( x ) det ( F t ( x )) d x . (3.3)Now let us explain how to compute the shape derivative of (cid:74) . Denoting G ( T t , F t ) : = ˆ Ω ( f ◦ T t )( x ) det ( F t ( x )) d x , (3.4)the chain rule gives (formally) dd t (cid:74) ( Ω t ) (cid:12)(cid:12)(cid:12)(cid:12) t = = dd t G ( T t , F t ) (cid:12)(cid:12)(cid:12)(cid:12) t = = (cid:129) d Gd T t d T t d t + d Gd F t d F t d t (cid:139)(cid:12)(cid:12)(cid:12)(cid:12) t = .Using that dT t d t ( x ) = V ( x ) and dF t d t ( x ) = ∂ V ( x ) , we get for the shape derivative d (cid:74) ( Ω )( V ) = dd t (cid:74) ( Ω t ) (cid:12)(cid:12)(cid:12)(cid:12) t = = (cid:129) d Gd T t V + d Gd F t ∂ V (cid:139)(cid:12)(cid:12)(cid:12)(cid:12) t = .This is the form we use for defining the first order shape derivative in NGSolve . Note that aLipschitz vector field is differentiable almost everywhere and hence ∂ V ( x ) is defined almosteverywhere and bounded.Given the function f ( x , x ) = ( x − ) / a + ( x − ) / b − R with a = b = / a and R = utomated Shape Optimisation in NGSolve f = ( ( X [ ] − / ∗∗ + (1.3 ∗ (X [ ] − ∗∗ − ∗∗ F = Id (2) G_f = f ∗ Det ( F ) ∗ dx Here, we introduce the symbol F and assign to it the value of the identity matrix in line 42. Thisallows us to differentiate with respect to F . Then we define the function G of (3.4) in line 43.The shape derivative is a bounded linear functional on a space of vector fields. We introduce avector-valued finite element space VEC and define the object representing the shape derivative dJOmega_f as a linear functional on
VEC . In line 48, we differentiate with respect to the spatialvariables in the direction given by V . Note that X is the coefficient function we introduced inline 13. In line 49, we deal with the differentiation with respect to F . VEC = VectorH1 ( mesh , o r d e r =
1, d i r i c h l e t = " " ) V = VEC . T e s t F u n c t i o n ( ) dJOmega_f = LinearForm (VEC) dJOmega_f += G_f . D i f f (X , V) dJOmega_f += G_f . D i f f ( F , grad (V) )
Remark 3.2.
Defining ξ t : = det ( F t ) and using dd t ξ t | t = = div V , it holds d Gd F t d F t d t (cid:12)(cid:12)(cid:12)(cid:12) t = = d Gd ξ t d ξ t d F t d F t d t (cid:12)(cid:12)(cid:12)(cid:12) t = = d Gd ξ t d ξ t d t (cid:12)(cid:12)(cid:12)(cid:12) t = = d Gd ξ t div V (cid:12)(cid:12)(cid:12)(cid:12) t = = ˆ Ω f div V d x .Therefore, we obtain for the first order shape derivative the well-known formula D (cid:74) ( Ω )( V ) = ˆ Ω ∇ f · V + f div V d x .Finally if Ω is smooth enough (for instance C ), then the shape derivative is given by D (cid:74) ( Ω )( V ) = ˆ ∂ Ω f V · n d s , (3.5)where n denotes the outward pointing normal along ∂ Ω . For Ω and f as in the previous section we consider (cid:74) bnd ( Ω ) = ˆ ∂ Ω f ( x ) d x . (3.6)Then we get (cid:74) bnd ( Ω t ) = ˆ ∂ Ω t f ( x (cid:48) ) d s x (cid:48) = ˆ ∂ Ω ( f ◦ T t )( x ) det ( F t ( x )) | F t ( x ) −(cid:62) n ( x ) | d s x , (3.7)see e.g. [
32, Prop. 2.47 ] , with the outer unit normal vector n and | · | denoting the Euclideannorm. Again, the shape derivative can be computed as the total derivative of this expressionwith respect to the parameter t . In NGSolve , the only difference lies in the necessity to use thetrace of the gradient of a test vector field V .0 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel G_f_bnd = f ∗ Det ( F ) ∗ Norm( Inv ( F ) . t r a n s ∗ s p e c i a l c f . normal (2) ) ∗ ds dJOmega_f_bnd = LinearForm (VEC) dJOmega_f_bnd += G_f_bnd . D i f f (X , V) dJOmega_f_bnd += G_f_bnd . D i f f ( F , grad (V) . Trace ( ) )
For second order shape derivatives, we consider perturbations of the form T s , t ( x ) = ( Id + sV + tW )( x ) , x ∈ R d ,for s , t ≥ Ω s , t : = T s , t ( Ω ) . Definition 3.3.
The second order shape derivative of a shape function (cid:74) at Ω in direction ( V , W ) ∈ C ( R d ) d × C ( R d ) d is defined by D (cid:74) ( Ω )( V )( W ) = d dsd t (cid:74) ( Ω s , t ) (cid:12)(cid:12)(cid:12)(cid:12) s = t = . (3.8) Remark 3.4.
We remark that if (cid:74) is smooth enough, the second order derivative as defined in(3.8) is symmetric by definition: D (cid:74) ( Ω )( V )( W ) = D (cid:74) ( Ω )( W )( V ) . (3.9)We stress that this derivative is not the same as the shape derivative obtained by repeated shapedifferentiation, that is, it does not coincide with (see, e.g., [
6, Chap. 9, Sec. 6 ] ) d (cid:74) ( Ω )( V )( W ) : = lim t (cid:38) D (cid:74) ( T Wt ( Ω ))( V ) − D (cid:74) ( Ω )( V ) t (3.10)which is in general asymmetric. However, in NGSolve we compute directly the second derivativeas defined in (3.8). However, this derivative is only symmetric if ∂ V W = d (cid:74) ( Ω )( V )( W ) = D (cid:74) ( Ω )( V )( W ) + D (cid:74) ( Ω )( ∂ V W ) . (3.11)In NGSolve , when repeating the shape differentiation procedure introduced in Section 3.1, wecompute directly the second order shape derivative as defined in (3.8). Here, we exploit that
GridFunction s are independent of the spatial coordinates, see also Section 2.2.Similarly to the computations of the first derivative, we use the notation F s , t : = ∂ T s , t = I + s ∂ V + t ∂ W . For the shape function defined in (3.1) we get d dsd t (cid:74) ( Ω s , t ) (cid:12)(cid:12)(cid:12)(cid:12) s = t = = d dsd t ˆ Ω s , t f ( x ) d x (cid:12)(cid:12)(cid:12)(cid:12) s = t = = d dsd t ˆ Ω ( f ◦ T s , t )( x ) det ( F s , t ( x )) d x (cid:12)(cid:12)(cid:12)(cid:12) s = t = . utomated Shape Optimisation in NGSolve G ( T s , t , F s , t ) = ´ Ω ( f ◦ T s , t )( x ) det ( F s , t ( x )) d x , we get d dsd t (cid:74) ( Ω s , t ) (cid:12)(cid:12)(cid:12)(cid:12) s = t = = d dsd t G ( T s , t , F s , t ) (cid:12)(cid:12)(cid:12)(cid:12) s = t = = dds (cid:18) d Gd T s , t d T s , t d t + d Gd F s , t d F s , t d t (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) s = t = .Using that d T s , t dsd t = d F s , t dsd t =
0, we get further d dsd t (cid:74) ( Ω s , t ) (cid:12)(cid:12)(cid:12)(cid:12) s = t = = dds (cid:18) d Gd T s , t (cid:19) d T s , t d t + dds (cid:18) d Gd F s , t (cid:19) d F s , t d t (cid:12)(cid:12)(cid:12)(cid:12) s = t = = (cid:18) d Gd T s , t d T s , t ds + d Gd F s , t d T s , t d F s , t ds (cid:19) d T s , t d t + (cid:18) d Gd T s , t d F s , t d T s , t ds + d Gd F s , t d F s , t ds (cid:19) d F s , t d t (cid:12)(cid:12)(cid:12)(cid:12) s = t = .(3.12)Formula (3.12) is used for the automatic derivation of the second order shape derivative in NGSolve . Using dT s , t ds ( x ) = V ( x ) , dT s , t d t ( x ) = W ( x ) and dF s , t ds ( x ) = ∂ V ( x ) , dF s , t d t ( x ) = ∂ W ( x ) , weget d dsd t (cid:74) ( Ω s , t ) (cid:12)(cid:12)(cid:12)(cid:12) s = t = = (cid:18) d Gd T s , t V + d Gd F s , t d T s , t ∂ V (cid:19) W + (cid:18) d Gd T s , t d F s , t V + d Gd F s , t ∂ V (cid:19) ∂ W (cid:12)(cid:12)(cid:12)(cid:12) s = t = .(3.13) Remark 3.5.
We remark that the formula (3.13) can be evaluated explicitly and reads D (cid:74) ( Ω )( V , W ) = ˆ Ω ∇ f V · W + ∇ f · W div V + ∇ f · V div W + f div V div W − f ∂ V (cid:62) : ∂ W d x .Formula (3.13) can be implemented in NGSolve as follows: d2JOmega_f = B i l i n e a r F o r m (VEC) W = VEC . T r i a l F u n c t i o n ( ) d2JOmega_f += (G_f . D i f f (X , W) + G_f . D i f f ( F , grad (W) ) ) . D i f f (X , V) d2JOmega_f += (G_f . D i f f (X , W) + G_f . D i f f ( F , grad (W) ) ) . D i f f ( F , grad (V) )
Notice that since W is a trial function it is not affected by the differentiation with respect to X ,see Section 2.2. In the same fashion, second order derivatives of boundary integrals of the form(3.6) can be computed. d2JOmega_f_bnd = B i l i n e a r F o r m (VEC) d2JOmega_f_bnd += (G_f_bnd . D i f f (X , W) + G_f_bnd . D i f f ( F , grad (W) . Trace ( ) ) ) . D i f f (X , V) d2JOmega_f_bnd += (G_f_bnd . D i f f (X , W) + G_f_bnd . D i f f ( F , grad (W) . Trace ( ) ) ) . D i f f ( F , grad (V) . Trace ( ) )
In this section, we describe the automatic derivation of the shape derivative for the followingtype of equality constrained shape optimisation problems:min ( Ω , u ) J ( Ω , u ) (4.1)2 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel subject to ( Ω , u ) ∈ (cid:65) × X solves e ( Ω , u ) =
0, (4.2)where e : (cid:65) × X → X ∗ with e ( Ω , · ) : X ( Ω ) → X ( Ω ) ∗ represents an abstract PDE constraint with X = ∪ Ω ∈(cid:65) X ( Ω ) being the union of Banach spaces X ( Ω ) and (cid:65) a set of admissible shapes. Forany given Ω ∈ (cid:65) we assume the PDE constraint (4.2) to admit a unique solution which wedenote by u Ω . Moreover, let (cid:74) ( Ω ) : = J ( Ω , u Ω ) denote the reduced cost functional. One way toeliminate the constraint (4.2) is to introduce a Lagrangian (cid:76) ( Ω , u , p ) : = J ( Ω , u ) + 〈 e ( Ω , u ) , p 〉 . (4.3)Now an initial shape Ω is perturbed by a family of transformations T t , resulting in a new shape Ω t : = T t ( Ω ) . Transforming back to the initial shape Ω leads to the Lagrangian: G ( t , u , p ) : = (cid:76) ( T t ( Ω ) , Φ t ( u ) , Φ t ( p )) , u , p ∈ X ( Ω ) , (4.4)where Φ t : X ( Ω ) → X ( Ω t ) is a bijective mapping. Here the transformation Φ t depends on thedifferential operator involved. For instance • if X ( Ω ) = H ( Ω ) , then Φ t ( u ) = u ◦ T − t , • if X ( Ω ) = H ( curl , Ω ) , then Φ t ( u ) = ∂ T −(cid:62) t ( u ◦ T − t ) , • if X ( Ω ) = H ( div, Ω ) , then Φ t ( u ) = ( ∂ T t ) ∂ T t ( u ◦ T − t ) .Now the shape differentiability of (4.1)–(4.2) is reduced to proving that (see [ ] ) D (cid:74) ( Ω )( V ) = dd t G ( t , u t , 0 ) | t = = ∂ t G ( u , p ) , (4.5)where u t : = u t ◦ T t and u t ∈ X ( Ω t ) solves e ( Ω t , u t ) =
0. The verification of (4.5) can be ac-complished by different methods, see e.g. [ ] for an overview. However, we remark that (4.5)holds for a large class of nonlinear PDE constrained shape optimisation; see [ ] .The rest of this section is organised as follows: We introduce a model problem, which isthe minimisation of a tracking-type cost functional subject to Poisson’s equation in Section 4.1.We illustrate how the first and second order shape derivative for this PDE-constrained modelproblem can be obtained in NGSolve in Sections 4.2 and 4.3. Finally, we also briefly discuss theextension to partial differentiation equations on surfaces.
We will illustrate the derivation of the first and second order shape derivative for the minimi-sation of a tracking-type cost functional subject to Poisson’s equation on the unknown domain Ω . Let d = f , u d ∈ H ( R d ) and (cid:65) ⊂ (cid:80) ( R d ) a set of admissible shapes. We consider theproblem min ( Ω , u ) J ( Ω , u ) = ˆ Ω | u − u d | d x (4.6a) utomated Shape Optimisation in NGSolve ( Ω , u ) ∈ (cid:65) × H ( Ω ) solve 〈 e ( Ω , u ) , ψ 〉 : = ˆ Ω ∇ u · ∇ ψ d x − ˆ Ω f ψ d x = ψ ∈ H ( Ω ) . (4.6b)The Lagrangian is given by (cid:76) ( Ω , ϕ , ψ ) : = ˆ Ω | ϕ − u d | d x + ˆ Ω ∇ ϕ · ∇ ψ d x − ˆ Ω f ψ d x . (4.7)Given an admissible shape Ω , a vector field V ∈ C ( R d ) d and t > Ω t : = ( Id + t V )( Ω ) the perturbed domain. Therefore the parametrised Lagrangian is given by G ( t , ϕ , ψ ) : = (cid:76) ( T t ( Ω ) , ϕ ◦ T − t , ψ ◦ T − t ) , ϕ , ψ ∈ H ( Ω ) . (4.8)Changing variables yields G ( t , ϕ , ψ ) = ˆ Ω | ϕ − u td | det ( F t ) d x + ˆ Ω ( F −(cid:62) t ∇ ϕ ) · ( F −(cid:62) t ∇ ψ ) det ( F t ) d x − ˆ Ω f t ψ det ( F t ) d x (4.9) = : ˜ G ( T t , F t , ϕ , ψ ) ,where u td = u d ◦ T t and f t = f ◦ T t . Recall that, for a given Ω ∈ (cid:65) , u Ω denotes the correspondingunique solution to (4.6b) and (cid:74) ( Ω ) the reduced cost functional, (cid:74) ( Ω ) : = J ( Ω , u Ω ) . Let u t ∈ H ( Ω ) the solution of the perturbed state equation brought back to the original domain Ω , thatis, u t ∈ H ( Ω ) is the unique solution to ∂ ψ G ( t , u t , 0 )( ψ ) = ψ ∈ H ( Ω ) . (4.10)Note that, for u t defined by (4.10), it holds (cid:74) ( Ω t ) = G ( t , u t , ψ ) for all ψ ∈ H ( Ω ) and thereforealso D (cid:74) ( Ω )( V ) = dd t G ( t , u t , ψ ) for all ψ ∈ H ( Ω ) .It can easily be shown that (4.5) holds and thus the shape derivative in the direction of avector field V ∈ C ( R ) d is given by D (cid:74) ( Ω )( V ) = ∂ t G ( u , p ) ,where p ∈ H ( Ω ) denotes the adjoint state and is defined as the unique solution p ∈ H ( Ω ) to ∂ ϕ G ( u , p )( ˆ ϕ ) = ϕ ∈ H ( Ω ) , (4.11)or explicitly ˆ Ω ∇ ˆ ϕ · ∇ p d x = − ˆ Ω ( u − u d ) ˆ ϕ d x for all ˆ ϕ ∈ H ( Ω ) . (4.12)4 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
By the discussion above, the first order shape derivative is given by ∂ t G ( u , p ) with G definedin (4.9) and u and p the unique solutions to the boundary value problems (4.6b) and (4.12),respectively.Writing ˜ G ( T t , F t ) : = ˜ G ( T t , F t , u , p ) = G ( t , u , p ) we obtain in analogy to the unconstrainedproblem D (cid:74) ( Ω )( V ) = dd t (cid:74) ( Ω t ) (cid:12)(cid:12)(cid:12)(cid:12) t = = (cid:129) d Gd T t V + d Gd F t ∂ V (cid:139)(cid:12)(cid:12)(cid:12)(cid:12) t = .We can compute explicitly d Gd F t | t = ∂ V = ˆ Ω div ( V )( u − u d ) − ( ∂ V + ∂ V (cid:62) ) ∇ u · ∇ p − div ( V ) ∇ u · ∇ p − f p div ( V ) d x ,(4.13) d Gd T t | t = V = ˆ Ω − ( u − u d ) ∇ u d · V − ∇ f · V p d x . (4.14)Now we are in a position to compute the first order shape derivative for the PDE-constrainedshape optimisation problem (4.6) in NGSolve . After solving the state equation as shown inSection 2.1, the adjoint equation can be solved as follows. ud = X [ ] ∗ (1 − X [ ] ) ∗ X [ ] ∗ (1 − X [ ] ) d e f Cost ( u ) : r e t u r n (u − ud ) ∗∗ ∗ Det ( F ) ∗ dx gfp = G r i d F u n c t i o n ( f e s ) dCostdu = LinearForm ( f e s ) dCostdu += Cost ( gfu ) . D i f f ( gfu , w) dCostdu . Assemble ( ) gfp . vec . data = − a . mat . I n v e r s e ( f e s . F r e e D o f s ( ) , i n v e r s e = " s p a r s e c h o l e s k y " ) . T ∗ dCostdu . vec Draw( gfp , mesh , " a d j o i n t " )
We can now define the Lagrangian (4.9) such that the shape derivative can be obtained by thesame procedure as in the unconstrained setting. Note that lines 82–83 coincide with lines 48–49. d e f Equation (u ,w) : r e t u r n ( ( Inv ( F ) . t r a n s ∗ grad ( u ) ) ∗ ( Inv ( F ) . t r a n s ∗ grad (w) ) − f ∗ w) ∗ Det ( F ) ∗ dx G_pde = Cost ( gfu ) + Equation ( gfu , gfp ) dJOmega_pde = LinearForm (VEC) dJOmega_pde += G_pde . D i f f (X , V) dJOmega_pde += G_pde . D i f f ( F , grad (V) ) utomated Shape Optimisation in NGSolve Let us introduce the notation G V , W ( s , t , ϕ , ψ ) = ˆ Ω ( F −(cid:62) s , t ∇ ϕ ) · ( F −(cid:62) s , t ∇ ψ ) det ( F s , t ) d x − ˆ Ω f ◦ T s , t ψ det ( F s , t ) d x (cid:124) (cid:123)(cid:122) (cid:125) = : 〈 E V , W ( s , t ) ϕ , ψ 〉 (4.15) + ˆ Ω | ϕ − u d ◦ T s , t | det ( F s , t ) d x (cid:124) (cid:123)(cid:122) (cid:125) = : J V , W ( s , t ) , (4.16)where T s , t ( x ) = x + sV ( x ) + tW ( x ) and F s , t : = ∂ T s , t . We observe that (cid:74) ( T s , t ( Ω )) = G V , W ( s , t , u s , t , p s , t ) (4.17)with ( u s , t , p s , t ) ∈ H ( Ω ) × H ( Ω ) being the solution to ∂ p G V , W ( s , t , u s , t , 0 )( ϕ ) = ϕ ∈ H ( Ω ) , ∂ u G V , W ( s , t , u s , t , p s , t )( ψ ) = ψ ∈ H ( Ω ) (4.18)for s , t ≥
0. In case t = u s : = u s , t | t = and p s : = p s , t | t = and similarly for t = s = u : = u s , t | s = t = and p : = p s , t | s = t = . Therefore, consecutive differentiation of (4.17) firstwith respect to t at zero and then with respect to s at zero yields D (cid:74) ( Ω )( V )( W ) = d dsd t G V , W ( s , t , u s , t , p s , t ) | s = t = = dds ∂ t G V , W ( s , 0, u s , p s ) | s = = ∂ s ∂ t G V , W (
0, 0, u , p ) + ∂ u ∂ t G V , W (
0, 0, u , p )( ∂ s u ) + ∂ p ∂ t G V , W (
0, 0, u , p )( ∂ s p ) ,(4.19)where ∂ s u ∈ H ( Ω ) solves the material derivative equation ∂ u ∂ p G V , W (
0, 0, u , 0 )( ψ )( ∂ s u ) = − ∂ s ∂ p G V , W (
0, 0, u , 0 )( ψ ) for all ψ ∈ H ( Ω ) , (4.20)or equivalently 〈 ∂ u E V , W (
0, 0 )( ∂ s u ) , ψ 〉 = −〈 ∂ s E V , W (
0, 0 ) u , ψ 〉 for all ψ ∈ H ( Ω ) . (4.21)Similarly the function ∂ s p ∈ H ( Ω ) solves the material derivative equation ∂ p ∂ u G V , W (
0, 0, u , p )( ψ )( ∂ s p ) = − ∂ u G V , W (
0, 0, u , p )( ψ )( ∂ s u ) − ∂ s ∂ u G V , W (
0, 0, u , p )( ψ ) (4.22)for all ψ ∈ H ( Ω ) . Formally (4.20) and (4.22) can be written as an operator equation (cid:129) ∂ u G V , W (
0, 0, u , p ) ∂ p ∂ u G V , W (
0, 0, u , p ) ∂ u ∂ p G V , W (
0, 0, u , p ) (cid:139) (cid:129) ∂ s u ∂ s p (cid:139) = − (cid:129) ∂ s ∂ u G V , W (
0, 0, u , p ) ∂ s ∂ p G V , W (
0, 0, u , 0 ) (cid:139) . (4.23)So to evaluate the second derivative (4.19) in some direction ( V , W ) we have to solve the system(4.23).This is realised in NGSolve by setting up a combined finite element space which we denoteby X2 . We define trial and test functions as well as grid functions representing the deformationvector fields V and W , which we initialise with some functions.6 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel X2 = FESpace ( [ f e s , f e s ] ) dsu , dsp = X2 . T r i a l F u n c t i o n ( ) uTest , pT es t = X2 . T e s t F u n c t i o n ( ) gfV = G r i d F u n c t i o n (VEC) gfW = G r i d F u n c t i o n (VEC) gfV . S e t ( ( X [ ] ∗ X [ ] ∗ X [ ] ∗ exp (X [ ] ) , X [ ] ∗ X [ ] ∗ X [ ] ∗ exp (X [ ] ) ) ) gfW . S et ( ( X [ ] ∗ X [ ] ∗ X [ ] ∗ exp (X [ ] ) , X [ ] ∗ X [ ] ∗ X [ ] ∗ exp (X [ ] ) ) ) We define a 2 × × shapeHessLag2 = B i l i n e a r F o r m ( X2 ) shapeGradLag2 = LinearForm ( X2 ) shapeHessLag2 += ( G_pde . D i f f ( gfu , uTest ) ) . D i f f ( gfu , dsu ) shapeHessLag2 += ( G_pde . D i f f ( gfu , uTest ) ) . D i f f ( gfp , dsp ) shapeHessLag2 += ( G_pde . D i f f ( gfp , p Te st ) ) . D i f f ( gfu , dsu ) shapeGradLag2 += ( G_pde . D i f f ( gfu , uTest ) ) . D i f f ( F , grad ( gfV ) ) shapeGradLag2 += ( G_pde . D i f f ( gfu , uTest ) ) . D i f f (X , gfV ) shapeGradLag2 += ( G_pde . D i f f ( gfp , p Te st ) ) . D i f f ( F , grad ( gfV ) ) shapeGradLag2 += ( G_pde . D i f f ( gfp , p Te st ) ) . D i f f (X , gfV ) We can solve this combined system for ∂ s u and ∂ s p and access and visualise the two compo-nents in the following way: gfCombined2 = G r i d F u n c t i o n ( X2 ) shapeHessLag2 . Assemble ( ) shapeGradLag2 . Assemble ( ) gfCombined2 . vec . data = shapeHessLag2 . mat . I n v e r s e ( X2 . F r e e D o f s ( ) , i n v e r s e = " umfpack" ) ∗ shapeGradLag2 . vec gfdsu = G r i d F u n c t i o n ( f e s ) g f d s p = G r i d F u n c t i o n ( f e s ) gfdsu . vec . data = gfCombined2 . components [ ] . vec g f d s p . vec . data = gfCombined2 . components [ ] . vec Draw( gfdsu , mesh , " dsu " )
Draw( gfdsp , mesh , " dsp " )
In order to obtain the second order shape derivative in the direction given by ( V , W ) , it remainsto evaluate the term (4.19). We define the three terms of (4.19) as bilinear forms, assemblethem and perform vector-matrix-vector multiplications: w1 = f e s . T r i a l F u n c t i o n ( ) q1 = f e s . T r i a l F u n c t i o n ( ) shapeHess11 = B i l i n e a r F o r m (VEC) shapeHess11 += ( G_pde . D i f f ( F , grad (W) ) + G_pde . D i f f (X , W) ) . D i f f ( F , grad (V) ) utomated Shape Optimisation in NGSolve shapeHess11 += ( G_pde . D i f f ( F , grad (W) ) + G_pde . D i f f (X , W) ) . D i f f (X , V) shapeHess11 . Assemble ( ) shapeHess12 = B i l i n e a r F o r m ( t r i a l s p a c e = f e s , t e s t s p a c e = VEC) shapeHess12 += ( G_pde . D i f f ( F , grad (V) ) + G_pde . D i f f (X , V) ) . D i f f ( gfu , w1) shapeHess12 . Assemble ( ) shapeHess13 = B i l i n e a r F o r m ( t r i a l s p a c e = f e s , t e s t s p a c e = VEC) shapeHess13 += ( G_pde . D i f f ( F , grad (V) ) + G_pde . D i f f (X , V) ) . D i f f ( gfp , q1 ) shapeHess13 . Assemble ( ) av = gfV . vec . C r e a t e V e c t o r ( ) av . data = shapeHess11 . mat ∗ gfV . vec adsu = gfV . vec . C r e a t e V e c t o r ( ) adsu . data = shapeHess12 . mat ∗ gfdsu . vec adsp = gfV . vec . C r e a t e V e c t o r ( ) adsp . data = shapeHess13 . mat ∗ g f d s p . vec d2J = I n n e r P r o d u c t (gfW . vec , av ) + I n n e r P r o d u c t (gfW . vec , adsu ) + I n n e r P r o d u c t (gfW .vec , adsp )
The automated shape differentiation is not restricted to partial differential equations on domains Ω , but is readily extended to surface PDEs. We consider a two dimensional closed surface M ⊂ R and denote by n the normal field along M . Let u d ∈ H ( R d ) be given and define J ( M , u ) = ˆ M | u − u d | d s , (4.24)where u ∈ H ( M ) solves the surface equation ˆ M ∇ M u · ∇ M ϕ + u ϕ d s = ˆ M f ϕ d s for all ϕ ∈ H ( M ) , (4.25)where ∇ M ϕ denotes the tangential gradient of ϕ ; see [
6, p.493, Def.5.1 ] . We assume that thefunction f ∈ H ( R ) is given. The Lagrangian is given by (cid:76) ( M , ϕ , ψ ) : = ˆ M | ϕ − u d | d s + ˆ M ∇ M ϕ · ∇ M ψ + ϕψ d s − ˆ M f ψ d s . (4.26)As in the previous section we fix an admissible shape M and let M t : = ( Id + t V )( M ) be a smallperturbation of M by means of a vector field V ∈ C ( R d ) d for t > G ( t , ϕ , ψ ) : = (cid:76) ( T t ( M ) , ϕ ◦ T − t , ψ ◦ T − t ) , ϕ , ψ ∈ H ( M ) . (4.27)Define the density ω ( F t ) : = det ( F t ) | F −(cid:62) t n | . Changing variables and using ( ∇ M t ϕ ) ◦ T t = B ( F t ) ∇ M ( ϕ ◦ T t ) , B ( F t ) = (cid:18) I − F −(cid:62) t n | F −(cid:62) t n | ⊗ F −(cid:62) t n | F −(cid:62) t n | (cid:19) F −(cid:62) t (4.28)8 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel yields G ( t , ϕ , ψ ) = ˆ M | ϕ − u td | ω ( F t ) d s + ˆ M (( B ( F t ) ∇ ϕ ) · ( B ( F t ) ∇ ψ ) + ϕψ ) ω ( F t ) d s − ˆ M f t ψ ω ( F t ) d s , (4.29)where u td = u d ◦ T t and f t = f ◦ T t .Writing ˜ G ( T t , F t ) : = G ( t , u , p ) we obtain in analogy to the domain case D (cid:74) ( Ω )( V ) = (cid:129) d Gd T t V + d Gd F t ∂ V (cid:139)(cid:12)(cid:12)(cid:12)(cid:12) t = . (4.30)We can compute explicitly d Gd F t | t = V = ˆ M div M ( V )( u − u d ) − ( ∂ M V + ∂ M V (cid:62) ) ∇ M u · ∇ M p + div M ( V )( ∇ M u · ∇ M p + up ) − f p div M ( V ) d s , (4.31) d Gd T t | t = ∂ V = ˆ M − ( u − u d ) ∇ u d · V − ∇ f · V p d s , (4.32)where ∂ M V denotes the tangential Jacobian of V and div M ( V ) : = ∂ M V : I the tangential diver-gence, which is defined as the trace of the tangential Jacobian; see [
6, p.495 ] .The implementation is analogous to the previous sections. We will only illustrate first orderderivatives here. We first define the geometry of the unit sphere, create a surface mesh anddefine a finite element space on the surface mesh: from netgen . c s g import ∗ from netgen . meshing import ∗ from n g s o l v e . i n t e r n a l import v i s o p t i o n s from n g s o l v e import ∗ g e o _ s u r f = CSGeometry ( ) sphere = Sphere ( Pnt ( 0 , 0 , 0 ) , 1 ) . bc ( " o u t e r " ) g e o _ s u r f . Add( sphere ) mesh_surf = Mesh ( g e o _ s u r f . GenerateMesh ( p e r f s t e p s e n d = MeshingStep . MESHSURFACE,o p t s t e p s 2 d = = mesh_surf . Curve (3) f e s _ s u r f = H1( mesh_surf , o r d e r = Next we define the transformed cost function and partial differential equation needed for settingup the Lagrangian (4.29). Here, we again make use of a symbolic object F to which we assignthe identity matrix. We define the tangential determinant ω and the matrix B defined in (4.28)as functions of the deformation gradient F t . X = C o e f f i c i e n t F u n c t i o n ( ( x , y , z ) ) func = C o e f f i c i e n t F u n c t i o n (X [ ] ∗ X [ ] ∗ X [ ] ) F = Id (3) tangDet = Det ( F ) ∗ Norm( Inv ( F ) . t r a n s ∗ s p e c i a l c f . normal (3) ) Bmat = ( Id (3) − / Norm( Inv ( F ) . t r a n s ∗ s p e c i a l c f . normal (3) ) ∗∗ ∗ OuterProduct ( Inv ( F ) .t r a n s ∗ s p e c i a l c f . normal (3) , Inv ( F ) . t r a n s ∗ s p e c i a l c f . normal (3) ) ) ∗ Inv ( F ) . t r a n s utomated Shape Optimisation in NGSolve d e f E q u a t i o n _ s u r f (u ,w) : r e t u r n ( ( Bmat ∗ grad ( u ) . Trace ( ) ) ∗ ( Bmat ∗ grad (w) . Trace ( ) ) + u ∗ w − func ∗ w) ∗ tangDet ∗ ds d e f C o s t _ s u r f ( u ) : r e t u r n u ∗∗ ∗ tangDet ∗ ds Now we can define the bilinear form and solve the state equation. Here, the right hand side ofthe equation is included in the bilinear form and the boundary value problem – although linear– is solved by Newton’s method (which terminates after only one iteration) for convenience. u _ s u r f , w_surf = f e s _ s u r f . TnT ( ) a = B i l i n e a r F o r m ( f e s _ s u r f ) a += E q u a t i o n _ s u r f ( u _ s u r f , w_surf ) g f u _ s u r f = G r i d F u n c t i o n ( f e s _ s u r f ) s o l v e r s . Newton ( a , g f u _ s u r f , p r i n t i n g = F a l s e )
Draw( g f u _ s u r f , mesh_surf , " g f u _ s u r f " )
The adjoint equation is solved as usual: l f c o s t _ s u r f = LinearForm ( f e s _ s u r f ) l f c o s t _ s u r f += C o s t _ s u r f ( g f u _ s u r f ) . D i f f ( g f u _ s u r f , w_surf ) l f c o s t _ s u r f . Assemble ( ) i n v a = a . mat . I n v e r s e ( f e s _ s u r f . F r e e D o f s ( ) , i n v e r s e = " s p a r s e c h o l e s k y " ) g f p _ s u r f = G r i d F u n c t i o n ( f e s _ s u r f ) g f p _ s u r f . vec . data = − i n v a . T ∗ l f c o s t _ s u r f . vec Draw( g f p _ s u r f , mesh_surf , " g f p _ s u r f " )
The shape derivative is obtained as in the case of PDEs posed on volumes by the evaluation of(4.30):
G _ s ur f = C o s t _ s u r f ( g f u _ s u r f ) + E q u a t i o n _ s u r f ( g f u _ s u r f , g f p _ s u r f )
VEC3d = VectorH1 ( mesh_surf , o r d e r = V3d = VEC3d . T e s t F u n c t i o n ( ) dJOmega_surf = LinearForm (VEC3d) dJOmega_surf += G _ s ur f . D i f f (X , V3d ) + G _ s ur f . D i f f ( F , Grad ( V3d ) . Trace ( ) )
In the previous sections we used the automatic differentiation capabilities of
NGSolve to alle-viate the shape differentiation procedure. However, so far we still had to include some knowl-edge about the problems at hand. So far, it was necessary to define the objective function orLagrangian G in the correct way, accounting for the correct transformation rules between per-turbed and unperturbed domain. In this section, we will show that also this step can be auto-mated since all necessary information is already included in the functional setting. The fullyautomated shape differentiation is incorporated by the command0 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
DiffShape(...) .In particular, in the fully automated setting it is enough to set up the cost function or Lagrangianfor the unperturbed setting. For a shape function of the type (3.1) we can define the shapederivative of the cost function in the following way:
G_f_0 = f ∗ dx dJOmega_f_0 = LinearForm (VEC) dJOmega_f_0 += G_f_0 . D i f f S h a p e (V)
Note that there is no term of the form
Det(F) showing up in line 186. Here, the transformationof the domain is taken care of automatically. It can be checked that this really gives the sameresult as dJOmega_f defined in lines 48–49. dJOmega_f . Assemble ( ) dJOmega_f_0 . Assemble ( ) d i f f e r e n c e V e c = dJOmega_f . vec . C r e a t e V e c t o r ( ) d i f f e r e n c e V e c . data = dJOmega_f . vec − dJOmega_f_0 . vec p r i n t ( " |dJOmega_f − dJOmega_f_0| = " , Norm( d i f f e r e n c e V e c ) ) The above code gives the output |dJOmega_f - dJOmega_f_0| = 1.571008573810619e-17 which confirms our claim. The same holds true for second order shape derivatives. The lines58–59 can be replaced by a repeated call of
DiffShape(...) : d2JOmega_f_0 = B i l i n e a r F o r m (VEC) d2JOmega_f_0 += G_f_0 . D i f f S h a p e (V) . D i f f S h a p e (W)
Again, it can be verified that d2JOmega_f_0 coincides with the previously defined quantity d2JOmega_f . Note that slightly different results may occur due to different integration rulesused. This can be cured by enforcing an integration rule of higher order for
G_f , i.e. by replacingthe symbol dx in the definition of G_f with dx(bonus_intorder=2) .In the more general setting of PDE-constrained shape optimisation, the procedure is verysimilar. Here the idea exploited in the implementation of the command
DiffShape(...) is tojust differentiate the general expression (4.4) with respect to the parameter t . The transforma-tions Φ t appearing in (4.4), which depend on the functional setting of the PDE, are identifiedautomatically from the finite element space from which the corresponding functions originate.The shape derivative of lines 82–83 can be obtained by the following code. d e f Cost_0 ( u ) : r e t u r n (u − ud ) ∗∗ ∗ dx d e f Equation_0 (u ,w) : r e t u r n ( grad ( u ) ∗ grad (w) − f 1 ∗ w) ∗ dx G_pde_0 = Cost_0 ( gfu ) + Equation_0 ( gfu , gfp ) dJOmega_pde_0 = LinearForm (VEC) dJOmega_pde_0 += G_pde_0 . D i f f S h a p e (V) utomated Shape Optimisation in NGSolve gfu and gfp represent the solutions to the state and adjoint equation, respectively, andmust have been computed previously. The bilinear form shapeHess11 used in Section 4.3 (seelines 121–122) can be obtained similarly: shapeHess11_0 = B i l i n e a r F o r m (VEC) shapeHess11_0 += G_pde_0 . D i f f S h a p e (W) . D i f f S h a p e (V)
The same holds true for boundary integrals
G_f_bnd_0 = f ∗ ds dJOmega_f_bnd_0 = LinearForm (VEC) dJOmega_f_bnd_0 += G_f_bnd_0 . D i f f S h a p e (V) and surface PDEs d e f C o s t _ s u r f _ 0 ( u ) : r e t u r n u ∗∗ ∗ ds d e f E q u a t i o n _ s u r f _ 0 (u ,w) : r e t u r n ( grad ( u ) . Trace ( ) ∗ grad (w) . Trace ( ) + u ∗ w − func ∗ w) ∗ ds G_surf_0 = C o s t _ s u r f _ 0 ( g f u _ s u r f ) + E q u a t i o n _ s u r f _ 0 ( g f u _ s u r f , g f p _ s u r f ) dJOmega_surf_0 = LinearForm (VEC3d) dJOmega_surf_0 += G_surf_0 . D i f f S h a p e ( V3d ) as well as their respective second order derivatives.
Remark 5.1.
We remark that the fully automated differentiation using
DiffShape(...) shouldbe seen to complement the semi-automated shape differentiation techniques introduced in Sec-tions 3 and 4 rather than to replace them. Using the semi-automated differentiation, the user hasthe possibility to, on the one hand, keep control over the involved terms, and on the other handalso to adjust the shape differentiation to their custom problems which may be non-standard. Asan example where the semi-automated differentiation may be beneficial compared to the fullyautomated differentiation we mention the case of time-dependent PDE constraints consideredin a space-time setting when a shape deformation is only desired in the spatial coordinates. Ofcourse, when one is interested in the shape derivative for a more standard problem, the fullyautomated way appears to be more convenient and less error prone.
We verify the expressions that we obtained in a semi-automatical or fully automatical way forthe first and second order shape derivatives by looking at the Taylor expansions of the perturbedshape functionals. We illustrate our findings in two examples in R . On the one hand, we con-sider a shape function as introduced in (3.1) with an additional boundary integral as in (3.6),henceforth denoted by (cid:74) ; on the other hand, we consider the PDE-constrained shape optimi-sation problem defined by (4.6), the reduced form of which will be denoted by (cid:74) ( Ω ) . More2 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel precisely, we consider (cid:74) ( Ω ) = ˆ Ω f ( x ) d x + ˆ ∂ Ω f ( x ) d s , (6.1) (cid:74) ( Ω ) = ˆ Ω | u Ω − u d | d x where u Ω solves (4.6b). (6.2)In the case of (cid:74) , we used the function f ( x , x ) = (cid:128) + (cid:198) x + x (cid:138) (cid:128) − (cid:198) x + x (cid:138) andfor (cid:74) , we used u d ( x , x ) = x ( − x ) x ( − x ) and f ( x , x ) = x ( − x ) + x ( − x ) forthe function f in the PDE constraint (4.6b).For the test of the first order shape derivatives D (cid:74) i ( Ω )( V ) we choose a fixed shape Ω and avector field V ∈ C ( R ) and observe the quantity δ ( (cid:74) i , t ) : = |(cid:74) i (( Id + t V )( Ω )) − (cid:74) i ( Ω ) − t D (cid:74) i ( Ω )( V ) | , (6.3)for t (cid:38)
0. Likewise, for the second order shape derivative, we consider the remainder δ ( (cid:74) i , t ) : = (cid:12)(cid:12)(cid:12)(cid:12) (cid:74) i (( Id + t V )( Ω )) − (cid:74) i ( Ω ) − t D (cid:74) i ( Ω )( V ) − t D (cid:74) i ( Ω )( V )( V ) (cid:12)(cid:12)(cid:12)(cid:12) (6.4)as t (cid:38)
0. By the definition of first and second order shape derivatives, it must hold that δ ( (cid:74) i , t ) = (cid:79) ( t ) and δ ( (cid:74) i , t ) = (cid:79) ( t ) as t (cid:38)
0. (6.5)This behavior can be observed in Figure 3(a) for (cid:74) and in Figure 3(b) for (cid:74) , where we used V ( x , x ) = ( x x e x , x x e x ) in both cases. -3 -2 -1 -18 -16 -14 -12 -10 -8 -6 -4 -3 -2 -1 -12 -10 -8 -6 -4 -2 (a) (b)Figure 3: Taylor test for functions (cid:74) and (cid:74) .The experiments for shape function (cid:74) was conducted on a mesh consisting of 13662 vertices,26946 elements and with polynomial order 2 (resulting in 54269 degrees of freedom), and theexperiment for (cid:74) with 95556 vertices and 190062 elements and polynomial degree 1 (95556degrees of freedom). We conducted these experiments for a number of different problems withdifferent vector fields V , in particular with different PDE constraints and boundary conditions,and obtained similar results in all instances provided a sufficiently fine mesh was used. utomated Shape Optimisation in NGSolve In this section we discuss how to use optimisation algorithms in conjunction with the automatedshape differentiation explained in the previous sections. The starting point of our discussion isa fixed initial shape Ω . Then we consider the mapping V (cid:55)→ g ( V ) : = (cid:74) (( Id + V )( Ω )) (6.6)defined on a suitable space of vector fields Θ ⊂ C ( D ) d . Since the mapping g is defined on anopen subset Θ of the Banach space C ( D ) d we can employ standard algorithms to minimise g over Θ . The only constraint we must impose is that Id + V remains invertible, which can bedifficult in practice. We observe that for V , W ∈ Θ we have ∂ g ( V )( W ) = D (cid:74) (( Id + V )( Ω ))( W ◦ ( Id + V ) − ) . (6.7) The gradient of ∂ g ( V ) in a Hilbert space H ⊂ C ( D ) d is defined by ∂ g ( V )( W ) = ( ∇ H g ( V ) , W ) H for all W ∈ H . (6.8)Typical choices for H are H = H ( D ) d , ( W , V ) H : = ˆ D ∂ W : ∂ V + V · W d x , (6.9) H = H ( D ) d , ( W , V ) H : = ˆ D (cid:34) ( W ) : (cid:34) ( V ) + V · W d x , (cid:34) ( V ) : = ( ∂ V + ∂ V (cid:62) ) , (6.10) H = H ( D ) d , ( W , V ) H : = ˆ D (cid:34) ( W ) : (cid:34) ( V ) + V · W + γ CR (cid:66) V · (cid:66) W d x , where γ CR > (cid:66) : = (cid:129) − ∂ x ∂ y ∂ y ∂ x (cid:139) . (6.12)The last choice corresponds to a penalised Cauchy-Riemann gradient and results in a gradientwhich is approximately conformal and hence preserves good mesh quality. We refer to [ ] fora detailed description. Let Ω be an initial shape and let H ⊂ C ( D ) d be a Hilbert space. Then a basic shape optimisationalgorithm reads as follows.4 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
Algorithm 1 gradient algorithm Input: domain Ω , n = N max > ε > γ ≥ Output: optimal shape Ω ∗ while n ≤ N max and |∇(cid:74) ( Ω n ) | > ε do if (cid:74) (( Id − α ∇(cid:74) ( Ω n ))( Ω n )) < (cid:74) ( Ω n ) − γα |∇(cid:74) ( Ω n ) | then Ω n + ← ( Id − α ∇(cid:74) ( Ω n ))( Ω n ) n ← n + increase α else reduce α end if end while We present and explain the numerical realisation of Algorithm 1 in
NGSolve for the caseof a PDE-constrained shape optimisation problem in two space dimensions. The simpler caseof an unconstrained shape optimisation problem or the case of three space dimensions can berealised by small modifications of the presented code.First of all, we mention that we realise shape modifications in
NGSolve by means of de-formation vector fields without actually modifying the coordinates of the underlying finite el-ement grid. Recall the vector-valued finite element space
VEC over a given mesh as intro-duced in code line 44. We define a vector-valued
GridFunction with the name gfset whichwill represent the current shape. We initialise it with some vector-valued coefficient function V ( x , x ) = ( x x , x x ) (cid:62) and obtain the deformed shape ( Id + V )( Ω ) by the command mesh.SetDeformation(gfset) : g f s e t = G r i d F u n c t i o n (VEC)
Draw( g f s e t , mesh , " g f s e t " )
S e t V i s u a l i z a t i o n ( deformation = True ) g f s e t . Se t ( ( X [ ] ∗ X [ ] ∗ X [ ] , X [ ] ∗ X [ ] ∗ X [ ] ) ) mesh . SetDeformation ( g f s e t ) Redraw ( )
Any operation involving the mesh such as integration or assembling of matrices is now carriedout for the deformed configuration. The deformation can be unset by the command mesh.UnsetDeformation() . Integrating the constant function over the mesh in the perturbedand unperturbed setting, p r i n t ( I n t e g r a t e ( 1 , mesh ) ) mesh . UnsetDeformation ( ) p r i n t ( I n t e g r a t e ( 1 , mesh ) ) gives the output respectively.In the course of the optimisation algorithm the state equation as well as the adjoint equationhave to be solved for every new shape. We define the following function, which computes thestate and adjoint state for a linear PDE constraint: utomated Shape Optimisation in NGSolve d e f solvePDE ( ) : a . Assemble ( ) L . Assemble ( ) dCostdu . Assemble ( ) i n v a = a . mat . I n v e r s e ( f e s . F r e e D o f s ( ) , i n v e r s e = " s p a r s e c h o l e s k y " ) gfu . vec . data = i n v a ∗ L . vec gfp . vec . data = − i n v a . T ∗ dCostdu . vec The shape derivative dJOmega for some problem at hand can be defined as illustrated in Sections4.1 and Section 5. Finally, we need to define the shape gradient, which is the solution to aboundary value problem of the form (6.8). We choose the bilinear form defined in (6.11) with γ CR = d e f eps ( u ) : r e t u r n 1 / ∗ ( grad ( u ) + grad ( u ) . t r a n s ) aX = B i l i n e a r F o r m (VEC)
W, V = VEC . TnT ( ) aX += I n n e r P r o d u c t ( eps (W) , eps (V) ) ∗ dx + I n n e r P r o d u c t (W, V) ∗ dx aX += ∗ ( grad (W) [ ] − grad (W) [ ] ) ∗ ( grad (V) [ ] − grad (V) [ ] ) ∗ dx aX += ∗ ( grad (W) [ ] + grad (W) [ ] ) ∗ ( grad (V) [ ] + grad (V) [ ] ) ∗ dx Now we can run Algorithm 1 for problem (4.6): alpha = a l p h a _ i n c r _ f a c t o r = gamma = − Nmax = e p s i l o n = − isConverged = F a l s e g f s e t . Se t ( ( 0 , 0 ) ) gfX = G r i d F u n c t i o n (VEC) gfsetTemp = G r i d F u n c t i o n (VEC) solvePDE ( )
Jnew = I n t e g r a t e ( Cost ( gfu ) , mesh )
J o l d = Jnew f o r k i n range (Nmax) : mesh . SetDeformation ( g f s e t ) aX . Assemble ( ) dJOmega_pde . Assemble ( ) invaX = aX . mat . I n v e r s e (VEC . F r e e D o f s ( ) , i n v e r s e = " s p a r s e c h o l e s k y " ) gfX . vec . data = invaX ∗ dJOmega_pde . vec currentNormGFX = Norm( gfX . vec ) while True : i f currentNormGFX < e p s i l o n : isConverged = True break P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel gfsetTemp . vec . data = g f s e t . vec − alpha ∗ gfX . vec mesh . SetDeformation ( gfsetTemp ) solvePDE ( ) Jnew = I n t e g r a t e ( Cost ( gfu ) , mesh ) mesh . UnsetDeformation ( ) i f Jnew < J o l d − gamma ∗ alpha ∗ currentNormGFX ∗∗ J o l d = Jnew g f s e t . vec . data = gfsetTemp . vec alpha ∗ = a l p h a _ i n c r _ f a c t o r break e l s e : alpha = alpha / Redraw ( b l o c k i n g = True )
Mesh movement and mesh optimisation
As an alternative to realizing the deformations via mesh.SetDeformation(...) , where the underlying mesh is not modified, one could alsojust move every mesh node in the direction of the given descent vector field by changing itscoordinates. This can be realised by invoking the following method: d e f moveNGmesh2D( d i s p l , mesh ) : f o r p i n mesh . ngmesh . P o i n t s ( ) : mip = mesh ( p [ ] , p [ ] ) v = d i s p l ( mip ) p [ ] += v [ ] p [ ] += v [ ] mesh . ngmesh . Update ( ) Here, the displacement vector field displ , which is of type
GridFunction , is evaluated foreach mesh node and, subsequently, the mesh nodes are updated. At the end of the procedure, themesh structure needs to be updated, see line 291. Note that the evaluation of
GridFunction srequires a mapped integration point mip of the mesh which is created in line 287.One advantage of this strategy is that a distorted mesh can easily be repaired by a call of themethod mesh.ngmesh.OptimizeMesh2d() followed by mesh.ngmesh.Update() . Figure 4shows a distorted mesh and the result of a call of mesh.ngmesh.OptimizeMesh2d() . The particular choice H = H ( D ) d and ( V , W ) H : = D (cid:74) ( Ω )( V )( W ) , (6.13)leads to Newton’s method. We refer to [
1, 8, 22, 24 ] where shape Newton methods were usedpreviously and to [
14, Chapter 2 ] and [
18, Chapter 5 ] for Newton’s method in an optimal controlsetting. This bilinear form is only positive semi-definite on H ( D ) d since D (cid:74) ( Ω )( V )( W ) = V , W with V = W = ∂ Ω . Moreover, from the structure theorem for second shape derivativesproved in [ ] we know that at a stationary point Ω , that is, D (cid:74) ( Ω )( V ) = V ∈ C ( D ) d ,we have D (cid:74) ( Ω )( V )( W ) = (cid:96) Ω ( V · n , W · n ) , (6.14) utomated Shape Optimisation in NGSolve mesh.ngmesh.OptimizeMesh2d() .where (cid:96) Ω : C ( ∂ Ω ) × C ( ∂ Ω ) → R is a bilinear function. Hence we also have D (cid:74) ( Ω )( V )( W ) = V , W such that V · n = W · n =
0. As a result the gradient ( ∇(cid:74) ( Ω ) , V ) H = D (cid:74) ( Ω )( V ) for all V ∈ H ( D ) d (6.15)according to (6.13) is not uniquely determined. To get around this difficulty, the shape Hessianis often regularised by an H term, i.e. (6.13) is replaced by D (cid:74) ( Ω )( V )( W ) + δ ˆ Ω ∂ V : ∂ W + V · W d x , (6.16)see, e.g. [ ] , which, however, impairs the convergence speed of Newton’s method. Alternative regularisation strategy.
Here, we propose the following strategy: We regularisethe shape Hessian only on the boundary ∂ Ω and only in tangential direction, i.e., we choose ( V , W ) H : = D (cid:74) ( Ω )( V )( W ) + δ ˆ ∂ Ω ( V · τ )( W · τ ) (6.17)with a regularisation parameter δ . To exclude the part of the kernel corresponding to interiordeformations, we solve the (regularised) Newton equation (6.15) only on the boundary ∂ Ω .This is realised by setting Dirichlet boundary conditions for all degrees of freedom except thoseon the boundary. VEC2 = VectorH1 ( mesh , o r d e r =
1, d i r i c h l e t = " c i r c l e " ) aX = B i l i n e a r F o r m (VEC) aX += G_f_0 . D i f f S h a p e (W) . D i f f S h a p e (V) aX += ∗ I n n e r P r o d u c t (W, s p e c i a l c f . t a n g e n t i a l (2) ) ∗ I n n e r P r o d u c t (V , s p e c i a l c f .t a n g e n t i a l (2) ) ∗ ds aX . Assemble ( ) invAX = aX . mat . I n v e r s e (~VEC2 . F r e e D o f s ( ) , i n v e r s e = " umfpack " ) gfX_bnd = G r i d F u n c t i o n (VEC) gfX_bnd . vec . data = invAX ∗ dJOmega_f_0 . vec P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
As a result, we get a shape gradient ˜ ∇(cid:74) ( Ω ) which is nonzero only on the boundary. We extendthis vector field to the interior by solving an additional boundary value problem (of linearisedelasticity type), where we use the deformation given by ˜ ∇(cid:74) ( Ω ) as Dirichlet boundary condi-tions. d e f g e t E x t e n s i o n ( gfX_bnd , f r e e d o f s , g f X _ e x t ) : u , v = VEC . TnT ( ) aX_ext = B i l i n e a r F o r m (VEC) aX_ext += I n n e r P r o d u c t ( grad ( u ) + grad ( u ) . t r a n s , grad ( v ) ) ∗ dx + I n n e r P r o d u c t (u , v ) ∗ dx g f X _ e x t . S et ( gfX_bnd ) aX_ext . Assemble ( ) r = gfX_bnd . vec . C r e a t e V e c t o r ( ) r . data = ( − ∗ aX_ext . mat ∗ g f X _ e x t . vec g f X _ e x t . vec . data += aX_ext . mat . I n v e r s e ( f r e e d o f s = f r e e d o f s ) ∗ r g e t E x t e n s i o n ( gfX_bnd , VEC2 . F r e e D o f s ( ) , gfX ) g f s e t . Se t ( ( 0 , 0 ) ) g f s e t . vec . data = g f s e t . vec − ∗ gfX . vec The Newton algorithm reads as follows.
Algorithm 2
Newton algorithm Input: domain Ω , n = N max > ε > Output: optimal shape Ω ∗ while n ≤ N max and |∇(cid:74) ( Ω n ) | > ε do solve (6.15) to get ∇(cid:74) ( Ω n ) Ω n + ← ( Id − ∇(cid:74) ( Ω n ))( Ω n ) n ← n + end while6.2.4 Newton’s method for PDE-constrained problems We consider the PDE-constrained model problem of Section 4.1 which is subject to the Poissonequation. The unregularised Newton system reads D (cid:74) ( Ω )( V )( W ) = − D (cid:74) ( Ω )( V ) for all V ∈ H ( Ω ) . (6.18)In Subsection 4.3 we discussed how the second order derivative can be evaluated. Recalling that d (cid:74) ( Ω )( V ) = ∂ s G V ,0 (
0, 0, u , p ) we see that (4.23) and (4.19) lead to ∂ s ∂ t G V , W (
0, 0, u , p ) ∂ u ∂ t G V , W (
0, 0, u , p ) ∂ p ∂ t G V , W (
0, 0, u , p ) ∂ s ∂ u G V , W (
0, 0, u , p ) ∂ u G V , W (
0, 0, u , p ) ∂ p ∂ u G V , W (
0, 0, u , p ) ∂ s ∂ p G V , W (
0, 0, u , p ) ∂ u ∂ p G V , W (
0, 0, u , p ) ˜ V ∂ s u ∂ s p = − ∂ t G V , W (
0, 0, u , p ) .(6.19)The component ˜ V then represents the direction which we use for the shape Newton optimisationstep. The matrix in (6.19) can be realised in NGSolve by using a combined finite element space X3 consisting of three components as follows: utomated Shape Optimisation in NGSolve X3 = FESpace ( [ VEC , f e s , f e s ] ) PHI , u1 , p1 = X3 . T r i a l F u n c t i o n ( )
PSI , uTest1 , pTest1 = X3 . T e s t F u n c t i o n ( ) shapeHessLag3 = B i l i n e a r F o r m ( X3 ) shapeHessLag3 += G_pde_0 . D i f f S h a p e ( PHI ) . D i f f S h a p e ( PSI ) shapeHessLag3 += G_pde_0 . D i f f S h a p e ( PSI ) . D i f f ( gfu , u1 ) shapeHessLag3 += G_pde_0 . D i f f S h a p e ( PSI ) . D i f f ( gfp , p1 ) shapeHessLag3 += G_pde_0 . D i f f ( gfu , uTest1 ) . D i f f S h a p e ( PHI ) shapeHessLag3 += ( G_pde_0 . D i f f ( gfu , uTest1 ) ) . D i f f ( gfu , u1 ) shapeHessLag3 += ( G_pde_0 . D i f f ( gfu , uTest1 ) ) . D i f f ( gfp , p1 ) shapeHessLag3 += G_pde_0 . D i f f ( gfp , pTest1 ) . D i f f S h a p e ( PHI ) shapeHessLag3 += ( G_pde_0 . D i f f ( gfp , pTest1 ) ) . D i f f ( gfu , u1 ) The right hand side of (6.19) can be defined as follows: shapeGradLag3 = LinearForm ( X3 ) shapeGradLag3 += ( − ∗ G_pde_0 . D i f f S h a p e ( PSI )
Recall that the system (6.15) has a nontrivial kernel as discussed in Section 6.2.3. This problemcan be circumvented by proceeding like in the unconstrained case. We add a regularisation onlyon the boundary, d e l t a = shapeHessLag3 += d e l t a ∗ I n n e r P r o d u c t ( PHI , s p e c i a l c f . t a n g e n t i a l (2) ) ∗ I n n e r P r o d u c t( PSI , s p e c i a l c f . t a n g e n t i a l (2) ) ∗ ds and exclude the interior degrees of freedom in the first row and column of the 3 × VEC2 = VectorH1 ( mesh , o r d e r = = " . ∗ " ) freeDofsCombined = B i t A r r a y (VEC2 . ndof + ∗ f e s . ndof ) f o r i i n range (VEC2 . ndof ) : freeDofsCombined [ i ] = not VEC2 . F r e e D o f s ( ) [ i ] f o r i i n range ( f e s . ndof ) : freeDofsCombined [ VEC2 . ndof + i ] = f e s . F r e e D o f s ( ) [ i ] freeDofsCombined [ VEC2 . ndof + f e s . ndof + i ] = f e s . F r e e D o f s ( ) [ i ] and solving the regularised system using these free dofs: gfCombined3 = G r i d F u n c t i o n ( X3 ) shapeHessLag3 . Assemble ( ) shapeGradLag3 . Assemble ( ) gfCombined3 . vec . data = shapeHessLag3 . mat . I n v e r s e ( f r e e d o f s = freeDofsCombined ,i n v e r s e = " umfpack " ) ∗ shapeGradLag3 . vec The newton direction is then given as the first of the three components of the obtained solution.
V t i l d e _ b n d = G r i d F u n c t i o n (VEC)
V t i l d e = G r i d F u n c t i o n (VEC)
V t i l d e _ b n d . vec . data = gfCombined3 . components [ ] . vec P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel g e t E x t e n s i o n ( V t i l d e_ b n d , VEC2 . F r e e D o f s ( ) , V t i l d e ) g f s e t . vec . data = g f s e t . vec + ∗ V t i l d e . vec
In this section we will apply the automated shape differentiation and all numerical algorithmsintroduced in the preceding sections in numerical examples.
In this section, we revisit problem (3.1) introduced in Section 3, i.e. the problem of finding ashape Ω such that the cost function (cid:74) ( Ω ) = ´ Ω f ( x ) d x is minimised. First order methods
We illustrate our first order methods in a problem which was also con-sidered in [ ] and reproduce the results obtained there. We choose the function f ( x , x ) = (cid:128)(cid:113) ( x − a ) + b x − (cid:138) (cid:128)(cid:113) ( x + a ) + b x − (cid:138) · (cid:128)(cid:113) b x + ( x − a ) − (cid:138) (cid:128)(cid:113) b x + ( x + a ) − (cid:138) − ε (6.20)with a = , b = ε = { ( x , x ) ∈ R : f ( x , x ) < } which is depicted in Figure 5 (right). We start our optimisation algorithm withthe unit disk, Ω = B ( ) as an initial design. Note that the optimal design cannot be reachedby means of shape optimisation using boundary perturbations. However, we expect the outercurve of the optimal shape to be reached very closely.We apply Algorithm 1 with the shape gradient ∇(cid:74) associated to the H inner product (6.9),to the bilinear form of linearised elasticity (6.10) and including the additional Cauchy-Riemannterm (6.11). We chose the algorithmic parameters γ = e − ε = e −
7, a mesh consisting of2522 vertices and 4886 elements and a globally continuous vector-valued finite element space
VEC of order 3. The results can be seen in Figures 6, 7 and 8, respectively.
Second order method
Since Newton’s method converges quadratically only in a neighbor-hood of the optimal solution, we choose a simpler optimal design here. We choose f ( x , x ) = x a + y b − a and b . We choose a = b = / a and again start the optimisation with the unit disk as initial shape. Figure 9 showsthe initial and optimised design after only six iterations of Algorithm 2 with ( · , · ) H chosen asin (6.17) with δ = δ =
100 and (6.16) with δ = δ was chosen empirically to get convergence as fast as possible. The experimentswere conducted on a finite element mesh consisting of 2522 nodes and 4886 triangular elementswith a finite element space VEC of order 3, with the algorithmic parameter ε = − . utomated Shape Optimisation in NGSolve Ω and optimal domain Ω ∗ for problem (3.1) with f chosen accordingto (6.20).Figure 6: Results of problem (3.1) with f as in (6.20) and the shape gradient associated to the H inner product (6.9).2 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
Figure 7: Results of problem (3.1) with f as in (6.20) and the shape gradient associated to theelasticity bilinear form (6.10).Figure 8: Results of problem (3.1) with f as in (6.20) and the shape gradient associated to theelasticity bilinear form with Cauchy-Riemann term (6.11).Figure 9: Numerical results for problem (6.17) with f as in (6.21) using second order method.Left: Initial design. Center: Optimised design after six iterations using (6.15) / (6.17). Right: Ob-jective value (cid:74) and norm of shape gradient (cid:107)∇(cid:74) ( Ω ) (cid:107) in the course of second order optimisationusing (6.16) with δ = δ = utomated Shape Optimisation in NGSolve Iterations -9 -8 -7 -6 -5 -4 J Iterations -5 -4 -3 -2 -1 (a) (b)Figure 10: Convergence behaviour for shape optimisation problem (4.6) with proposed regular-isation strategies (6.17) and (6.16) as well as first order method with constant step size α = (cid:74) . (b) Behaviour of norm of shape gradient (cid:107)∇(cid:74) ( Ω ) (cid:107) . In this section, we revisit the model problem introduced in Section 4.1 with f ( x , x ) = x ( − x ) + x ( − x ) and u d ( x , x ) = x ( − x ) x ( − x ) . Note that the data is chosen in such away that, for Ω ∗ = (
0, 1 ) it holds (cid:74) ( Ω ∗ ) = Ω ∗ is a global minimiser of (cid:74) . We showresults obtained by first and second order shape optimisation methods exploiting automateddifferentiation.We ran the optimisation algorithm in three versions. On the one hand, we applied a first ordermethod with constant step size α =
1. On the other hand, we applied two second order methodswith the two different regularisation strategies for the shape Hessian in (6.15) introduced in(6.16) and (6.17). We chose the regularisation parameters δ empirically such that the methodperforms as well as possible. In the case of (6.16) we chose δ = δ =
1. The experiments were conducted on a finite element mesh consisting of 4886 elementswith 2522 vertices and polynomial degree 1. In Figure 10, we can observe the decrease of theobjective function as well as of the norm of the shape gradient over 200 iterations for thesethree algorithmic settings.Figure 10 shows the initial design as well as the design after 200 iterations of the secondorder method with regularisation strategy (6.17). Note that the improved design is very closeto the global solution Ω ∗ = (
0, 1 ) . The initial design was chosen as the disk of radius centeredat the point (cid:0) , (cid:1) (cid:62) . The objective value was reduced from 5.297 · − to 1.0317 · − . Here, we illustrate the applicability of the automated shape differentiation and optimisation inthe more realistic and more complicated setting of nonlinear elasticity in two space dimensionsusing a Saint Venant–Kirchhoff material with Young’s modulus E = ν = P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
Figure 11: Shape optimisation for problem (4.6). Left: Initial design. Right: Improved designafter 200 iterations of second order algorithm with regularisation as proposed in (6.17). Objec-tive value was reduced from 5.297 · − to 1.0317 · − . Color shows solution of constrainingPDE (4.6b).left parts of the boundary, Γ l = { } × ( ) and Γ l = { } × (
0, 0.12 ) , respectively, and issubject to a surface force g N = ( − ) (cid:62) on Γ r = { } × ( ) . The initial geometry with3 holes is depicted in Figure 12 (a). Let Γ l : = Γ l ∪ Γ l and H Γ l ( Ω ) the subspace of H ( Ω ) withvanishing trace on Γ l . The displacement u ∈ H Γ l ( Ω ) under the surface force g N is given as thesolution to the boundary value problem ˆ Ω S ( u ) : ∇ v d x = ˆ Γ r g N · v d s (6.22)for all v ∈ H Γ l ( Ω ) . Here, S ( u ) denotes the Saint Venant–Kirchhoff stress tensor S ( u ) = ( I + ∇ u ) (cid:149) λ Tr (cid:129) ( C ( u ) − I ) (cid:139) I + µ ( C ( u ) − I ) (cid:152) , (6.23)where C ( u ) = ( I + ∇ u ) (cid:62) ( I + ∇ u ) and I is the identity matrix, see also [
2, Sec. 8 ] , and λ and µ denote the Lamé constants, λ = E ν ( + ν )( − ν ) , µ = E ( + ν ) . (6.24)We minimise the functional J ( Ω , u ) = ˆ Ω S ( u ) : ∇ u d x + α ˆ Ω x (6.25)with α = [
2, Sec.8 ] . Nevertheless, application of the automated shape differentiation and optimisation yieldsa significant improvement of the initial design. The highly nonlinear PDE constraint (6.22) issolved by Newton’s method. In order to have good starting values, a load stepping strategy is utomated Shape Optimisation in NGSolve alpha = 0.1 (as an initial value), alpha_incr_factor =1 (i.e. no increase), gamma = 1e-4 and epsilon = 1e-7 . Moreover, we used (6.9) with anadditional Cauchy-Riemann term as in (6.11) with weight γ CR =
10. The objective value wasreduced from 3.125 to 2.635 (volume term from 1.290 to 1.096) in 15 iterations of Algorithm1. The results were obtained on a mesh consisting of 10614 elements and 5540 vertices usingpiecewise linear, globally continuous finite elements.
In this section, we consider the problem of finding the optimal shape of a scattering object. Moreprecisely, we consider the minimisation of the functional ˆ Γ r uu d s (6.26)subject to the Helmholtz equation with impedance boundary conditions on the outer boundary,Find u ∈ H ( Ω , C ) : ˆ Ω [ ∇ u · ∇ ¯ w − ω u ¯ w (cid:3) d x − i ω ˆ Γ u ¯ w d s = ˆ Ω f ¯ w (6.27)for all w ∈ H ( Ω , C ) . Here, w denotes the complex conjugate of a complex-valued function w , ω denotes the wave number, i denotes the complex unit and the function f on the right handside is chosen as f ( x , x ) = · e − (( x − ) +( x − ) ) , (6.28)see Figure 13(a). Furthermore Ω = B (( ) (cid:62) , 1 ) \ B (( ) (cid:62) , 0.15 ) denotes the domain ofinterest, Γ = { ( x , x ) : x + x = } the outer boundary and Γ r = { ( x , x ) : x + x = x ≥ } P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel (a) (b) (c)Figure 13: (a) Geometry with right hand side. (b) Initial shape of scatter. (c) Optimised shapeof scatter. (a) (b)Figure 14: (a) Absolute value of state u for initial configuration. (b) Absolute value of state u for optimised configuration.the right half of the outer boundary. Here, only the inner boundary ∂ Ω \ Γ is subject to the shapeoptimisation. Thus, the aim of this model problem is to find a shape of the scattering object suchthat the waves are reflected away from Γ r .Figure 13 (b) and (c) show the initial and final shape of the scattering object, respectively.Figure 14 shows the norm of the state for the initial configuration (circular shape of scatteringobject) and for the optimised configuration. The objective value was reduced from 3.44 · − to3.31 · − . The forward simulations were performed using piecewise linear finite elements ona triangular grid consisting of with 34803 degrees of freedom. The optimisation stopped after12 iterations. In this section, we consider the setting of three-dimensional nonlinear magnetostatics in H ( curl , D ) as it appears in the simulation of electrical machines. Let D ⊂ R denote the computational do-main, which consists of ferromagnetic material, air regions and permanent magnets. Our aim isto minimise the functional ˆ Ω g | curl u · n − B nd | d x , (6.29)where Ω g denotes the air gap region of the machine, n denotes an extension of the normalvector to the interior of Ω g , B nd : Ω g → R is a given smooth function and u ∈ H ( curl , D ) is the utomated Shape Optimisation in NGSolve ˆ D ν Ω ( | curl u | ) curl u · curl w + δ u · w d x = ˆ Ω m M · curl w d x (6.30)for all w ∈ H ( curl , D ) . Here, Ω ⊂ D denotes the union of the ferromagnetic parts of the electricalmachine, Ω m denotes the permanent magnets subdomain and ν Ω = χ Ω ( x ) ˆ ν ( | curl u | ) + χ D \ Ω ( x ) ν (6.31)denotes the magnetic reluctivity, which is a nonlinear function ˆ ν inside the ferromagnetic regionsand equal to a constant ν elsewhere. Further, δ > M : D → R denotes the magnetisation in the permanent magnets. The nonlinear function ˆ ν satisfies a Lipschitz condition and a strong monotonicity condition such that problem (6.30) iswell-posed. We refer the reader to [
11, Sec. 6 ] for a more detailed description of the problemand [ ] for a 2D version of the same problem.As mentioned in Section 4, the transformation Φ t used in (4.4) depends on the differentialoperator. For the curl -operator, we have Φ t ( u ) = ∂ T −(cid:62) t ( u ◦ T − t ) and ( curl ( Φ t ( u ))) ◦ T t = ( ∂ T t ) ∂ T t curl ( u ) , (6.32)see e.g. [
20, Section 3.9 ] . Thus, the variational equation (6.30) can be defined as follows. from math import p i nu0 = / (4 ∗ p i ) d e l t a = F = Id (3) c1 = / Det ( F ) ∗ F c2 = Inv ( F ) . t r a n s d e f E q u a t i o n I r o n (u ,w) : r e t u r n ( nuIron (Norm( c1 ∗ c u r l ( u ) ) ) ∗ ( c1 ∗ c u r l ( u ) ) ∗ ( c1 ∗ c u r l (w) ) + d e l t a ∗ ( c2 ∗ u ) ∗ (c2 ∗ w) ) ∗ Det ( F ) ∗ dx ( " i r o n " ) d e f E q u a t i o n A i r (u ,w) : r e t u r n ( nu0 ∗ ( c1 ∗ c u r l ( u ) ) ∗ ( c1 ∗ c u r l (w) ) + d e l t a ∗ ( c2 ∗ u ) ∗ ( c2 ∗ w) ) ∗ Det ( F ) ∗ dx ( " a i r |a i r g a p " ) d e f EquationMagnets (u ,w) : r e t u r n ( nu0 ∗ ( c1 ∗ c u r l ( u ) ) ∗ ( c1 ∗ c u r l (w) ) + d e l t a ∗ ( c2 ∗ u ) ∗ ( c2 ∗ w) − magn ∗ ( c1 ∗ c u r l (w) ) ) ∗ Det ( F ) ∗ dx ( " magnets " ) d e f Equation (u ,w) : r e t u r n E q u a t i o n I r o n (u ,w) + E q u a t i o n A i r (u ,w) + EquationMagnets (u ,w)
Here, the computational domain consists of a subdomain representing the ferromagnetic partof the machine ( “iron” ) and a subdomain comprising the permanent magnets ( “magnets” );the union of all air subdomains, including the air gap between rotating and fixed part of themachine, is given by “air|airgap” . Moreover, nuIron denotes the nonlinear reluctivity func-tion ˆ ν and magn contains the magnetization direction of the permanent magnets. Likewise, theobjective function can be defined as follows,8 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel (a) (b)Figure 15: (a) Initial design of electrical machine. (b) Optimised design. -2 80-1010 -4 curl(u)*n in air gap (initial design) angle z-component
020 -1-20 -2 80-1010 -4 curl(u)*n in air gap (optimized design) angle z-component
020 -1-20 -2 80-1010 -4 desired curve in air gap angle z-component
020 -1-20 (a) (b) (c)Figure 16: Improvement of curl u · n as a function of z and the angle ϕ for a fixed radius r compared to desired curve B nd . (a) curl u · n for initial configuration. (b) curl u · n for optimisedconfiguration after 10 iterations. (c) Desired curve B nd in polar coordinates as function of z andthe angle ϕ for a fixed radius r . d e f Cost ( u ) : r e t u r n ( I n n e r P r o d u c t ( c1 ∗ c u r l ( u ) ,n2D) − Bd) ∗ Det ( F ) ∗ dx ( " a i r g a p " ) where n2D and Bd represent the extension of the normal vector to the interior of the air gapand the desired curve, respectively. For the definition of all quantities, we refer the reader to thecode which was submitted together with this manuscript. The shape differentiation as well asthe optimisation loop now works in the same way as in the previous examples. Figure 15 showsthe initial design of the motor as well as the optimised design obtained after 11 iterations ofAlgorithm 1 with γ =
0. The experiment was conducted using a tetrahedral finite element meshconsisting of 13440 vertices, 57806 elements and Nédélec elements of order 2 (resulting in atotal of 323808 degrees of freedom). The objective value was reduced from 2.5944 · − to4.565 · − in the course of the first order optimisation algorithm after 11 iterations. It can beseen from Figure 16 that the difference between the quantity curl ( u ) · n and the desired curve B nd inside the air gap decreases significantly. utomated Shape Optimisation in NGSolve Iterations -9 -8 -7 -6 -5 -4 -3 J -7 -6 -5 -4 -3 -2 History of surface Laplacian (a) (b)Figure 17: (a) Initial geometry for shape optimisation with respect to surface PDE (4.24)–(4.25).(b) History of objective value and norm of shape gradient using a first order algorithm with linesearch.
Finally, we also show the application of a shape optimisation algorithm to a problem constrainedby a surface PDE. We solve problem (4.24)–(4.25) with u d = f ( x , x , x ) = x x x andinitial shape M = S the unit sphere in R . We applied a first order algorithm with a line search.Figure 17 shows the initial geometry as well as the decrease of the objective function and ofthe norm of the shape gradient. The objective value was reduced from 7.08 · − to 9.88 · − .Figure 18 shows the final design which was obtained after 575 iterations from two differentperspectives. The experiment was conducted using a surface mesh with 332 vertices and 660faces and polynomial degree 3 (resulting in 2972 degrees of freedom).(a) (b)Figure 18: (a) Final design after 575 iterations. (b) Different view of (a).0 P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel
Conclusion
We showed how to obtain first and second order shape derivatives for unconstrained and PDE-constrained shape optimization problems in a semi-automatic and fully automatic way in thefinite element software package
NGSolve . We verified the proposed method numerically byTaylor tests and by showing its successful application to several shape optimisation problems.We believe that this intuitive approach can help research scientists working in the field of shapeoptimisation to further improve numerical methods on the one hand, and product engineersworking with
NGSolve to design devices in an optimal fashion on the other hand.
Acknowledgements
Michael Neunteufel has been funded by the Austrian Science Fund (FWF) project W1245.
Replication of results
The python scripts which were used for the results presented in this paper are available withthe arxiv submission. All computations were performed using NGSolve version V6.2.2004.
References [ ] G. Allaire, E. Cancès, and J. L. Vié. Second-order shape derivatives along normal trajec-tories, governed by Hamilton-Jacobi equations.
Structural and Multidisciplinary Optimiza-tion , 54(5):1245–1266, 2016. [ ] G. Allaire, F. Jouve, J., and A.-M. Toader. Structural optimization using sensitivity analysisand a level-set method.
Journal of Computational Physics , 194(1):363 – 393, 2004. [ ] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells. Unified form language.
ACM Transactions on Mathematical Software , 40(2):1–37, 2014. [ ] M. Berggren. A unified discrete-continuous sensitivity analysis method for shape optimiza-tion. In
Applied and numerical partial differential equations , volume 15 of
Comput. MethodsAppl. Sci. , pages 25–39. Springer, New York, 2010. [ ] M. C. Delfour and J.-P. Zolésio.
Shapes and geometries , volume 22 of
Advances in Designand Control . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA,second edition, 2011. Metrics, analysis, differential calculus, and optimization. [ ] M. C. Delfour and J. P. Zolésio.
Shapes and geometries . Society for Industrial and AppliedMathematics, 2011. [ ] J. S. Dokken, S. K. Mitusch, and S. W. Funke. Automatic shape derivatives for transientPDEs in FEniCS and Firedrake. arXiv e-prints , page arXiv:2001.10058, 2020. utomated Shape Optimisation in NGSolve [ ] K. Eppler, H. Harbrecht, and R. Schneider. On convergence in elliptic shape optimization.
SIAM Journal on Control and Optimization , 46(1):61–83, 2007. [ ] F. Feppon, G. Allaire, F. Bordeu, J. Cortial, and C. Dapogny. Shape optimization of a cou-pled thermal fluid-structure problem in a level set mesh evolution framework.
SeMA ,76(3):413–458, 2019. [ ] P. Gangl, U. Langer, A. Laurain, H. Meftahi, and K. Sturm. Shape optimization of an elec-tric motor subject to nonlinear magnetostatics.
SIAM Journal on Scientific Computing ,37(6):B1002–B1025, 2015. [ ] P. Gangl and K. Sturm. Asymptotic analysis and topological derivative for 3D quasi-linearmagnetostatics, 2019. arXiv:1908.10775. [ ] D. A. Ham, L. Mitchell, A. Paganini, and . Wechsung, F. Automated shape differentiationin the unified form language.
Structural and Multidisciplinary Optimization , 60(5):1813–1820, 2019. [ ] M. Hintermüller and A. Laurain. Electrical impedance tomography: from topology toshape.
Control and Cybernetics , 37(4):913–933, 2008. [ ] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich.
Optimization with PDE constraints .Springer, New York, 2009. [ ] R. Hiptmair, A. Paganini, and S. Sargheini. Comparison of approximate shape gradients.
BIT , 55(2):459–485, 2015. [ ] D. Hömberg and J. Sokolowski. Optimal shape design of inductor coils for surface hard-ening.
SIAM Journal on Control and Optimization , 42(3):1087–1117, 2003. [ ] J. A. Iglesias, K. Sturm, and F. Wechsung. Two-dimensional shape optimization with nearlyconformal transformations.
SIAM Journal on Scientific Computing , 40(6):A3807–A3830,2018. [ ] K. Ito and K. Kunisch.
Lagrange multiplier approach to variational problems and applica-tions . Society for Industrial and Applied Mathematics, Philadelphia, PA, 2008. [ ] A. Laurain. A level set-based structural optimization code using FEniCS.
Structural andMultidisciplinary Optimization , 58(3):1311–1334, 2018. [ ] P. Monk.
Finite element methods for Maxwell’s equations . Numerical Mathematics andScientific Computation. Clarendon Press, 2003. [ ] A. Novruzi and M. Pierre. Structure of shape derivatives.
Journal of Evolution Equations ,2(3):365–382, 2002. [ ] A. Novruzi and J. R. Roche. Newton’s method in shape optimisation: A three-dimensionalcase.
Bit Numerical Mathematics , 40(1):102–120, 2000.2
P. Gangl, J. Schöberl, K. Sturm and M. Neunteufel [ ] A. Paganini, S. Sargheini, R. Hiptmair, and Ch. Hafner. Shape optimization of microlenses.
Optics Express , 23(10):13099, 2015. [ ] A. Paganini and K. Sturm. Weakly normal basis vector fields in RKHS with an applicationto shape Newton methods.
SIAM Journal on Numerical Analysis , 57(1):1–26, 2019. [ ] Anton Schiela and Julian Ortiz. Second order directional shape derivatives, March 2017. [ ] S. Schmidt. A two stage CVT / eikonal convection mesh deformation approach for largenodal deformations. arXiv e-prints , page arXiv:1411.7663, 2014. [ ] S. Schmidt. Weak and strong form shape Hessians and their automatic generation.
SIAMJournal on Scientific Computing , 40(2):C210–C233, 2018. [ ] S. Schmidt, C. Ilic, V. Schulz, and N. Gauger. Three-dimensional large-scale aerodynamicshape optimization based on shape calculus.
AIAA Journal , 51(11):2615–2627, 2013. [ ] S. Schmidt, C. Ilic, V. Schulz, and N. R. Gauger. Airfoil design for compressible inviscidflow based on shape calculus.
Optimization and Engineering , 12(3):349–369, 2011. [ ] J. Schöberl. C ++
11 implementation of finite elements in NGSolve. Technical Report 30,Institute for Analysis and Scientific Computing, Vienna University of Technology, 2014. [ ] V. H. Schulz. A riemannian view on shape optimization.
Foundations of ComputationalMathematics , 14(3):483–501, 2014. [ ] J. Sokołowski and J.-P. Zolésio.
Introduction to shape optimization , volume 16 of
SpringerSeries in Computational Mathematics . Springer-Verlag, Berlin, 1992. Shape sensitivityanalysis. [ ] K. Sturm. Minimax Lagrangian approach to the differentiability of nonlinear PDE con-strained shape functions without saddle point assumption.
SIAM Journal on Control andOptimization , 53(4):2017–2039, 2015. [ ] K. Sturm. Shape differentiability under non-linear PDE constraints. In
New trendsin shape optimization , volume 166 of
Internat. Ser. Numer. Math. , pages 271–300.Birkhäuser //