Fully and semi-automated shape differentiation in NGSolve

In this paper, we present a framework for automated shape differentiation in the finite element software NGSolve. Our approach combines the mathematical Lagrangian approach for differentiating PDE-constrained shape functions with the automated differentiation capabilities of NGSolve. The user can decide which degree of automatisation is required, thus allowing 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 unconstrained model problems as well as for more realistic problems that are constrained by different types of partial differential equations. We consider linear as well as nonlinear problems and also problems which are posed on surfaces. In numerical experiments, we verify the accuracy of the computed derivatives via a Taylor test. Finally, we present first- and second-order shape optimisation algorithms and illustrate them for several numerical optimisation examples ranging from nonlinear elasticity to Maxwell’s equations. Supplementary Information The online version contains supplementary material available at 10.1007/s00158-020-02742-w.


Introduction
Numerical simulation and shape optimisation tools to solve the problems have become an integral part in the design process of many products.Starting out from an initial design, nonparametric shape optimisation techniques based on first and second order shape derivatives can assist in finding shapes of a product which are optimal with respect to a given objective function.Examples include the optimal design of aircrafts [28,29], optimal inductor design [16], optimisation of microlenses [23], the optimal design of electric motors [10], applications from mechanical engineering [2,19], multiphysics problems [9] or electrical impedance tomography (EIT) in medical sciences to name only a few [13].
Shape optimisation algorithms are based on the concept of shape derivatives.Let ⊂ (R d ) a set of admissible shapes and : → R a shape function.Given an admissible shape Ω ∈ and a sufficiently smooth vector field V , we define the perturbed domain Ω t := (Id + t V )(Ω) for a small perturbation parameter t > 0. The shape derivative is defined as In most practically relevant applications, the objective functional depends on the shape of a (sub-)domain via the solution of a partial differential equation (PDE).Thus, one is facing a problem of PDE-constrained shape optimisation of the form min (Ω,u)∈ ×X J(Ω, u) s.t.(Ω, u) ∈ × X : e(Ω; u, v) = 0 for all v ∈ X . (1.2) Here, the second line represents the constraining boundary value problem posed on a Hilbert space X , which we assume to be uniquely solvable for all admissible Ω ∈ .Denoting the unique solution for a given Ω ∈ by u Ω , we introduce the notation for the reduced functional (Ω) := 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 [33] for an overview of different approaches.In the following we focus on computing the so-called volume form of the shape derivative which in a finite element context is known to give a better approximation compared to the boundary form; see [4,15].The convergence of shape optimisation algorithms can be speeded up by using second order shape derivatives.Given two sufficiently smooth vector fields V, W and an admissible shape Ω ∈ , let Ω s,t := (Id + sV + tW )(Ω) the perturbed domain.Then, the second order shape derivative is defined as . ( 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 error prone, several authors have employed automatic differentiation (AD) tools, see e.g.[27] and [12] for two approaches based on the Unified Form Language (UFL) [3].In [12], the authors present a fully automated shape differentiation software which uses the transformation properties on the finite element level.In [27] (see also the earlier work [26]) the automated derivatives are computed using UFL.The strategy of [12] and [27] differ in that, for the latter, the software computes an unsymmetric shape Hessian since it involves the term D (Ω)(∂ V W ).
Optionally the software allows to make the shape Hessian symmetric by requiring ∂ V W = 0. Let us also mention [7] where automated shape derivatives for transient PDEs in FEniCS and Firedrake are presented.
In this paper we present an alternative framework for AD of PDE constrained problems of type (1.2).There exist several approaches for the rigorous derivation of the shape derivative of PDE-constrained shape functionals, see [34] for an overview.The main idea, however, is always similar.After transforming the perturbed setting back to the original domain, shape differentiation in the direction of a given vector field reduces to the differentiation with respect to the scalar parameter t which now enters via the corresponding transformation and its gradient.It is shown in [33] that the shape derivative for a nonlinear PDE-constrained shape optimisation problem can be computed as the derivative of the Lagrangian with respect to the perturbation parameter.We will illustrate this systematic procedure for a number of different applications and utilise symbolic differentiation provided by the finite element software package NGSolve [30] to obtain the 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 value problems.The aim of this paper is to extend NGSolve by the possibility of semi-automatic and fully 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 and the cost function, • a semi-automated setting which offers a highly customizable user interface, but requires mathematical 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 will also lay the foundation for the following sections.In Section 3 we present a first unconstrained shape optimisation problem and show how to solve it in NGSolve.For this purpose we show how to compute the first and second order shape derivative in a semi-automated way.Section 4 extends the preceding section by incorporating a PDE constraint.The strategy is illustrated by means of a simple Poisson equation.We also show how to treat the computation of shape derivatives when the PDE is defined on surfaces.While the semi-automated shape differentiation presented in Sections 3 and 4 requires mathematical background knowledge, in Section 5 we show how the shape derivatives can be computed in a fully automated fashion.In the last section of the paper we verify the computed formulas by a Taylor test, discuss optimisation algorithms and present several numerical optimisation examples including nonlinear elasticity, Maxwell's equations and Helmholtz's equation.

A brief introduction to NGSolve
In this section, we give a brief overview over the main concepts of the finite element software NGSolve [30].We first describe the main principles for numerically solving boundary value problems in NGSolve before focusing on its built-in automatic differentiation capabilities.In the subsequent sections of this paper, these ingredients will be combined to implement the shape derivative of unconstrained and PDE-constrained shape optimisation problems in an automated way.

Solving PDEs with finite elements in NGSolve
In this section, we illustrate the syntax of NGSolve using the python programming language for the Poisson equation with homogeneous Dirichlet conditions as a model problem.We refer the reader to the online documentation https://ngsolve.org/docu/latest/for a more detailed description of the many features.
The weak form of the model problem on a domain We consider a ball of radius 1 2 in two space dimensions centered at the point (0.5, 0.5) , i.e.Ω = B((0.5,0.5) , 0.5), and the right hand side is defined by 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.The first line imports all modules from the package NGSolve.The second line includes the SplineGeometry function which enables us to define a mesh via a geometric description, in our case a circle centered at (0.5, 0.5) of radius 0.5.Finally the mesh is generated in line 7 and in line 8 we specify that we want to use a curved finite element mesh for a more accurate approximation of the geometry.
Next in line 9 we define an H 1 conforming finite element space of polynomial degree 3 and include Dirichlet boundary conditions on the boundary of the domain ∂ Ω (referenced by the string "circle" that we assigned in line 5).On this space we define a trial function u in line 11 and a test function w in line 12.These are purely symbolic objects which are used to define boundary value problems in weak form.For a more compact presentation later on, we define a coefficient function X which combines the three spatial components: Now, the left and right hand side of problem (2.1) can be conveniently defined as a bilinear or linear form, respectively, on the finite elements space fes by the following lines.
We assemble the system matrix coming from the bilinear form a and the load vector coming from L and solve the corresponding system of linear equations.Here, gfu is defined as a GridFunction over the finite element space fes.The Dirichlet conditions are incorporated into the direct solution of the linear system and the numerical solution is drawn in the graphical user interface.The numerical solution is depicted in Figure 1.

Automatic Differentiation in 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 in direction v := dir, that is, When calling the Diff command for expr, the expression tree of expr is gone through node by node, and for each node the corresponding differentiation rules such as product rule or chain rule are applied.When a node represents the variable with respect to which the differentiation is carried out, it is replaced by the direction dir of differentiation.
Figure 2 shows the differentiation of the expression expr= 2x*x+3y with respect to x into the direction given by v: The output of print(expr) will give the following output: Diff u wrt x: ConstantCF, val = 0 Diff w wrt x: ConstantCF, val = 0 Diff gf wrt x: ConstantCF, val = 0

Semi-automatic shape differentiation without constraints
We will illustrate the steps to be taken in order to obtain the shape derivative of a shape function in a semi-automatic way for a simple shape optimisation problem.For Ω ⊂ R d bounded and open and a continuously differentiable function f ∈ C 1 (R d ), we consider the shape differentiation of the shape function Clearly the minimiser of over all measurable sets in R d is given by We also refer to [25] for the computations of first and second order variations of type (3.1)where Ω is a submanifold of R d .

First order shape derivative
Given a vector field V ∈ C 0,1 (R d ) d , we define the transformation Definition 3.1.The first order shape derivative of a shape function at

Shape differentiation of unconstrained volume integrals
Using the transformation y = T t (x) and the notation F t := ∂ T t = I + t∂ V for the Jacobian of the transformation T t , we get for as in (3.1), Now let us explain how to compute the shape derivative of .Denoting the chain rule gives (formally) Using that , we get for the shape derivative This is the form we use for defining the first order shape derivative in NGSolve.Note that a Lipschitz vector field is differentiable almost everywhere and hence ∂ V (x) is defined almost everywhere and bounded.Given the function f Here, we introduce the symbol F and assign to it the value of the identity matrix in line 42.This allows 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 a vector-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 spatial variables in the direction given by V. Note that X is the coefficient function we introduced in line 13.In line 49, we deal with the differentiation with respect to F.
Therefore, we obtain for the first order shape derivative the well-known formula Finally if Ω is smooth enough (for instance C 1 ), then the shape derivative is given by where n denotes the outward pointing normal along ∂ Ω.

Shape differentiation of unconstrained boundary integrals
For Ω and f as in the previous section we consider

Second order shape derivatives
For second order shape derivatives, we consider perturbations of the form for s, t ≥ 0 and define Ω s,t := T s,t (Ω).
Definition 3.3.The second order shape derivative of a shape function at Remark 3.4.We remark that if is smooth enough, the second order derivative as defined in (3.8) is symmetric by definition: (3.9) We stress that this derivative is not the same as the shape derivative obtained by repeated shape differentiation, that is, it does not coincide with (see, e.g., [6, Chap.9, Sec.6]) which is in general asymmetric.However, in NGSolve we compute directly the second derivative as defined in (3.8).However, this derivative is only symmetric if ∂ V W = 0 since In NGSolve, when repeating the shape differentiation procedure introduced in Section 3.1, we compute directly the second order shape derivative as defined in (3.8).Here, we exploit that GridFunctions 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.12) Formula (3.12) is used for the automatic derivation of the second order shape derivative in NGSolve.
Remark 3.5.We remark that the formula (3.13) can be evaluated explicitly and reads Formula (3.13) can be implemented in NGSolve as follows: 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.

