A Unified Approach to Shape and Topological Sensitivity Analysis of Discretized Optimal Design Problems

We introduce a unified sensitivity concept for shape and topological perturbations and perform the sensitivity analysis for a discretized PDE-constrained design optimization problem in two space dimensions. We assume that the design is represented by a piecewise linear and globally continuous level set function on a fixed finite element mesh and relate perturbations of the level set function to perturbations of the shape or topology of the corresponding design. We illustrate the sensitivity analysis for a problem that is constrained by a reaction–diffusion equation and draw connections between our discrete sensitivities and the well-established continuous concepts of shape and topological derivatives. Finally, we verify our sensitivities and illustrate their application in a level-set-based design optimization algorithm where no distinction between shape and topological updates has to be made.


Introduction
Numerical methods for the design optimization of technical systems are of great interest in science and engineering.Applications include the optimization of mechanical structures [25,2], electromagnetic devices [16,6], fluid flow [18], heat dissipation [19] and many more.There exist several different approaches to computational design optimization.On the one hand, shape optimization techniques based on the mathematical concept of shape derivatives [13] can modify boundaries and material interfaces in a smooth way, but typically cannot alter the topology of a design.An exception being the level set method for shape optimization [2] where the design is represented by the zero level set of a design function whose evolution is guided by shape sensitivities via a transport equation.While this approach allows for merging of components, it lacks a nucleation mechanism and is often coupled with the topological derivative concept [26,24], see e.g.[9,1].In the class of density-based topology optimization methods [7], a design is represented by a density function ρ(x) that is allowed to attain any value in the interval [0, 1].Then, regions with ρ(x) = 0 and ρ(x) = 1 are interpreted as occupied by material 1 and 2, respectively, while intermediate density values 0 < ρ(x) < 1 are penalized in order to obtain designs that are almost "black-and-white".One advantage of density based methods is that the system response depends continuously on ρ and the standard notions of derivatives in vector spaces can be applied.Here, interfaces are typically not crisp and there is no measure of optimality with respect to shape variations at the interface.Finally we mention the level-set algorithm for topology optimization introduced in [4], where the design is guided solely by the topological derivative, which however is not defined on the material interfaces.As a consequence, the final designs cannot be shown to be optimal with respect to shape variations at the interface.This aspect has been thorougly analyzed in [5] where the authors draw a connection to density-based methods and, for two particular problem classes, propose an interpolation scheme which relates the derivative with respect to the density function, to topological and shape derivatives in the interior and on the interface, respectively.The goal of this paper is to unify the concepts of topological and shape perturbations and to treat design optimization problems by a unified sensitivity, called the topological-shape derivative.
In this way, we aim at combining topological sensitivity information (related to the topological derivative) in the interior of each subdomain and shape sensitivity information (related to the shape derivative) at the material interface.While the topological derivative is defined as the sensitivity of a design-dependent cost function with respect to the introduction of a small hole or inclusion of different material, the shape derivative is defined as the cost function's sensitivity with respect to a transformation of the domain.In order to unify these two concepts, we consider a domain description by means of a continuous level set function which attains positive values in one of two subdomains and negative values in the other.Then a perturbation of the level set function in the interior of a subdomain can be related to a topological perturbation, and a perturbation close to the material interface can be seen as a perturbation of the shape of the domain.We remark that this point of view is in alignment with the concept of dilations of points and curves as introduced in [11], see also [12].In principle, this idea was already followed in [8], however only for the case of shape optimization and not in combination with topology optimization.In [20] the author represents domains by level set functions and relates shape and topological derivatives of shape functionals to derivatives with respect to the level set function in a continuous setting without PDE constraints.In contrast to this, we consider PDE-constrained problems, but our analysis is performed on the discrete level, i.e. we follow the paradigm "discretize-then-optimize" for our sensitivity analysis with respect to a level set function.The rest of this paper is organized as follows: In Section 2, we introduce the model problem and the classical concepts of topological and shape derivative in the continuous setting.After presenting the discretized setting in Section 3, we proceed to compute the numerical topological-shape derivative of our discretized model problem in Section 4. In Section 5 we compare the computed sensitivities with the sensitivities obtained by discretizing the continuous formulas.Finally we verify our computed formulas and present optimization results in Section 6 before giving a conclusion in Section 7.

Model problem and continuous setting
where Here, A denotes a set of admissible subsets of D, and the data given.The weak formulation of the PDE constraint reads with admits a unique solution which we denote by u(Ω).We introduce the reduced cost function g(Ω) := g(Ω, u(Ω)).