Semi-automatic shape differentiation with PDE constraints
In this section, we describe the automatic derivation of the shape derivative for the following type of equality constrained shape optimisation problems: where e : × X → X * with e(Ω, •) : X (Ω) → X (Ω) * represents an abstract PDE constraint with X = ∪ Ω∈ X (Ω) being the union of Banach spaces X (Ω) and a set of admissible shapes.For any given Ω ∈ we assume the PDE constraint (4.2) to admit a unique solution which we denote by u Ω .Moreover, let (Ω) := J(Ω, u Ω ) denote the reduced cost functional.One way to eliminate the constraint (4.2) is to introduce a Lagrangian 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: where Here the transformation Φ t depends on the differential operator involved.For instance Now the shape differentiability of (4.1)-(4.2) is reduced to proving that (see [34]) 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 accomplished by different methods, see e.g.[34] for an overview.However, we remark that (4.5) holds for a large class of nonlinear PDE constrained shape optimisation; see [33].
The rest of this section is organised as follows: We introduce a model problem, which is the 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 model problem can be obtained in NGSolve in Sections 4.2 and 4.3.Finally, we also briefly discuss the extension to partial differentiation equations on surfaces.

PDE-constrained model problem
We will illustrate the derivation of the first and second order shape derivative for the minimisation of a tracking-type cost functional subject to Poisson's equation on the unknown domain ⊂ (R d ) a set of admissible shapes.We consider the problem min The Lagrangian is given by Given an admissible shape Ω, a vector field the perturbed domain.Therefore the parametrised Lagrangian is given by Changing variables yields Recall that, for a given Ω ∈ , u Ω denotes the corresponding unique solution to (4.6b) and (Ω) the reduced cost functional, (Ω) := J(Ω, u Ω ).Let u t ∈ H 1 0 (Ω) the solution of the perturbed state equation brought back to the original domain Ω, that is, Note that, for u t defined by (4.10), it holds . It can easily be shown that (4.5) holds and thus the shape derivative in the direction of a vector field where p ∈ H 1 0 (Ω) denotes the adjoint state and is defined as the unique solution p ∈ H 1 0 (Ω) to or explicitly