Classical topological derivative
Let ω ⊂ R d with 0 ∈ ω.For a point z ∈ Ω ∪ (D \ Ω), let ω ε := z + εω denote a perturbation of the domain around z of (small enough) size ε and of shape ω.The continuous topological derivative of the shape function g = g(Ω) is defined by Note that the topological derivative is not defined for points z ∈ ∂Ω on the material interface.For problem (2) we obtain for z ∈ D \ Ω whereas for z ∈ Ω see, e.g.[3].

Classical shape derivative
We recall the definition of the classical shape derivative as well as its formula for our model problem (1)- (2).Given an admissible shape Ω ∈ A and a smooth vector field V ∈ C ∞ c (D) that is compactly supported in D, we define the perturbed domain for a small perturbation parameter t > 0 where id : R d → R d denotes the identity operator.The classical shape derivative of g at Ω with respect to V is then given by if this limit exists and the mapping V → dg(Ω)(V ) is linear and continuous.Under suitable assumptions it can be shown that this shape derivative admits the tensor representation for some tensors [21].Here, ∂V denotes the Jacobian of the vector field V .The structure theorem of Hadamard-Zolésio [13, pp. 480-481] states that under certain smoothness assumptions the shape derivative of a shape function g with respect to a vector field V can always be written as an integral over the boundary of a scalar function L times the normal component of V , i.e., where n denotes the unit normal vector pointing out of Ω.For problem (2) one obtains [21] S where I ∈ R d,d denotes the identity matrix, and Here, S Ω,in 1 and S Ω,out 1 denote the restrictions of the tensor S Ω 1 to Ω and D \ Ω, respectively.Furthermore, for two column vectors a, b ∈ R d , a ⊗ b = ab ∈ R d×d denotes their outer product, τ denotes the tangential vector and p ∈ H 1 0 (D) is the solution to the adjoint equation Moreover, motivated by the definition of the topological derivative (4) with the volume of the difference of the perturbed and unperturbed domains in the denominator, we introduce the alternative definition of a shape derivative with the symmetric difference of two sets A B := (A \ B) ∪ (B \ A).Note that the volume of the symmetric difference in (13) can be written as The proof is given in Appendix A.1.From Lemma 2.1, we immediately obtain the following relation between dg(Ω)(V ) and dg(Ω)(V ).
Corollary 2.2.Suppose that g is shape differentiable at Ω and that Ω and V are smooth and Proof.This follows immediately from the definition of dg(Ω)(V ) by Lemma 2.1 since