First order shape derivative
By the discussion above, the first order shape derivative is given by ∂ t G(0, u, p) with G defined in (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 unconstrained problem We can compute explicitly Now we are in a position to compute the first order shape derivative for the PDE-constrained shape optimisation problem (4.6) in NGSolve.After solving the state equation as shown in Section 2.1, the adjoint equation can be solved as follows.We can now define the Lagrangian (4.9) such that the shape derivative can be obtained by the same procedure as in the unconstrained setting.Note that lines 82-83 coincide with lines 48-49.

Second order shape derivative
Let us introduce the notation where T s,t (x) = x + sV (x) + tW (x) and F s,t := ∂ T s,t .We observe that for s, t ≥ 0. In case t = 0 we write u s := u s,t | t=0 and p s := p s,t | t=0 and similarly for t = s = 0 we write u := u s,t | s=t=0 and p := p s,t | s=t=0 .Therefore, consecutive differentiation of (4.17) first with respect to t at zero and then with respect to s at zero yields where ∂ s u 0 ∈ H 1 0 (Ω) solves the material derivative equation or equivalently Similarly the function ∂ s p 0 ∈ H 1 0 (Ω) solves the material derivative equation . Formally (4.20) and (4.22) can be written as an operator equation 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 denote by X2.We define trial and test functions as well as grid functions representing the deformation vector fields V and W , which we initialise with some functions.
We define a 2×2 block bilinear form as well as a 2×1 block linear form which will represent the left and right hand sides of (4.23), respectively.The operator equation in (4.23) can be conveniently defined by differentiating the Lagrangian with respect to the corresponding variables.
#l i n e 2 shapeGradLag2 += ( G_pde .D i f f ( gfp , p T e s t ) ) .D i f f ( F , grad ( gfV ) ) shapeGradLag2 += ( G_pde .D i f f ( gfp , p T e s t ) ) .D i f f (X , gfV ) We can solve this combined system for ∂ s u 0 and ∂ s p 0 and access and visualise the two components in the following way: 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 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 remains to evaluate the term (4.19).We define the three terms of (4.19) as bilinear forms, assemble them and perform vector-matrix-vector multiplications:

PDEs on surfaces
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 3  and denote by n the normal field along M .Let u d ∈ H 1 (R d ) be given and define where u ∈ H 1 (M ) solves the surface equation where ∇ M ϕ denotes the tangential gradient of ϕ; see [6, p.493, Def.5.1].We assume that the function f ∈ H 1 (R 3 ) is given.The Lagrangian is given by As in the previous section we fix an admissible shape M and let M t := (Id + t V )(M ) be a small perturbation of M by means of a vector field V ∈ C 1 (R d ) d for t > 0 small.The parametrised Lagrangian is given by Define the density ω(F t ) := det(F t )|F − t n|.Changing variables and using where , u, p) we obtain in analogy to the domain case We can compute explicitly where ∂ M V denotes the tangential Jacobian of V and div M (V ) := ∂ M V : I the tangential divergence, 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 order derivatives here.We first define the geometry of the unit sphere, create a surface mesh and define a finite element space on the surface mesh: #s e t up and s o l v e s t a t e e q u a t i o n 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: The shape derivative is obtained as in the case of PDEs posed on volumes by the evaluation of (4.30):

Fully automated shape differentiation
In the previous sections we used the automatic differentiation capabilities of NGSolve to alleviate the shape differentiation procedure.However, so far we still had to include some knowledge about the problems at hand.So far, it was necessary to define the objective function or Lagrangian G in the correct way, accounting for the correct transformation rules between perturbed and unperturbed domain.In this section, we will show that also this step can be automated since all necessary information is already included in the functional setting.The fully automated shape differentiation is incorporated by the command DiffShape(...).
In particular, in the fully automated setting it is enough to set up the cost function or Lagrangian for the unperturbed setting.For a shape function of the type (3.1) we can define the shape derivative of the cost function in the following way: Note that there is no term of the form Det(F) showing up in line 186.Here, the transformation of the domain is taken care of automatically.It can be checked that this really gives the same result as dJOmega_f defined in lines 48-49.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 rules used.This can be cured by enforcing an integration rule of higher order for G_f, i.e. by replacing the 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 very similar.Here the idea exploited in the implementation of the command DiffShape(...) is to just differentiate the general expression (4.4) with respect to the parameter t.The transformations Φ t appearing in (4.4), which depend on the functional setting of the PDE, are identified automatically 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.
as well as their respective second order derivatives.Remark 5.1.We remark that the fully automated differentiation using DiffShape(...) should be seen to complement the semi-automated shape differentiation techniques introduced in Sections 3 and 4 rather than to replace them.Using the semi-automated differentiation, the user has the possibility to, on the one hand, keep control over the involved terms, and on the other hand also to adjust the shape differentiation to their custom problems which may be non-standard.As an example where the semi-automated differentiation may be beneficial compared to the fully automated differentiation we mention the case of time-dependent PDE constraints considered in a space-time setting when a shape deformation is only desired in the spatial coordinates.Of course, when one is interested in the shape derivative for a more standard problem, the fully automated way appears to be more convenient and less error prone.

Code verification
We verify the expressions that we obtained in a semi-automatical or fully automatical way for the first and second order shape derivatives by looking at the Taylor expansions of the perturbed shape functionals.We illustrate our findings in two examples in R 2 .On the one hand, we consider a shape function as introduced in (3.1) with an additional boundary integral as in (3.6), henceforth denoted by 1 ; on the other hand, we consider the PDE-constrained shape optimisation problem defined by (4.6), the reduced form of which will be denoted by 2 (Ω).More precisely, we consider where u Ω solves (4.6b).(6.2) In the case of 1 , we used the function f (x 1 , x 2 ) = 0.5 and for 2 , we used For the test of the first order shape derivatives D i (Ω)(V ) we choose a fixed shape Ω and a vector field V ∈ C 0,1 (R 2 ) 2 and observe the quantity for t 0. Likewise, for the second order shape derivative, we consider the remainder as t 0. By the definition of first and second order shape derivatives, it must hold that This behavior can be observed in Figure 3(a) for 1 and in Figure 3(b) for 2 , where we used ) in both cases.The experiments for shape function 1 was conducted on a mesh consisting of 13662 vertices, 26946 elements and with polynomial order 2 (resulting in 54269 degrees of freedom), and the experiment for 2 with 95556 vertices and 190062 elements and polynomial degree 1 (95556 degrees of freedom).We conducted these experiments for a number of different problems with different 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.

Optimisation algorithms
In this section we discuss how to use optimisation algorithms in conjunction with the automated shape differentiation explained in the previous sections.The starting point of our discussion is a fixed initial shape Ω.Then we consider the mapping defined on a suitable space of vector fields Θ ⊂ C 0,1 (D) d .Since the mapping g is defined on an open subset Θ of the Banach space C 0,1 (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 be difficult in practice.We observe that for V, W ∈ Θ we have ). (6.7)

Gradient computation
The gradient of ∂ g(V ) in a Hilbert space H ⊂ C 0,1 (D) d is defined by Typical choices for H are ) ) 12) The last choice corresponds to a penalised Cauchy-Riemann gradient and results in a gradient which is approximately conformal and hence preserves good mesh quality.We refer to [17] for a detailed description.

Basic algorithm
Let Ω be an initial shape and let H ⊂ C 0,1 (D) d be a Hilbert space.Then a basic shape optimisation algorithm reads as follows.