The continuous topological-shape derivative
Here and in the following, we assume that the domain Ω is described by a level-set function φ : For given φ, let Ω(φ) denote the unique domain defined by (18a)-(18c).In this section, in contrast to the setting in Section 2.2, we perturb Ω indirectly by perturbing φ such that φ ε = O ε φ for some operator O ε : C 0 (D) → C 0 (D) depending on ε ≥ 0 with the property Ω(O 0 φ) = Ω(φ).Later on, in the discrete setting, we will distinguish between two different types of perturbation operators O ε corresponding to shape or topological perturbations of Ω. Let, from now on, J (φ) := g(Ω(φ)) denote the reduced cost function as a function of the level set function φ.This way, a continuous topological-shape derivative can be defined as Note that this sensitivity depends on the choice of the perturbation operator O ε , which can represent either a shape perturbation or a topological perturbation.We will mostly be concerned with its discrete counterpart, which will be introduced in Section 4. Note that, in the case of shape perturbations, due to the scaling |Ω(φ ε ) Ω(φ)| instead of ε in the denominator the shape derivative is modified and does not correspond to (8) but rather to (13).
Relation to literature The sensitivity of shape functions with respect to perturbations of a level set function (representing a shape) was investigated in [20] for the case without PDE constraints.There, the author considers smooth level set functions and rigorously computes the Gâteaux (semi-)derivative in the direction of a smooth perturbation of the level set function, both for the case of shape and topological perturbations.In the case of shape perturbations, it is shown that the Gâteaux derivative coincides with the shape derivative (8) with respect to a suitably chosen vector field.On the other hand, a resemblance between the notions of Gâteaux derivative and topological derivative is shown, yet the Gâteaux derivative may vanish or not exist in cases where the topological derivative is finite.Evidently, this discrepancy results from the fact that the denominator in the definition of the Gâteaux derivative is always of order one whereas it is of the order of the space dimension in the topological derivative.
While the analysis for shape and topological perturbations is carried out separately in [20], a more unified approach is followed in [11,12].In these publications, the idea is to consider sensitivities with respect to domain perturbations that are obtained by the dilation of lower-dimensional objects.
Here, given a set E ⊂ R d of dimension k ≤ d, the dilated set of radius r is given by E r = {x ∈ R d : d E (x) ≤ r} where d E (x) denotes a distance of a point x to a set E. For instance, when E is chosen as a single point, the dilated set is just a ball of radius r around that point and performing a sensitivity analysis with respect to the volume of the dilated object leads to the topological derivative.On the other hand, when E is chosen as the boundary of a domain, E r can be defined using a signed distance function and corresponds to a uniform expansion of the domain.Then, a similar procedure leads to the shape derivative with respect to a uniform expansion in normal direction (i.e.V = n in ( 7)-( 8)).
In [11], a sensitivity analysis for various choices of E is carried out with respect to the volume of the perturbation.We note, however, that arbitrary shape perturbations are not covered and would require an extension of the theory.Comparing [20] and [11], we observe that in the former paper only smooth perturbations of a level set function are admissible whereas, in the latter approach, domain perturbations by dilations can be interpreted as perturbations of level set functions by a (non-smooth) distance function.Finally, we mention [8] where a domain is represented by a discretized level set function and a shape sensitivity analysis is carried out with respect to a perturbation of the level set values close to the boundary.This procedure can be interpreted as an application of the idea of dilation to discretized shape optimization problems.As the authors point out, this kind of shape sensitivity analysis is more natural compared to the standard approach based on domain transformations when employed in a level-set framework; an observation also made in [20,Sec. 3].The authors show numerical results for the shape optimization of an acoustic horn, but do not consider topological perturbations in this work.
As it can be seen from ( 19), our approach is related to the dilation concept since we also consider the sensitivity with respect to the volume of the domain perturbation Ω t Ω.In the following, we will investigate the topological-shape derivative in a discretized setting.Similarly to [8], we will consider shape sensitivity analysis with respect to level set values on mesh nodes close to the boundary.Moreover, we will also be able to deal with topological perturbations and treat shape and topological updates in a unified way by a discretized version of (19), called the numerical topological-shape derivative.

Numerical setting
In this section we consider the discretization of (2).Let T be a given finite element mesh covering D with M nodes {x k } M k=1 and N triangular elements {τ l } N l=1 .We introduce the index set element indices of elements τ l where x k is a node of τ l , Moreover, is the index set of all node indices of nodes x k in τl .Furthermore, we introduce the one-ring of a node These sets are illustrated in Figure 1.Let P 1 = {a + bx 1 + cx 2 : a, b, c ∈ R} denote the space of affine linear polynomials in two space dimensions and S 1 h (D) the space of piecewise affine linear and globally continuous functions on T , with the hat basis functions ϕ i ∈ S 1 h (D) which satisfy ϕ i (x j ) = δ ij , i, j = 1, . . ., M .The discretization of problem (2) leads to the discretized optimization problem min with the solution vector u ∈ R M and A = M + K. Here, the mass matrices M, M ∈ R M ×M , the stiffness matrix K ∈ R M ×M , and the right-hand-side vector f ∈ R M depend on the shape Ω and are given by On the reference element For an element τ l ∈ T , we denote the global vertex indices of its three vertices by l 1 , l 2 , l 3 and assume them to be numbered in counter-clockwise orientation.Then, the respective local finite element matrices and the local right-hand-side vector for element τ l are given by and the mapping Φ l between τ R and τ l and its Jacobian J l are given by

Numerical topological-shape derivative
Given the discretization introduced in Section 3, in contrast to the continuous topological-shape derivative, the numerical topological-shape derivative is only defined at the nodes of the finite element mesh T .For a given piecewise linear level set function φ ∈ S 1 h (D) let Ω(φ) be defined by (18a)-(18c) and is the finite element function corresponding to the solution of (23b).Note that, in this section, Ω(φ) is polygonal since φ ∈ S 1 h (D).The topological-shape derivative at node x k ∈ T is defined by Here, given φ ∈ S 1 h (D), the respective sets are defined by Whenever, the level set function φ is clear from the context, we will drop the argument and write T − , T + , S for brevity.Furthermore, is the positive discrete topological perturbation operator defined by its action Shape derivative: Initial domain in green, perturbed domain in orange.Interface is moved to the left.
Topological derivative: Shape derivative at a topologically perturbed level-set function (dashed blue) Finally, the discrete shape perturbation operator S k,ε : Remark 4.1.Note that the discrete perturbation operators defined above only change the nodal value of the finite element function