end if 11: end while
We present and explain the numerical realisation of Algorithm 1 in NGSolve for the case of a PDE-constrained shape optimisation problem in two space dimensions.The simpler case of an unconstrained shape optimisation problem or the case of three space dimensions can be realised by small modifications of the presented code.
First of all, we mention that we realise shape modifications in NGSolve by means of deformation vector fields without actually modifying the coordinates of the underlying finite element grid.Recall the vector-valued finite element space VEC over a given mesh as introduced in code line 44.We define a vector-valued GridFunction with the name gfset which will represent the current shape.We initialise it with some vector-valued coefficient function ) 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 ( d e f o r m a t i o n=True ) g f s e t .S e t ( ( X

S e t D e f o r m a t i o n ( g f s e t ) Redraw ( )
Any operation involving the mesh such as integration or assembling of matrices is now carried out for the deformed configuration.The deformation can be unset by the command mesh.UnsetDeformation().Integrating the constant function over the mesh in the perturbed and 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 1.7924529046862627 0.7854072970684544 respectively.
In the course of the optimisation algorithm the state equation as well as the adjoint equation have to be solved for every new shape.We define the following function, which computes the state and adjoint state for a linear PDE constraint: The shape derivative dJOmega for some problem at hand can be defined as illustrated in Sections 4.1 and Section 5. Finally, we need to define the shape gradient, which is the solution to a boundary value problem of the form (6.8).We choose the bilinear form defined in (6.11) with γ CR = 10:  Mesh movement and mesh optimisation As an alternative to realizing the deformations via mesh.SetDeformation(...), where the underlying mesh is not modified, one could also just move every mesh node in the direction of the given descent vector field by changing its coordinates.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 ( ) : mesh .ngmesh .Update ( ) Here, the displacement vector field displ, which is of type GridFunction, is evaluated for each mesh node and, subsequently, the mesh nodes are updated.At the end of the procedure, the mesh structure needs to be updated, see line 291.Note that the evaluation of GridFunctions requires 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 the method mesh.ngmesh.OptimizeMesh2d() followed by mesh.ngmesh.Update().Figure 4 shows a distorted mesh and the result of a call of mesh.ngmesh.OptimizeMesh2d().

Newton's method for unconstrained problems
The particular choice H = H 1 0 (D) d and leads to Newton's method.We refer to [1,8,22,24] where shape Newton methods were used previously and to [14, Chapter 2] and [18, Chapter 5] for Newton's method in an optimal control setting.This bilinear form is only positive semi-definite on Moreover, from the structure theorem for second shape derivatives proved in [21] we know that at a stationary point Ω, that is, Figure 4: Before and after mesh optimisation by mesh.ngmesh.OptimizeMesh2d(). where As a result the gradient according to (6.13) is not uniquely determined.To get around this difficulty, the shape Hessian is often regularised by an H 1 term, i.e. (6.13) is replaced by see, e.g.[27], which, however, impairs the convergence speed of Newton's method.

Alternative regularisation strategy.
Here, we propose the following strategy: We regularise the shape Hessian only on the boundary ∂ Ω and only in tangential direction, i.e., we choose with a regularisation parameter δ.To exclude the part of the kernel corresponding to interior deformations, 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 those on the boundary.As a result, we get a shape gradient ∇ (Ω) which is nonzero only on the boundary.We extend this vector field to the interior by solving an additional boundary value problem (of linearised elasticity type), where we use the deformation given by ∇ (Ω) as Dirichlet boundary conditions.The Newton algorithm reads as follows.

Newton's method for PDE-constrained problems
We consider the PDE-constrained model problem of Section 4.1 which is subject to the Poisson equation.The unregularised Newton system reads In Subsection 4.3 we discussed how the second order derivative can be evaluated.Recalling that d (Ω)(V ) = ∂ s G V,0 (0, 0, u, p) we see that (4.23) and (4.19) lead to 19) The component Ṽ then represents the direction which we use for the shape Newton optimisation step.The matrix in (6.19) can be realised in NGSolve by using a combined finite element space X3 consisting of three components as follows: The right hand side of (6.19) can be defined as follows: Recall that the system (6.15) has a nontrivial kernel as discussed in Section 6.2.3.This problem can be circumvented by proceeding like in the unconstrained case.We add a regularisation only on the boundary, and solving the regularised system using these free dofs: The newton direction is then given as the first of the three components of the obtained solution.g f s e t .vec .data = g f s e t .vec + 1 * V t i l d e .vec

Numerical shape optimisation of model problems
In this section we will apply the automated shape differentiation and all numerical algorithms introduced in the preceding sections in numerical examples.

A first shape optimisation problem
In this section, we revisit problem (3.1) introduced in Section 3, i.e. the problem of finding a shape Ω such that the cost function (Ω) = ´Ω f (x) dx is minimised.

First order methods
We illustrate our first order methods in a problem which was also considered in [17] and reproduce the results obtained there.We choose the function with a = 4 5 , b = 2 and ε = 0.001.Recall that the optimal shape is given by {(x 1 , x 2 ) ∈ R 2 : f (x 1 , x 2 ) < 0} which is depicted in Figure 5 (right).We start our optimisation algorithm with the unit disk, Ω 0 = B 1 (0) as an initial design.Note that the optimal design cannot be reached by means of shape optimisation using boundary perturbations.However, we expect the outer curve of the optimal shape to be reached very closely.
We apply Algorithm 1 with the shape gradient ∇ associated to the H 1 inner product (6.9), to the bilinear form of linearised elasticity (6.10) and including the additional Cauchy-Riemann term (6.11).We chose the algorithmic parameters γ = 1e − 4, ε = 1e − 7, a mesh consisting of 2522 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 neighborhood of the optimal solution, we choose a simpler optimal design here.We choose which yields an ellipse with the lengths of the two semi-axes a and b.We choose a = 1.3 and b = 1/a and again start the optimisation with the unit disk as initial shape.Figure 9 shows the initial and optimised design after only six iterations of Algorithm 2 with (•, •) H chosen as in (6.17) with δ = 100.A comparison of the convergence histories between the choice (6.17) with δ = 100 and (6.16) with δ = 0.5 is shown in the right picture of Figure 9.In both cases, the parameter δ was chosen empirically to get convergence as fast as possible.The experiments were conducted on a finite element mesh consisting of 2522 nodes and 4886 triangular elements with a finite element space VEC of order 3, with the algorithmic parameter ε = 10 −7 .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: Objective value and norm of shape gradient ∇ (Ω) in the course of second order optimisation using (6.16) with δ = 0.5 and (6.17) with δ = 100.

Shape optimisation subject to the Poisson equation
In this section, we revisit the model problem introduced in Section 4.
Note that the data is chosen in such a way that, for Ω * = (0, 1) 2 it holds (Ω * ) = 0 and thus Ω * is a global minimiser of .We show results obtained by first and second order shape optimisation methods exploiting automated differentiation.
We ran the optimisation algorithm in three versions.On the one hand, we applied a first order method with constant step size α = 1.On the other hand, we applied two second order methods with 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 method performs as well as possible.In the case of (6.16) we chose δ = 0.001 and in the case of (6.17) δ = 1.The experiments were conducted on a finite element mesh consisting of 4886 elements with 2522 vertices and polynomial degree 1.In Figure 10, we can observe the decrease of the objective function as well as of the norm of the shape gradient over 200 iterations for these three algorithmic settings.
Figure 10 shows the initial design as well as the design after 200 iterations of the second order method with regularisation strategy (6.17).Note that the improved design is very close to the global solution Ω * = (0, 1) 2 .The initial design was chosen as the disk of radius 1  2 centered at the point 1  2 , 1

Nonlinear elasticity
Here, we illustrate the applicability of the automated shape differentiation and optimisation in the more realistic and more complicated setting of nonlinear elasticity in two space dimensions using a Saint Venant-Kirchhoff material with Young's modulus E = 1000 and Poisson ratio ν = 0.3.We consider a two-dimensional cantilever which is clamped on the upper and lower left parts of the boundary, Γ 1 l = {0} × (0.88, 1) and Γ 2 l = {0} × (0, 0.12), respectively, and is subject to a surface force g N = (0, −100) on Γ r = {1} × (0.45, 0.55).The initial geometry with 3 holes is depicted in Figure 12  for all v ∈ H 1 Γ l (Ω) 2 .Here, S(u) denotes the Saint Venant-Kirchhoff stress tensor where C(u) = (I 2 + ∇u) (I 2 + ∇u) and I 2 is the identity matrix, see also [2, Sec.8], and λ and µ denote the Lamé constants,

.24)
We minimise the functional with α = 2.5 subject to (6.22) which amounts to maximising the structure's stiffness while bounding the allowed amount of material used.We remark that the well-posedness of (6.22) is not clear, see also the discussion in [2, Sec.8].Nevertheless, application of the automated shape differentiation and optimisation yields a significant improvement of the initial design.The highly nonlinear PDE constraint (6.22) is solved by Newton's method.In order to have good starting values, a load stepping strategy is employed, i.e., the load on the right hand side is gradually increased, the PDE is solved and the solution is used as an initial guess for the next load step.This is repeated until the full load is applied.With these ingredients at hand, Algorithm 1 (i.e.code lines 244-284) can be run.We chose the algorithmic parameters 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 an additional Cauchy-Riemann term as in (6.11) with weight γ CR = 10.The objective value was reduced from 3.125 to 2.635 (volume term from 1.290 to 1.096) in 15 iterations of Algorithm 1.The results were obtained on a mesh consisting of 10614 elements and 5540 vertices using piecewise linear, globally continuous finite elements.