Computation of the numerical topological-shape derivative for the area cost functional
Before we compute the numerical topological-shape derivative (25) for the full model problem (2), we consider the case of the pure volume cost function and neglect the PDE constraint, i.e., we set c 1 = 1, c 2 = 0 in (1).
For that purpose, we investigate for , S k,ε }.We note that for the computation of (30) only those "cut elements" are relevant which have a node x k , i.e.
where H(x) denotes the Heaviside step function and is the set of all indices of elements adjacent to x k which are intersected by the perturbed interface.Note that, for ε > 0 small enough, C k does not depend on the concrete value of ε.For an element + 1 The nodal values of the level-set functions are indicated by −, +.The interface is drawn in red.
τ l with l ∈ I x k we denote the three vertices in counter-clockwise orientation by x l1 , x l2 , x l3 and assume that x k = x l1 .Moreover we denote φ lj := φ(x lj ) and φlj := O k,ε φ(x lj ) for j = 1, 2, 3 and small enough ε.In the following, we will be interested in with δ k,ε a l defined in (30).We consider six different sets (see Figure 3 for an illustration) such that with a direct sum on the right hand side.We can thus split the sum in (31) into six parts, x k in the case of topological perturbations.with Therefore, For 4a) and we have to consider (35) with φ l1 = 0. Thus, we obtain and conclude for this case Moreover, for x k (see Figure 4b) and we have for which leads to Therefore, we obtain for this case Finally, for x k ∈ S and O k,ε = S k,ε we deduce from (35) lim ε 0 and from (38) lim ε 0 Configuration B We note that Configuration B can only occur for the case and Thus, lim ε 0 For the case l ∈ I B− x k we have and obtain lim ε 0 Configuration C Analogously as for Configuration B, we note that Configuration C can only occur for the case x k ∈ S and O k,ε = S k,ε .We have lim ε 0 and lim ε 0 Summarizing, we have shown the following result.
Theorem 4.2.For x k ∈ T − we have For x k ∈ T + we have For x k ∈ S we have Remark 4.3.The corresponding computations for the denominators in (25), i.e. |Ω(O k,ε φ) Ω(φ)|, are closely related to the computations presented in this section for (30).Denoting where d k ãl = |d k a l | with the formulas for d k a l given in (46)-(48).It is obvious that, for x k ∈ T − , d k a l < 0 and for x k ∈ T + , d k a l > 0.Moreover, note that, for x k ∈ S, it holds d k a l < 0 for all l ∈ C k .This yields that

Computation of the numerical topological-shape derivative via Lagrangian framework
Next, we consider the computation of the numerical topological-shape derivative of an optimization problem that is constrained by a discretized PDE.For that purpose, we set c 1 = 0 in (1) and J(φ, u) := g(Ω(φ), u).The discretized problem reads and we are interested in the sensitivity of J when the level set function φ representing the geometry is replaced by a perturbed level set function φ ε = O k,ε φ.The perturbed Lagrangian for (51)-( 52) with respect to a perturbation of φ reads where we use the abbreviated notation A ε := A φε , and f ε := f φε .Moreover, for ε ≥ 0, we define the perturbed state u ε as the solution to i.e. u ε is the solution to and the (unperturbed) adjoint state p as the solution to for the state u given, i.e. p solves Note that we use the notation u for u ε=0 .The numerical topological-shape derivative at the node x k can be computed as the limit With the help of the Lagrangian (53), we can rewrite the right hand side as where we used that u ε solves (54) for ε ≥ 0. Following the approach used in [17], we use the fundamental theorem of calculus as well as (55) to rewrite this as Thus the numerical topological-shape derivative can be obtained as the sum of three limits, where , where Mε := Mφε represents the matrix M defined in (24) with Ω = Ω(φ ε ), we get Moreover, and Lemma 4.5.There exist constants c > 0, ε > 0 such that for all ε ∈ (0, ε) Here, o = 1 in the case of a shape perturbation and o = 2 in the case of a topological perturbation.
Proof.Subtracting (54) for ε = 0 from the same equation with ε > 0, we get and thus, by the ellipticity of the bilinear form corresponding to A ε and the triangle inequality, there is a constant c > 0 such that for all ε > 0 small enough For the difference between the perturbed and unperturbed right hand sides we have The result follows from the boundedness of |ϕ i (x)| and |∇ϕ i (x)| together with (49) which implies the existence of c > 0 such that |Ω(φ ε ) Ω(φ)| ≤ cε o (cf.Remark 4.3).
From Lemma 4.5 it follows that the terms R 1 (u, p) and R 2 (u, p) vanish.We remark that this is in contrast to the continuous case, where asymptotic analysis shows that R 2 does not vanish.We will address this issue in more detail in Section 5. Thus, in the discrete setting we obtain where with o = 1 for x k ∈ S and o = 2 for x k ∈ T − ∪ T + .To obtain (62), we divided both numerator and denominator of (60) by ε o and used that the limit of the quotient coincides with the quotient of the limits provided both limits exist and the limit in the denominator does not vanish.Next we state the numerical topological-shape derivative of problem (2) for arbitrary constant weights c 1 , c 2 ≥ 0 in the cost function (1).
Theorem 4.6 (Numerical topological-shape derivative).For l = 1, . . ., N , let u l = [u l1 , u l2 , u l3 ] and p l = [p l1 , p l2 , p l3 ] be the nodal values for element τ l of the solution and the adjoint, and Moreover, u k = u(x k ), p k = p(x k ), and ûk = û(x k ).For x k ∈ T − the numerical topological derivative reads whereas for x k ∈ T + we have For x k ∈ S the numerical shape derivative reads where the entries of the element matrix d k m I l and of the element vector d k f I l are dependent on the local cut situation (cases I = A ± , B ± , C ± ) and are given in Appendix A.2.The values d k ã can be computed by (48) considering Remark 4.3.
Proof.We evaluate (62) for . Thus, o = 2 in (63).We note that We have for and with (37) we obtain Due to and conclude Furthermore, with it follows that Analogously to (69) we have and obtain In the present situation, d k ã is given by the absolute value of (37) (see also Remark 4.3), By inserting (68), (70), (72), (74), and (75) in (62), together with Corollary 4.4, we obtain the sought expression (64).Formula (65) can be obtained in an analogous way as (64).The formula in (66) follows directly from (62) together with Corollary 4.4.The values of d k m I l and d k f I l for all possible cut situations I ∈ {A + , A − , B + , B − , C + , C − } are given in Appendix A.2 and were computed using symbolic computer algebra tools.A mathematical derivation of these terms is omitted here for brevity.

Connection between continuous and discrete sensitivities
The topological-shape derivative introduced in ( 25) and computed for model problem (2) in Theorem 4.6 represents a sensitivity of the discretized problem (23).In this section, we draw some comparisons with the classical topological and shape derivatives defined on the continuous level before discretization.While the purpose of this paper is to follow the idea discretize-then-differentiate, we consider the other way here for comparison reasons.