Helmholtz equation
In this section, we consider the problem of finding the optimal shape of a scattering object.More precisely, we consider the minimisation of the functional ˆΓr uu ds (6.26) subject to the Helmholtz equation with impedance boundary conditions on the outer boundary, for all w ∈ H 1 (Ω, 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 hand side is chosen as f (x 1 , x 2 ) = 10 3 • e −9((x 1 −0.2) 2 +(x 2 −0.5) 2 ) , (6.28) see Figure 13(a).Furthermore Ω = B((0.5,0.5) , 1)\B((0.75, 0.5) , 0.15) denotes the domain of the right half of the outer boundary.Here, only the inner boundary ∂ Ω\Γ is subject to the shape optimisation.Thus, the aim of this model problem is to find a shape of the scattering object such that 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 scattering object) and for the optimised configuration.The objective value was reduced from 3.44 • 10 −3 to 3.31 • 10 −3 .The forward simulations were performed using piecewise linear finite elements on a triangular grid consisting of with 34803 degrees of freedom.The optimisation stopped after 12 iterations.

Application to Electrical Machine
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 3 denote the computational domain, which consists of ferromagnetic material, air regions and permanent magnets.denotes the magnetic reluctivity, which is a nonlinear function ν inside the ferromagnetic regions and equal to a constant ν 0 elsewhere.Further, δ > 0 is a small regularisation parameter and M : D → R 3 denotes the magnetisation in the permanent magnets.The nonlinear function ν satisfies a Lipschitz condition and a strong monotonicity condition such that problem (6.30) is well-posed.We refer the reader to [11,Sec. 6] for a more detailed description of the problem and [10] for a 2D version of the same problem.
As mentioned in Section 4, the transformation Φ t used in (4.4) depends on the differential operator.For the curl-operator, we have see e.g.[20, Section 3.9].Thus, the variational equation (6.30) can be defined as follows.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 part of 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 the machine, is given by "air|airgap".Moreover, nuIron denotes the nonlinear reluctivity function ν and magn contains the magnetization direction of the permanent magnets.Likewise, the objective function can be defined as follows, where n2D and Bd represent the extension of the normal vector to the interior of the air gap and the desired curve, respectively.For the definition of all quantities, we refer the reader to the code which was submitted together with this manuscript.The shape differentiation as well as the optimisation loop now works in the same way as in the previous examples.Figure 15 shows the initial design of the motor as well as the optimised design obtained after 11 iterations of Algorithm 1 with γ = 0.The experiment was conducted using a tetrahedral finite element mesh consisting of 13440 vertices, 57806 elements and Nédélec elements of order 2 (resulting in a total of 323808 degrees of freedom).The objective value was reduced from 2.5944 • 10 −8 to 4.565 • 10 −10 in the course of the first order optimisation algorithm after 11 iterations.It can be seen from Figure 16 that the difference between the quantity curl(u) • n and the desired curve B n d inside the air gap decreases significantly.(b) History of objective value and norm of shape gradient using a first order algorithm with line search.

Surface PDEs
Finally, we also show the application of a shape optimisation algorithm to a problem constrained by a surface PDE.We solve problem (4.24)-(4.25)with u d = 0, f (x 1 , x 2 , x 3 ) = x 1 x 2 x 3 and initial shape M = S 2 the unit sphere in R 3 .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 of the norm of the shape gradient.The objective value was reduced from 7.08 • 10 −4 to 9.88 • 10 −9 .
Figure 18 shows the final design which was obtained after 575 iterations from two different perspectives.The experiment was conducted using a surface mesh with 332 vertices and 660 faces and polynomial degree 3 (resulting in 2972 degrees of freedom).

Conclusion
We showed how to obtain first and second order shape derivatives for unconstrained and PDEconstrained shape optimization problems in a semi-automatic and fully automatic way in the finite element software package NGSolve.We verified the proposed method numerically by Taylor 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 shape optimisation to further improve numerical methods on the one hand, and product engineers working with NGSolve to design devices in an optimal fashion on the other hand.
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=0.2) ) mesh .Curve( 3 ) 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 ( )
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 " )

Figure 2 :
Figure 2: Illustration of Diff command for example expr= 2x*x+3y.(a) Expression tree for expr.(b) Expression tree for expression obtained by call of expr.Diff(x, dir).