Connections between continuous and discrete topological derivative
For comparison, we also illustrate the derivation of the continuous topological derivative according to (4) for problem (2).We use the same Lagrangian framework as introduced in Section 4.2, see also [17].Given a shape Ω ∈ A, a point z ∈ D \ ∂Ω, an inclusion shape ω ⊂ R d with 0 ∈ ω and ε ≥ 0, we define the inclusion ω ε = z + εω and the perturbed Lagrangian For simplicity, we only consider the latter case, i.e. z ∈ D \ Ω in the sequel.Noting that u ε , ε ≥ 0, is the solution to the perturbed state equation with parameter ε, the topological derivative can also be written as with the solution to the unperturbed adjoint state equation p.As in Section 4.2, this leads to the topological derivative consisting of the three terms where provided that these limits exist.It is straightforwardly checked that for z For the term R 2 (u, p), we obtain We have a closer look at the diffusion term In the continuous setting, we now define Next, one can show the weak convergence ∇K ε ∇K for K ∈ ḂL(R 2 ) being defined as the solution to the exterior problem Assuming continuity of ∇p around the point of perturbation z, it follows that It can be shown that the other terms in (77) vanish and thus R 2 (u, p) = R λ 2 (u, p).Finally, it follows from the analysis in [17,Sec. 5 Comparing the topological derivative formula obtained here with the sensitivity for interior nodes x k ∈ T − ∪ T + obtained in Section 4, we see that the term corresponding to R λ 2 (u, p), i.e., the term lim , vanishes in the discrete setting.This can be seen as follows: For u h ε , ε ≥ 0, and p h ∈ V h , we have the expansion in the finite element basis If we now plug in these discretized functions into (78) and consider a fixed mesh size h, we get on the other hand where we used the continuity of ε → u ε according to Lemma 4.5.Note that, since the mesh and the basis functions are assumed to be fixed and independent of ε, unlike in the continuous setting, here the continuity of the normal flux of the discrete solution u h ε across the interface ∂ω ε is not preserved.We mention that, when using an extended discretization technique that accounts for an accurate resolution of the material interface (e.g.XFEM [23] or CutFEM [10]), the corresponding discrete sensitivities would include a term corresponding to R λ 2 (u, p).

Connection between continuous and discrete shape derivative
The continuous shape derivative dg(Ω)(V ) for a PDE-constrained shape optimization problem given a shape Ω ∈ A and a smooth vector field V can also be obtained via a Lagrangian approach.For our problem (2), it can be obtained as dg(Ω)(V ) = ∂ t G(0, u, p) with , see [27] for a detailed description.In the continuous setting, the shape derivative reads in its volume form 1 and S Ω 0 given in ( 11) and ( 12), respectively.Under certain smoothness assumptions, it can be transformed into the Hadamard or boundary form )n • n given by where L λ is given by Here, ∇u in , ∇p in and ∇u out , ∇p out denote the restrictions of the gradients to Ω and D \ Ω, respectively, and n denotes the unit normal vector pointing out of Ω.Note that, when using a finite element discretization which does not resolve the interface such that the gradients of the discretized state and adjoint variable are continuous, i.e. ∇u h,in = ∇u h,out and ∇p h,in = ∇p h,out , (83) becomes We now discretize the continuous shape derivative formula (81) for the vector field V (k) that is obtained from the perturbation of the level set function φ in (only) node x k , k ∈ S fixed.For that purpose we fix φ ∈ S 1 h (D) and the corresponding domain Ω(φ).Note that we consider V (k) to be supported only on the discretized material interface ∂Ω(φ) ∩ D. We begin with the case of the pure volume cost function by setting c 2 = 0. Proposition 5.2.Let c 2 = 0 and x k ∈ S fixed.Let V (k) the vector field that corresponds to a perturbation of the value of φ at position x k .Then where d k a l is as in (48).
Proof.For c 2 = 0 we also have p = 0 and thus L = c 1 , i.e., we are in the case of pure volume minimization.From (81) we know that dg(Ω(φ))(V (k) ) = c 1 ∂Ω(φ) V (k) • n dS x .First of all, we note that the vector field V (k) corresponding to a perturbation of φ at node x k is only nonzero in elements τ l for l ∈ C k with C k as defined in (32).Thus, the shape derivative reduces to We compute the vector field V (k) and normal vector n explicitly depending on the cut situation.
Recall the sets x k introduced in (34), see also Figure 3.Given two points p and q and their respective level set values a and b of different sign, ab < 0, we denote the root of the linear interpolating function by and note that d da y(p, q, a, b) = − b (b−a) 2 (q − p).We begin with Configuration A. For an element index l ∈ I x k , we denote the vertices of element τ l by x l1 , x l2 , x l3 and assume their enumeration in counter-clockwise order with x k = x l1 .The corresponding values of the given level set function φ are denoted by φ l1 , φ l2 , φ l3 , respectively.We parametrize the line connecting the two roots of the perturbed level set function along the edges by and obtain the vector field corresponding to the perturbation of φ k = φ l1 along the line τ l ∩ ∂Ω(φ) as Introducing the notation d ki,kj := |y(x l k , x li , φ l k , φ li ) − y(x l k , x lj , φ l k , φ lj )| for the length of the interface in element τ l , the normed tangential vector along τ l ∩ ∂Ω(φ) and the normed normal vector pointing out of Ω(φ) are given by tA where R denotes a 90 degree counter-clockwise rotation matrix, Finally, by elementary computation we obtain for l ∈ I and the same formula with a different sign for l ∈ I A− x k .Proceeding analogously, we obtain for l ∈ I