2 F
3, b = 1/a and R = 0.5, we implement the transformed cost function (3.3) as follows: f = ( ( X[0] −0.5) / 1 .3 ) * * 2+(1.3 * (X[1] −0.5) ) * * 2 − 0.5 * * = I d ( 2 ) # s y m b o l i c i d e n t i t y m a t r i x G_f = f * Det ( F ) * dx # F o n l y a c t s as a dummy v a r i a b l e VEC = VectorH1 ( mesh , o r d e r =1, d i r i c h l e t=" " ) #v e c t o r i a l FE s p a c e o f o r d e r 1 V = VEC .T e s t F u n c t i o n ( ) Cost ( u ) : r e t u r n ( u−ud ) * * 2 * Det ( F ) * dx #s o l v e a d j o i n t e q u a t i o n g f p = 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 ( ) g f p .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 " ) d e f Equation ( u ,w) : r e t u r n ( ( I n v ( F ) .t r a n s * grad ( u ) ) * ( I n v ( F ) .t r a n s * grad (w) ) − f * w) * Det ( F ) * dx G_pde = Cost ( gfu ) + E qu at io n ( gfu , g f p ) 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) ) 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 * g f d s u .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 ) 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 ( ) s p h e r e = Sphere ( Pnt ( 0 , 0 , 0 ) , 1 ) .bc ( " o u t e r " ) g e o _ s u r f .Add( s p h e r e ) 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 =3,maxh=0.2) ) mesh_surf .Curve ( 3 ) f e s _ s u r f = H1( mesh_surf , o r d e r = 3) Next we define the transformed cost function and partial differential equation needed for setting up the Lagrangian (4.29).Here, we again make use of a symbolic object F to which we assign the 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 [ 0 ] * X [ 1 ] * X [ 2 ] ) F = I d ( 3 ) tangDet = Det ( F ) * Norm( I n v ( F ) .t r a n s * s p e c i a l c f .normal ( 3 ) ) Bmat = ( I d ( 3 ) − 1/Norm( I n v ( F ) .t r a n s * s p e c i a l c f .normal ( 3 ) ) * * 2 * OuterProduct ( I n v ( F ) .t r a n s * s p e c i a l c f .normal ( 3 ) , I n v ( F ) .t r a n s * s p e c i a l c f .normal ( 3 ) ) ) * I n v ( F ) .t r a n s 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 * * 2 * tangDet * ds Now we can define the bilinear form and solve the state equation.Here, the right hand side of the 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.
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 " ) 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 lines 58-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) d e f Cost_0 ( u ) : r e t u r n ( u−ud ) * * 2 * 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 , g f p ) dJOmega_pde_0 = LinearForm (VEC) dJOmega_pde_0 += G_pde_0 .D i f f S h a p e (V) Here, gfu and gfp represent the solutions to the state and adjoint equation, respectively, and must have been computed previously.The bilinear form shapeHess11 used in Section 4.3 (see lines 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 * * 2 * 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 ( )
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 g f p .vec .data = −i n v a .T * dCostdu .vec d e f eps ( u ) : r e t u r n 1/2 * ( 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 ( ) # d e f i n e t r i a l f u n c t i o n W and t e s t f u n c t i o n V 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 += 10 * ( grad Now we can run Algorithm 1 for problem (4.6): alpha = 1 a l p h a _ i n c r _ f a c t o r = 1 .2 gamma = 1e−4 Nmax = 100 e p s i l o n = 1e−7 i s C o n v e r g e d = F a l s e g f s e t .S e 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 .S e t D e f o r m a t i o n ( 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 ) w h i l e True : i f currentNormGFX < e p s i l o n : i s C o n v e r g e d = True break gfsetTemp .vec .data = g f s e t .vec − alpha * gfX .vec mesh .S e t D e f o r m a t i o n ( gfsetTemp ) solvePDE ( ) Jnew = I n t e g r a t e ( Co st ( gfu ) , mesh ) mesh .UnsetDeformation ( ) i f Jnew < J o l d − gamma * alpha * currentNormGFX * * 2 : 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 / 2 Redraw ( b l o c k i n g=True )

VEC2=
VectorH1 ( mesh , o r d e r =1, d i r i c h l e t = " c i r c l e " ) #a u x i l i a r y s p a c e f o r boundary c o n d i t i o n s 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 += 100 * 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

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 e t ( gfX_bnd ) aX_ext .Assemble ( ) r = gfX_bnd .vec .C r e a t e V e c t o r ( ) r .data = (−1) * 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 .S e t ( ( 0 , 0 ) ) g f s e t .vec .data = g f s e t .vec − 1 * gfX .vec

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 ) #b l o c k ( 1 , 1 )

d e l t a = 1 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×3 block system.This can be realised by setting Dirichlet boundary conditions for the interior degrees of freedom, i.e. by dealing with the free degrees of freedom, # copy o f VEC with D i r i c h l e t boundary c o n d i t i o n s on whole boundary : VEC2 = VectorH1 ( mesh , o r d e r = 1 , d i r i c h l e t = " .* " ) freeDofsCombined = B i t A r r a y (VEC2 .ndof + 2 * 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 ] 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 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 ld e _ b n d .vec .data = gfCombined3 .components [ 0 ] .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 )

Figure 6 :
Figure 6: Results of problem (3.1) with f as in(6.20)and the shape gradient associated to the H 1 inner product (6.9).

Figure 7 :
Figure 7: Results of problem (3.1) with f as in (6.20) and the shape gradient associated to the elasticity bilinear form (6.10).

Figure 8 :
Figure 8: Results of problem (3.1) with f as in (6.20) and the shape gradient associated to the elasticity bilinear form with Cauchy-Riemann term (6.11).

Figure 10 :
Figure 10: Convergence behaviour for shape optimisation problem (4.6) with proposed regularisation strategies (6.17) and (6.16) as well as first order method with constant step size α = 1.(a) Behaviour of objective function .(b) Behaviour of norm of shape gradient ∇ (Ω) .

Figure 12 :
Figure 12: Initial and optimised geometry of cantilever under vertical force on right hand side using St. Venant-Kirchhoff model in nonlinear elasticity.(a) Initial geometry.(b) Optimised geometry (reference configuration).(c) Optimised geometry (deformed configuration).
Our aim is to minimise the functionalˆΩg | curl u • n − B n d | 2 dx,(6.29)whereΩ g denotes the air gap region of the machine, n denotes an extension of the normal vector to the interior ofΩ g , B n d : Ω g → R 3 is a given smooth function and u ∈ H 0 (curl, D) is the solution to the boundary value problem ˆD ν Ω (| curl u|) curl u • curl w + δu • w dx = ˆΩm M • curl w dx (6.30)for all w ∈ H 0 (curl, D).Here, Ω ⊂ D denotes the union of the ferromagnetic parts of the electrical machine, Ω m denotes the permanent magnets subdomain and ν Ω = χ Ω (x)ν(| curl u|) + χ D\Ω (x)ν 0(6.31)

Figure 17 :
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 line search.
[32,e e.g.[32, Prop.2.47], with the outer unit normal vector n and | • | denoting the Euclidean norm.Again, the shape derivative can be computed as the total derivative of this expression with respect to the parameter t.In NGSolve, the only difference lies in the necessity to use the trace of the gradient of a test vector field V.