B+
x k and l ∈ I and further respectively.Again, the formulas for l ∈ I B− x k and l ∈ I C− x k just differ by a different sign.Finally, comparing the computed values with the formulas of d k a l (48) yields the claimed result.
In view of Proposition 5.2, the definition in (13) and Remark 4.3, we see that, in the case c 2 = 0, it holds which is in alignment with the first term of the formula in (66).Next, we consider the general PDE-constrained case where c 2 > 0.
Proposition 5.3.Let c 1 = 0 and x k ∈ S fixed.Let V (k) the vector field that corresponds to a perturbation of the value of φ at position x k .Then where we use the same notation as in Theorem 4.6.In particular, d k m I l and d k f I l depend on the cut situation, I ∈ {A + , A − , B + , B − , C + , C − } and are given explicitly in Appendix A.2.
Proof.Let an element index l ∈ C k fixed and u l = [u l1 , u l2 , u l3 ] , p l = [p l1 , p l2 , p l3 ] contain the nodal values of the finite element functions u h and p h corresponding to the three vertices x l1 , x l2 , x l3 of element l, respectively.Also here, the ordering is in counter-clockwise direction starting with x l1 = x k .We compute the shape derivative (81) with L given in (82) after discretization (i.e. after replacing the functions u, p, û by finite element approximations u h , p h , ûh ).In particular, the term L λ is approximated by (84).Depending on how the material interface ∂Ω(φ) cuts through element τ l , the normal component of the vector field V (k) along the line τ l ∩ ∂Ω(φ) is given in (85)-(87).For and I ∈ {A + , A − , B + , B − , C + , C − }, it can be seen by elementary yet tedious calculations that with d k m I l and d k f I l as given in Appendix A.2. Examplarily, we illustrate the calculation for the second of these terms for the cut situation I = A + , see Figure 3(a).Let u l,12 and u l,13 denote the values of the linear function u h | τ l at the intersection of the interface ∂Ω(φ) with the edges that connect the point x l1 with x l2 and x l1 with x l3 , respectively.Note the relations u l,12 = and u l,13 = . Analogously we define the values p l,12 and p l,13 .The function u h along the line τ l ∩ ∂Ω(φ) can now be written as ûh (s) = u l,12 + s(u l,13 − u l,12 ), s ∈ [0, 1] and we get The last equality is obtained by plugging in (85) and straightforward (yet tedious) calculation.Finally, since u h and p h are linear and the normal vector is constant on τ l ∩ ∂Ω(φ), we see that L λ h is constant and, using Proposition 5.2, we obtain Combining the findings of Propositions 5.2 and 5.3 and dividing by d k ã (defined in 4.3), we obtain the resulting formula for the alternative definition of the shape derivative as defined in ( 13)

Complex step derivative test
In order to overcome this drawback of the finite difference test, we next consider a test based on the complex step (CS) derivative [22].For that purpose, using Remark 4.3, let us first rewrite (91) as where o = 1 if x k ∈ S and o = 2 if x k ∈ T − ∪ T + .Moreover, assuming a higher order expansion of the form with some higher order sensitivities d 2 J (φ)(x k ), d 3 J (φ)(x k ) and assuming that (94) also holds for complex-valued ε, we can follow the idea of the complex step derivative [22]: Setting ε = ih in (94) with h > 0 and i the complex unit yields in the case o = 1 where x k ∈ S, and in the case o = 2 where x k ∈ T − ∪ T + .This means α1 α2 α 1 α 2 λ 1 λ 2 f 1 f 2 1 0.9 1 0.2 1 0.6 1 0.5 Table 1: Problem parameters for the numerical experiment Analogously to (92), we define the summed errors e CS S (h) and e CS T (h) by just replacing δJ ε (φ)(x k ) by δJ CS h (φ)(x k ) defined above.Figure 5(b) shows the errors e CS S and e CS T for a positive, decreasing sequence of h where we observe quadratic decay for both errors.While the error e CS S corresponding to the shape nodes x k ∈ S decays to machine precision, the error e CS T corresponding to the interior nodes x k ∈ T − ∪ T + deteriorates at some point due to the cancellation error ocurring when subtracting J (φ) from J (O k,ih φ) in (96).

Test based on hyper-dual numbers
In order to overcome this cancellation error also for the case of x k ∈ T − ∪ T + , we resort to hyperdual (HD) numbers as introduced in [14].Here, the idea is to consider numbers with three non-real components denoted by Assuming that expansion (94) holds up to order o + 1 also for such hyper-dual values of ε, we can set ε = hE 1 + hE 2 + 0E 1 E 2 for some h > 0. For o = 1, considering only the first non-real part (i.e., the E 1 -part) and exploiting that E 2 1 = 0, we obtain the equality for x k ∈ S. Similarly, with the same choice of ε, by considering only the E 1 E 2 -part of the expansion and exploiting for x k ∈ T − ∪ T + .In this case, the corresponding summed errors e HD S (h) and e HD T (h) vanish for arbitrary h ∈ R.This is also observed numerically since both (97) and (98) suffer neither from a truncation nor a cancellation error.Figure 5(c) shows that the the obtained results agree up to machine precision with the derivatives obtained by (64), (65), and the respective formula for the shape derivative (66).

Application of optimization algorithm to model problem
Finally we show the use of the numerical topological-shape derivative computed in Section 4 within a level-set based topology optimization algorithm.We first state the precise model problem, before introducing the algorithm and showing numerical results.

Problem setting
We consider the unit square D = [0, 1] 2 and minimize the objective function (1) with c 1 = 0 and c 2 = 1 subject to the PDE constraint in (2).The chosen problem parameters are shown in Table 1.We consider a mixed Dirichlet-Neumann problem by choosing This yields that Ω are two (approximated) circles with radii 0.2 and 0.1 respectively, see Figure 6.

Optimization algorithm
The optimization algorithm we use to solve the problem introduced in Section 6.2.1 is inspired by [4].
Definition 6.1.We say a level set function φ ∈ S 1 h (D) is locally optimal for the problem described by J if We introduce the generalized numerical topological-shape derivative G φ ∈ S 1 h (D) with With this definition, we immediately get the following optimality condition: Then φ is locally optimal in the sense of Definition 6.1.
Remark 6.6.In practice it turned out to be advantageous to include a smoothing step of the level set function.Thus, we chose the following update strategy: We first set , with the same notation as above before smoothing the level set function in T − (ψ) ∪ T + (ψ) by Finally, the level-set function is normalized and the next iterate is given by

Numerical results
As an initial design for the optimisation, we take the empty set, Ω = ∅.This is realized by choosing φ 0 = 1/ 1 L2(D) as the initial level set function.We use the algorithm outlined in Section 6.2.2 to update this level set function.We terminated the algorithm after the fixed number of 800 iterations.
The final as well as some intermediate configurations are illustrated in  for the five different levels of discretization.We observe that in all cases the two circles are recovered with high accuracy.In Figure 12 the evolution of the objective function as well as of the norm of the generalized numerical topologicalshape derivative is plotted.We observe that objective function decreases fast and after 800 iterations a reduction by a factor of approximately 10 −5 − 10 −8 could be achieved.Moreover, we observe that the norm of the topological-shape derivative decreases continuously, more and more approaching the optimality condition (102).

Conclusions
In this work we presented a new sensitivity concept, called the topological-shape derivative which is based on a level set representation of a domain.This approach allows for a unified sensitivity analysis for shape and topological perturbations, which we carried out for a discretized PDE-constrained design optimization problem in two space dimensions.For the discretization we used a standard first order finite element method which does not account for the interface position in the approximation

Let 2 D
D be a given, fixed, open and bounded hold-all domain and Ω ⊂ D an open and measurable subset.Let the boundary of D be decomposed into Γ D , Γ N ⊂ D with Γ D ∪Γ N = ∂D and Γ D ∩Γ N = ∅.In the present paper, we consider a topology optimisation problem with a tracking type cost function g(Ω, u) = c 1 |Ω| + c αΩ |u − û| 2 dx (1) where û ∈ H 1 (D) is a given desired state, and c 1 , c 2 ∈ R are given constants.The continuous topology optimization problem reads min Ω∈A g(Ω, u),

Figure 1 :
Figure 1: Illustration of the sets I x k , J τ l , and R x k .

Figure 2 :
Figure 2: Illustration of the shape derivative and the topological derivative.

Figure 4 :
Figure 4: Illustration of I A+ x k and I A−

Figure 5 :
Figure 5: (a) Results of the finite difference test.(b) Results obtained with the complex step derivative.(c) Results obtained with hyper-dual numbers.

Figure 6 :
Figure 6: The different meshes and corresponding sought shapes used in the numerical experiments.

Figure 7 :
Figure 7: Evolution of the level-set function for the 145 nodes mesh

Figure 8 :
Figure 8: Evolution of the level-set function for the 545 nodes mesh

Figure 9 :
Figure 9: Evolution of the level-set function for the 2113 nodes mesh

Figure 10 :
Figure 10: Evolution of the level-set function for the 8321 nodes mesh

Figure 11 :
Figure 11: Evolution of the level-set function for the 33025 nodes mesh