A variational formulation for computing shape derivatives of geometric constraints along rays

. In the formulation of shape optimization problems, multiple geometric constraint functionals involve the signed distance function to the optimized shape Ω. The numerical evaluation of their shape derivatives requires to integrate some quantities along the normal rays to Ω, a challenging operation to implement, which is usually achieved thanks to the method of characteristics. The goal of the present paper is to propose an alternative, variational approach for this purpose. Our method amounts, in full generality, to compute integral quantities along the characteristic curves of a given velocity ﬁeld without requiring the explicit knowledge of these curves on the spatial discretization; it rather relies on a variational problem which can be solved conveniently by the ﬁnite element method. The well-posedness of this problem is established thanks to a detailed analysis of weighted graph spaces of the advection operator β · ∇ associated to a C 1 velocity ﬁeld β . One novelty of our approach is the ability to handle velocity ﬁelds with possibly unbounded divergence: we do not assume div( β ) ∈ L ∞ . Our working assumptions are fulﬁlled in the context of shape optimization of C 2 domains Ω, where the velocity ﬁeld β = ∇ d Ω is an extension of the unit outward normal vector to the optimized shape. The eﬃciency of our variational method with respect to the direct integration of numerical quantities along rays is evaluated on several numerical examples. Classical albeit important implementation issues such as the calculation of a shape’s curvature and the detection of its skeleton are discussed. Finally, we demonstrate the convenience and potential of our method when it comes to enforcing maximum and minimum thickness constraints in structural shape optimization.

The recent achievements of shape and topology optimization techniques in predicting very efficient designs beyond the reach of the intuition of engineers have raised a tremendous enthusiasm in industry. One key issue about the practical use of these techniques is that they often lead to shapes which are too complicated to be assembled by traditional milling or casting processes, to the point that one of the greatest challenges faced by shape and topology optimization is about the understanding and the implementation of manufacturing constraints. Such constraints generally bring into play the geometry of shapes; for instance, minimum thickness constraints, maximum thickness constraints, minimum distance between members, casting constraints have been considered; see [62,16,38,50] or the recent survey [46]. Among the variety of methods proposed in the literature to enforce these constraints, a body of work has recently been devoted to formulate them by means of integral functionals involving the signed distance function d Ω to the optimized domain Ω, that is, the distance function to ∂Ω with a negative (resp. positive) sign inside (resp. outside) Ω; see Section 3.1 below, and [5,6,21,22,50] for the formulation of geometric constraints based on this concept. More precisely, the shape Ω is sought among all possible subsets of a fixed 'hold-all' domain D ⊂ R n (n = 2, 3 in applications) as the solution of a constrained minimization problem of the form: where J(Ω) is a measure of the performance of Ω (e.g. its elastic compliance), j : R → R is a given, smooth function, and P (Ω) is a geometric constraint. The above framework is quite appealing insofar as the signed distance function d Ω naturally lends itself to clear mathematical formulations of maximum, minimum thickness constraint functionals (and the other aforementioned geometric criteria), thus making (1.1) amenable to an efficient numerical treatment by means of standard (e.g. steepest-descent) optimization algorithms. However, several technical stages in its numerical implementation pose difficulties; the leading motivation of the present article is to introduce a new variational method that make these substantially simpler to implement. Let us provide a little more details about the main numerical issues raised by the implementation of (1.1), while staying at an explanatory level; see Sections 3.1 and 4 below for full details. The use of traditional optimization methods for the program (1.1) relies (in particular) on the knowledge of the shape derivative of the considered constraint functional Ω → P (Ω). Throughout this article, this notion is understood in the framework of the Hadamard method: the shape derivative corresponds to the Fréchet derivative of θ → P (Ω θ ) at θ = 0, where local variations of Ω are considered under the form Ω θ = (Id + θ)(Ω), using vector fields θ ∈ W 1,∞ (D, R n ) with sufficiently small norm (see Section 4 hereafter). It was shown in our previous work [2] that the shape derivative of P (Ω) as in (1.1) has the form P (Ω)(θ) = ∂Ω u θ · n dy, (1.2) where dy refers to the surface measure on ∂Ω, and the function u ∈ L 1 (∂Ω) is given by the formula In (1.3), n and κ i (y) stand for the unit normal vector to ∂Ω pointing outward Ω and the (n − 1) principal curvatures of ∂Ω respectively; the set ray(y) is the normal ray that originates from y ∈ ∂Ω and that stops at either the boundary ∂D of the hold-all domain, or at the skeleton (or medial axis) Σ of Ω; see Section 3.1 for a short review of these notions and Figure 1a for an illustration. The normal velocity θ = −un is then exploited by the optimization algorithm as a constraint gradient to build a descent update for the shape Ω [1]. The numerical component of this program raises the need to calculate the quantity u : ∂Ω → R defined in (1.3).
In principle, the formula (1.3) featuring integration along the normal rays to ∂Ω can be implemented efficiently. In particular, its numerical evaluations at several points y ∈ ∂Ω can be performed in parallel. However, the practical implementation when all the functions involved in the integrand of (1.3) are discretized on a computational mesh is not trivial: it requires the computation of (i) the ray trajectories on the discrete mesh (for instance by some variant of the method of characteristics), as exemplified on Figure 1b, and (ii) the (a) The shape optimization setting: shape Ω ⊂ D, skeleton Σ with normal rays (ray(y)) y∈∂Ω , local coordinate change η, outward normal vector n, center of curvature C, orthogonal projection p ∂Ω (x) of a point x ∈ D onto ∂Ω.
(b) Integration along rays requires the detection of the skeleton of the shape and the computation of the triangle path at the mesh level, that is the intersection points of each triangle with the ray. principal curvatures κ i of the boundary ∂Ω from the datum of its numerical definition. As we shall illustrate in our numerical investigations, the accurate calculation of these quantities -even when ∂Ω is endowed with a discretization as an explicit triangulated surface-is not trivial and is still an active research direction, already in the case of two space dimensions, and especially in the three-dimensional context [49,55,31].
The objective of this paper is to introduce a robust and easy to implement method that allows to compute integral quantities such as (1.3) by solving a variational problem for advection-like operators which alleviates the need for resorting to direct integration along rays. This yields an efficient numerical method to calculate shape derivatives of geometric constraints, which is very simple to implement using a standard finite element software, and which does not raise any additional algorithmic difficulty in 3-d than in 2-d.
Our numerical method for the calculation of (1.3) is based on the following variational problem, set over a suitable Hilbert space V ω : Find u ∈ V ω such that ∀v ∈ V ω , ∂Ω uvds + D\Σ ω(∇d Ω · ∇u)(∇d Ω · ∇v)dx = − D\Σ j (d Ω )vdx, (1.4) where ω > 0 is a rather arbitrary weight function (over which relevant assumptions shall be stated later on), and Σ is the closure of the skeleton set Σ of Ω (see Figure 1a and Section 3.1 below). Indeed, our key result is to show that (1.4) is well-posed, and that the trace of its unique solution u on ∂Ω is exactly (1.3).
A formal insight that allows to see this is the following: assume that for any v 0 ∈ C 0 (∂Ω), the function v satisfying v = v 0 on ∂Ω and taking constant values along the rays normal to ∂Ω (i.e. ∇d Ω · ∇v = 0 in D\Σ) belongs to the space V ω . Using v as a test function in the variational formulation (1.4) then yields: (1 + κ i (y)d Ω (z))ds   v 0 (y)dy, (1.5) where the last equality follows from a change of variables allowing to rewrite integration on D \ Σ as a nested integration on points y ∈ ∂Ω and on elements z in ray(y) (see (2.9) and (2.8) below). Since v 0 ∈ C 0 (∂Ω) can be chosen arbitrarily in (1.5), the identity (1.3) follows. These considerations are made rigorous in Section 2.5.
The variational formulation (1.4) makes it possible to compute integrals (1.3) along the normal rays to ∂Ω without the need to calculate these rays or the curvatures κ i explicitly on a discretization (i.e. a mesh) of the ambient space. This feature is especially convenient for shape optimization applications relying on the level set method, as described in Section 4; there, the gradient of the signed distance function ∇d Ω is easy to calculate on an unstructured mesh of the considered hold-all domain D from a P 1 approximation of d Ω . The variational formulation (1.4) can then be readily implemented in a finite element framework, even if the boundary Γ = ∂Ω is not meshed explicitly. Our method requires only the identification of the skeleton Σ. Although this task is notoriously difficult to carry out with a high accuracy, only a rough estimate of the location of Σ is needed for our purpose, upon a judicious choice of the weight ω in (1.4); see the numerical examples in Section 3. Last, the previous arguments work identically when ∇d Ω is replaced with an arbitrary C 1 vector field β: we obtain that a variational formulation analogous to (1.4) (given in (2.5)) allows to compute integrals quantities along the characteristics curves of β, which is subject to offer wider applications than shape optimization.
With these perspectives in mind, the present article is organized as follows. In Section 2, we carefully discuss the mathematical setting that guarantees the existence and uniqueness of a solution to the variational problem (1.4) (in fact its generalization to arbitrary vector fields β), and the justification of the key identity (1.3) satisfied by the trace of its solution. For that purpose, we provide a variational theory for the advection operator β · ∇ associated to arbitrary C 1 vector fields β on the weighted graph space V ω ; in particular, we examine the existence of traces on Γ and we prove an adapted Poincaré inequality for functions v ∈ V ω . Let us stress that related works about these matters (e.g. [29,32]) usually require less regularity on the vector fields β and do not rely on the existence of a flow η (see (2.1) below), but at the cost of rather strong boundedness assumptions on the divergence of β (typically div(β) ∈ L ∞ (D)), which do not hold for our shape optimization applications (2.7) and (2.8). Indeed, in the latter situation where β = ∇d Ω , the divergence div(∇d Ω ) typically blows up near Σ; hence the need for our different approach.
In Section 3, we investigate the numerical accuracy of our variational method for calculating integrals along characteristic curves in the shape optimization setting (2.7) where β = ∇d Ω . After a short review of the properties of the signed distance function, we compare the direct numerical integration along rays with the use of our variational method on several numerical examples where the value of (1.3) is analytically known. We also consider "practical" cases where some of the regularity assumptions imposed by our framework are not fulfilled. Most importantly, we illustrate how the selection of an appropriate weight ω in the variational formulation (1.4) allows to deal with the presence of cracks in the working domain, for instance the skeleton Σ in shape optimization when it is not explicitly meshed.
Eventually, in Section 4, we turn to practical shape optimization applications, in the particular case of two space dimensions, in order to demonstrate the simplicity and effectiveness of our method. Let us emphasize that the variational approach (2.5) is equally simple to implement in any space dimension while the 3-d implementation of geometric integration along rays would require much more efforts than in 2-d. Elaborating on our previous works [3,6,22] we demonstrate that the method proposed in this paper allows for a convenient and efficient implementation of maximum and minimum thickness constraints in structural design.
2. Weighted graph space of the advection operator β · ∇ for velocity fields of class C 1 In this mathematical section, we consider a slightly more general variational problem than (1.4): ∇d Ω is replaced by a rather arbitrary vector field β. This setting and some technical assumptions are described in Section 2.1. In order to obtain its well-posedness and the trace identity (1.3), we introduce in Section 2.2 suitable functional spaces V ω , in which the directional derivative β · ∇ naturally makes sense. Then, we show in Section 2.3 that under suitable hypotheses, C 1 functions are dense in V ω . In Section 2.4, we establish a trace theorem for functions in V ω and we provide a Poincaré-type inequality. Finally, we prove in Section 2.5 that the generalized version of the variational problem (1.4) is well-posed on V ω and we obtain an identity generalizing (1.3) in this context involving arbitrary fields β.

Preliminaries, notations and assumptions
Let U ⊂ R n , Γ ⊂ U and β : U → R n be respectively a (possibly non smooth) bounded open set, a hypersurface and a vector field of class C 1 . We do not require Γ to be a compact manifold ( Γ may differ from its closure Γ); in this case, we require Γ to be a manifold with boundary, which in particular prevents Γ from showing spiralling patterns near its ends-an assumption which is needed for technical reasons (see e.g. the proof of Lemma 3 below). Two examples of admissible settings are represented on Figure 2.
We assume that Γ is a Poincaré section or a stream surface for β, meaning that for any y ∈ Γ, there is a unique characteristic curve s → η(y, s) passing through y = η(y, 0) at time s = 0, and that lives in the domain U on some maximal interval s ∈ (ζ − (y), ζ + (y)). In other words, for any y ∈ Γ, (ζ − (y), ζ + (y)) is the maximum existence interval such that the solution s → η(y, s) of the ordinary differential equation    d ds η(y, s) = β(η(y, s)), η(y, 0) = y, remains in the domain U (note that by definition, ζ − ≤ 0 ≤ ζ + ). We assume that ζ − and ζ + are continuous, bounded functions on Γ, satisfying the following separation condition: The vector field β is required to be U -filling, in the sense that its related flow η realizes a C 1 diffeomorphism from the tensor product set W = {(y, s) | y ∈ Γ, s ∈ (ζ − (y), ζ + (y))}, (2.3) onto U (see [44], Chap. IV), where W is a manifold obtained as an open subset of the tensor product of Γ with the real line R (note that the "open" character of W comes from the continuity assumption on ζ − and ζ + ).
Finally, we denote by |Dη| the Jacobian of the local coordinate change η: where the Jacobian matrix of η reads ∇η = ∂ y η ∂ s η and ∂ y denotes the collection of derivatives with respect to the (n − 1) tangential coordinates of Γ.
Remark 1. Let us consider a particular case where U is a smooth bounded domain in R n , and Γ ⊂ ∂U is defined as the inlet boundary of the flow field β, i.e.
where n is the unit normal vector to ∂U , pointing outward U , the separation condition (2.2) exactly requires that the inflow and outflow boundaries be separated by a positive minimum distance, which is a rather classical assumption in the study of the advection operator, see e.g. Section 2.1.3 in [28]. In our case, assumption (2.2) is required in the proof of Proposition 2.
The main point of this section is to mathematically justify that given rather arbitrary function f : U → R and weight ω ∈ C 0 (U, R * + ), the trace of the solution u to the variational problem is given explicitly in terms of line integrals by the formula which yields a numerical method for computing (2.6) by solving (2.5). The shape optimization setting outlined in the introduction, which is exemplified on Figure 1a, reduces to the particular case where Ω ⊂ D are bounded Lipschitz domains of R n , and ζ − (y) and ζ + (y) are the distances at which ray(y) hits either the skeleton Σ of Ω or the boundary ∂D of the hold-all domain D. Indeed, in this situation, the local coordinate change η and its Jacobian |Dη| are explicitly given by (see [13] and Section 3.1): (1 + κ i (y)s), ∀y ∈ ∂Ω, ∀s ∈ (ζ − (y), ζ + (y)), (2.8) so that the function u of (1.3) coincides with the expression (2.6) for f = −j (d Ω ). Let us emphasize that, in this context, the open set U is not smooth because it features a "crack", namely the skeleton Σ (we call it a "cracked domain" in Figure 1a); in particular, U is not located on one side of its boundary. This "lack of smoothness" of U prevents from using many convenient results from functional analysis [59], and thus raises the need for several technical ingredients in the sequel, which otherwise would have been fairly classical.
Throughout the article, the considered measure on W is that induced by the standard product measure of the surface measure dy on Γ and of the Lebesgue measure ds on R. Thus, the space L 1 (W ) of integrable functions on W is defined as the space of measurable functions f : W → R such that Γ ζ+(y) ζ−(y) |f (y, s)|dsdy < +∞. Integration of f over W is then defined by: Under these notations, the classical change of variable formula between manifolds reads (see [44], Chap. XVI): where the Jacobian |Dη| is defined by (2.4).
Note that, above and in the sequel, the notations α −1 (y, s) = 1/α(y, s) and ω −1 (y, s) = 1/ω(y, s) are used for the inverses of the scalar weights α and ω respectively, whereas the notation η −1 : U → W shall stand for the reciprocal mapping of the diffeomorphism η. Hypothesis (H1) is essentially a regularity assumption, positivity being no surprise for ω to be an admissible weight in the variational formulation (2.5). Notice that no assumption is made about the behavior of ω near the boundary of U ; in particular, ω(x) may tend to 0 as x approaches ∂U . Hypothesis (H2) is an upper-boundedness assumption for α. The weights that we are going to consider in our practical applications in Section 3.3 will satisfy (H1) and (H2) almost automatically. Finally, (H3) is roughly a monotonicity constraint for the decay of α towards 0 as s → ζ ± (y). In practice, we will rely on the following lemma which provides a simple monotonicity condition under which the condition (H3) is fulfilled, indicating that the class of weights satisfying (H3) is large enough. Lemma 1. Let α ∈ C 0 (W, R * + ) be a weight of the product form: ∀(y, s) ∈ W, α(y, s) = f (y, s)g(y, s) (2.12) where f and g are two positive functions on W satisfying: (i) There exist two constants g − , g + ∈ R such that for all (y, s) ∈ W , 0 < g − ≤ g(y, s) ≤ g + .

2.2.
Definition of the graph space V ω of the advection operator β · ∇ We recall that the assumption β ∈ C 1 (U, R n ) implies that div(β) belongs to L ∞ loc (U ). Definition 1 (Derivative along characteristic curves). Let v ∈ L 1 loc (U ) be a locally integrable function on U . The directional derivative β · ∇v ∈ D (U ) is the distribution on U defined by Remark 4. The definition (2.15) of β · ∇ mimics the integration by part formula that holds for functions v of class C 1 . Considering test functions φ ∈ C ∞ c (U ) allows to avoid imposing classical regularity requirements on the open domain U such as that being locally located on one side of its boundary [59].
Definition 2 (Graph space of the operator β · ∇). Let ω ∈ C 0 (U, R * + ) be a positive weight on U . The weighted graph space V ω of the advection operator β · ∇ is defined by: V ω is a Hilbert space when it is endowed with the norm Remark 5. The completeness of V ω follows from very standard arguments, see e.g. [60,28]. It is also easy to see from the assumption ω ∈ C 0 (U, R * + ) that the inclusion L 2 ω (U ) ⊂ L 1 loc (U ) holds. Let us then introduce corresponding notions of derivative with respect to the s variable and of graph space on the set W : Definition 3 (Graph space of the derivative ∂ s ). Let α ∈ C 0 (W, R * + ) be a positive weight on W . The weak derivative ∂ s v ∈ D (W ) of a function v ∈ L 1 loc (W ) is the distribution defined by The weighted graph space of the operator ∂ s is defined by: . Remark 6. The defining identity (2.15) of the weak derivative β · ∇v also holds for test functions φ ∈ C 1 c (U ), as follows from a standard density argument. Likewise, (2.17) holds for test functions φ which are only continuously differentiable with respect to the s variable.
The motivation behind the introduction of the spaces V α is that they allow to transfer the difficulty of studying the "curvilinear" directional derivative β · ∇ over the domain U , to that of the more standard, "flat" derivative ∂ s over W . The price to pay is the need to account for the weight of the Jacobian |Dη| resulting from the change of variables. This point is made precise by the following technical lemma.

Conversely, one proves in a similar way that if
Remark 7. In the following we prove a density result, a trace theorem and a Poincaré inequality in V α and we then obtain the direct counterparts of these results in the setting of the graph space V ω thanks to the above Lemma 2. There exists actually a wide literature about weighted Sobolev spaces such as V α [34,42,43,60], in which analogous results are proved for weights α in the so-called Mückenhoupt class A 2 (see e.g. [30] for a definition of the class A p ). Our working assumptions however are of a different nature: for instance, the hypothesis (H3) essentially requires that the inverse weight α −1 belong to the Mückenhoupt class A 1 .

Density of functions of class C 1 in the weighted space V ω
We now examine the density of C ∞ (W )∩ V α in V α , whence we shall infer the density of C 1 (U )∩V ω in V ω -the loss of regularity between both statements coming from the fact that the coordinate change η is only of class C 1 on W .
Our study classically involves mollifying functions; since the space V α of interest contains functions defined on the manifold W , a little treatment is in order. In particular, we shall need the so-called Ahlfors regularity of Γ (see [24]); this is the purpose of the next lemma, whose proof is outlined for the convenience of the reader. Here and throughout the article, B(y, h) ⊂ R n stands for the open ball with center y and radius h. The measure of a Lebesgue measurable set A ⊂ R n is denoted by |A|.
Lemma 3 (Area covered by extrinsic balls on an embedded manifold). Let Γ ⊂ R n be a C 1 hypersurface such that either Γ is compact (Γ = Γ) or Γ is a compact manifold with boundary. Then there exists h 0 > 0 and constants m > 0 and M > 0 depending only on Γ such that for any 0 < h < h 0 , Proof. Since Γ is compact, there exist finitely many open subset V i ⊂ R n , i = 1, ..., N , such that Γ ⊂ i V i , and for each i = 1, ..., N , there exists a local coordinate chart φ i : is a convex open subset, and R n−1 + is the upper half-space of R n−1 . Let h 0 > 0 be a Lebesgue's number associated with this cover, that is: Now, given y ∈ V i , one has for any 0 < h < h 0 : where i is the index supplied by (2.22) and 1 A denotes the characteristic function of a subset A ⊂ R n , and |Dφ i | is the Jacobian associated to the change of variables between manifolds induced by φ i .
Let σ n−1 (∇φ i (x)) ≤ σ 1 (∇φ i (x)) be respectively the smallest and the largest singular values of the n × (n − 1) Jacobian matrix ∇φ i (x), and let 0 < σ − ≤ σ + be defined by: Applying the Taylor formula to φ i yields: Therefore, possibly shrinking the value of the constant h 0 supplied by (2.22) and taking x 0 = φ −1 i (y), we obtain: Finally, denoting by |B n−1 | the volume of the unit ball in R n−1 , we obtain: which completes the proof.
In what follows, the hypersurface Γ is the one introduced by Section 2.1. In the next lemma, we construct the kernels ρ h and ξ h which shall be used for mollification purposes in W .

Lemma 4.
For any h > 0, there exist two positive, smooth functions ρ h ∈ C ∞ c (R n ×R n , R) and ξ h ∈ C ∞ c (R, R) satisfying the conditions: (i) For all points y ∈ R n , supp(ρ h (y, ·)) ⊂ B(y, h).
(iii) There exist constants C > 0 and h 0 > 0 depending only on Γ such that is not a function of (y − z). The conditions (i) and (ii) of the statement are obviously satisfied by (2.24). The condition (iii) is a consequence of Lemma 3, which implies: so that (2.23) holds with C = 2 n−1 M/m. Eventually, a function ξ h satisfying (iv) is constructed from any positive function ξ ∈ C ∞ c (R, R) with compact support inside [−1, 1] and unit integral over R, by setting ξ h = h −1 ξ(·/h).
Definition 4 (Mollification on the tensor product manifold W ). For h > 0, let ρ h and ξ h be two functions as in the statement of Lemma 4. For any u ∈ L 1 loc (W ), (y, s) ∈ W and h > 0 sufficiently small (depending on (y, s)), the mollification of u is the function (2.25) Note that for a given open subdomain W 1 ⊂⊂ W and because Γ is compact, there exists h 0 > 0, depending on W 1 , sufficiently small so that (2.25) makes sense for any (y, s) ∈ W 1 and 0 < h < h 0 .
Lemma 5. The following properties hold true: Proof.
(i) This results from the Lebesgue dominated convergence theorem and the smoothness of ρ h and ξ h .
(ii) This is again a consequence of the Lebesgue dominated theorem and of an integration by parts: (iii) For a given subset W 1 ⊂⊂ W , let > 0 and h 1 > 0 be uniform continuity constants for φ on W 1 , i.e.
Then for any 0 < h < h 1 : Proof. We adapt the proof of Prop. 1.18, p.30 in [53]. Let W 1 ⊂⊂ W be an open subset of W and consider another open subset W 2 such that W 1 ⊂⊂ W 2 ⊂⊂ W . For > 0, let φ ∈ C 0 (W ) be such that (see [54]): The first two terms in the right-hand side of (2.27) are controlled by owing to the previous discussion; as for the last term, we get from the Cauchy-Schwarz inequality: Multiplying both sides by α and integrating over W 1 yields: , where C is the constant supplied by the condition (iii) in the statement of Lemma 4. By assumption, α −1 is a continuous function on W 2 , which implies by Lemma We conclude this subsection with the desired density result of C 1 functions in V ω : Proof. The proof of the density (i) of C ∞ (W ) ∩ V α in V α relies on a partition of unity argument and on the properties of Lemma 5 and Proposition 1, exactly along the lines of the proof of Theorem 5.15 in [53], to which the reader is referred for details. The density (ii) of C 1 (U ) in V ω follows then from the density of C ∞ (W ) ∩ V α in V α with α = ω • η|Dη| and by composition with the C 1 diffeomorphism η.
Remark 8. Corollary 1 does not imply the density of C 1 (U ) in V ω . A result of this kind would require careful regularity assumptions on ∂U and on the behavior of ω near ∂U .

Trace theorem and Poincaré inequality in V ω
In this section, the trace operator on Γ is defined and studied for functions in the weighted space V ω , or equivalently for functions in V α on Γ × {0} = {(y, 0) | y ∈ Γ}. In the sequel, with a little abuse of notations, the latter set Γ × {0} is identified with Γ.
Proposition 2 (Trace theorem). Let ω ∈ C 0 (U, R * + ) be a positive weight on U . The trace operator Proof. Introducing α = ω • η|Dη|, using the change of variables (2.9) and the density result of Corollary 1, it is enough to prove that there exists a constant C > 0 such that: Let us consider the following partition of Γ: and is the parameter featured in the separation condition (2.2). Then Vα , which implies (2.30) and therefore terminates the proof of Proposition 2.
Remark 9. The proof of Proposition 2 supplies the existence of the trace on Γ of an arbitrary function v ∈ V α , which we shall also denote by v |Γ .
For later purposes (see Section 2.5), we shall need the surjectivity of the above trace operator; this is the purpose of the next proposition.
To this end, consider ξ and K be as in the proof of Proposition 2. For an arbitrary function v 0 ∈ L 2 (Γ), we define v by the formula v(y, s) := v 0 (y)ξ(s) a.e. in W. Obviously, v(y, 0) = v 0 (y) and ∂ s v(y, s) = v 0 (y)∂ s ξ(s), whence the Cauchy-Schwarz inequality yields: Hence v is a function in V α such that v |Γ = v 0 , which is the desired result. (ii) If (H2) is satisfied, then for an arbitrary function v 0 ∈ L 2 (Γ), we simply define v by the formula: v(y, s) := v 0 (y) a.e. in W.
Then, clearly ∂ s v = 0; what's more, one has v ∈ L 2 α (W ) as follows from the following estimate: where g α is as in the statement of (H2). Hence the function v belongs to V α , and so v := v • η −1 is an element of V ω which has the desired properties owing to Lemma 2. 13 We now prove a Poincaré-type inequality in the spaces V ω under the additional assumptions (H2) and (H3) about the weight ω.

Well-posedness of the variational problem (2.5)
We are now in a position to state and prove the main result of this section.

Proposition 5.
Let ω ∈ C 0 (U, R * + ) be a positive weight on U , satisfying the assumptions (H1) to (H3). Then: (i) For any function f ∈ L 2 ω −1 (U ), there exists a unique solution u ∈ V ω to the variational problem Proof.
(i) The assumption f ∈ L 2 ω −1 (U ) ensures that v → U f vdx is a continuous linear form on V ω owing to the Cauchy-Schwartz inequality: Moreover, the Poincaré inequality of Proposition 4 ensures the coercivity on V ω of the bilinear form Hence, the classical Lax-Milgram theorem (see e.g. [32]) yields the existence and uniqueness of a solution u ∈ V ω to (2.37). (ii) Since ω satisfies (H2), it follows from Proposition 3 (ii) that for any v 0 ∈ C 0 (Γ), there exists a function v ∈ V ω such that v |Γ = v 0 and β · ∇v = 0. Taking v as an admissible test function in (2.37) and using once again the change of variables (2.9) yields: whence (2.38).
Remark 11. Problem (2.5) or (2.37) can be solved, for an arbitrary choice of weight ω satisfying (H1) to (H3), if the right hand side satisfies f ∈ L 2 ω −1 (U ). This holds true, for example when f belongs to L ∞ (U ), as soon as the weight ω satisfies: ω −1 ∈ L 1 (U ). The latter property is not a consequence of (H1) to (H3) as explained in Remark 14.
Remark 12. If (H2) is not satisfied, which is the case if for example ω blows up "too fast" near some part of the boundary of U , functions v ∈ V ω are expected to vanish near this part of ∂U and then V ω may not contain all functions which are constant along the characteristic curves of β.

Numerical methods for integration along normal rays
We now instantiate the general setting of Section 2 to the shape optimization context (2.7) outlined in the introduction. As we have mentioned, it relies heavily on the notion of signed distance function, which we briefly recall in the next Section 3.1 for the convenience of the reader. We then clarify, for comparison purposes, some of the algorithmic stages required by the direct integration of (2.6) along rays in Section 3.2-such as the numerical computation of the principal curvatures κ i of the considered shapes, and the delicate detection of their skeleton on unstructured meshes-. We discuss next in Section 3.3 the construction of suitable weights ω that allow to to solve accurately the variational problem (2.5) by means of P 1 conforming finite elements. Finally, in Section 3.4, we compare on several 2-d analytical examples the numerical calculations of (2.6) by means of our variational formulation (2.5) with those produced by direct integration of this formula along rays.

A short reminder about the signed distance function
In this section, we consider a fixed, bounded and Lipschitz 'hold-all' open domain D ⊂ R n , as well as a bounded Lipschitz open subdomain Ω ⊂ D. We recall the main definitions and properties of the signed distance function that will be relevant for the shape optimization applications of Section 4, referring to [13,25,26] for further details. where || · || denotes the Euclidean norm of R n . The signed distance function d Ω : D → R to the domain Ω is defined by: We then recall the definitions of the projection mapping onto ∂Ω and of the skeleton of Ω, which are illustrated on Figure 1a.
(i) The (non empty) set of projections Π ∂Ω (x) of a point x ∈ D onto ∂Ω is the set of points y ∈ ∂Ω achieving the minimum in (3.1), that is: (ii) The skeleton Σ ⊂ D of Ω (relative to D) is the set of points x ∈ D such that Π ∂Ω (x) is not a singleton.
(iii) For x ∈ D\Σ, the unique element in Π ∂Ω (x) is denoted by p ∂Ω (x) and is called the projection of x onto ∂Ω.
In the following reminders on the signed distance function d Ω , we suppose for simplicity that Ω is a domain of class C 2 (although this assumption may not be minimal).
Proposition 6 (Differentiability of d Ω ). The signed distance function d Ω is differentiable at any point x ∈ D\Σ, and it is not differentiable on Σ. The gradient ∇d Ω is an extension of the unit normal vector n to ∂Ω pointing outward Ω, in the sense that it satisfies the properties: ∀y ∈ ∂Ω, ∇d Ω (y) = n(y) and ∀x ∈ D\Σ, ||∇d Ω || = 1. where ζ − (y) and ζ + (y) are the maximum distances at which ray(y) hits either the skeleton Σ or the boundary ∂D of the hold-all domain: The functions ζ − and ζ + are continuous on ∂Ω (see [15] or [45]).
Let us now turn to properties related to the second-order regularity of the signed distance function d Ω , still under the assumption that the domain Ω is at least of class C 2 . In all what follows, we denote for any y ∈ ∂Ω by n(y) the outward normal to Ω. The tangential gradient ∂ y n(y) (with respect to the coordinates of ∂Ω) is a symmetric tensor whose eigendecomposition involves the principal curvatures (κ i (y)) 1≤i≤n−1 and the associated principal directions (τ i (y)) 1≤i≤n−1 of ∂Ω at y. Note that for any unit extension n : V → R n of n to some tubular neighborhood V ⊂ D of ∂Ω satisfying n(y) = n(y) for any y ∈ ∂Ω and || n(x)|| = 1 for any x ∈ V, it holds that ∇ n(y) = ∂ y n(y). The next proposition is related to the differentiability properties of the projection p ∂Ω , or equivalently to the second-order differentiability of d Ω (see [8,13,15,37,48]).
Proposition 7 (Differentiability of p ∂Ω ). The following properties hold true: (1) For any point x ∈ D\Σ, one has: 1 + κ i (p ∂Ω (x))d Ω (x) > 0, i = 1, ..., n − 1. 16 (2) The projection p ∂Ω is differentiable on D\Σ and Remark 13. The projection p ∂Ω is not differentiable on the boundary ∂Σ of the skeleton Σ. Indeed (see Figure 1a and [15,37]), a point C ∈ D belongs to the closure Σ either because it lies on Σ, or because it is a center of curvature of ∂Ω, i.e. there exists y ∈ Π ∂Ω (C) and i = 1, ..., n − 1 such that Since d Ω is Lipschitz, it follows from Rademacher's theorem (see [33], §3.1.2) that the skeleton Σ has zero Lebesgue measure. However, if Ω is only supposed to be of class C 2 , the closure Σ may have non zero Lebesgue measure (see [48] for an example). For Σ to be also of null Lebesgue measure, a little more regularity is required about Ω, namely, that it be at least of class C 3 [15,48].
From (3.7) and (3.8), we may calculate the Hessian matrix of d Ω and the Jacobian |Dη|. The last proposition of this section sums up the characteristics of the shape optimization setting in the perspective of the previous Section 2.

Proposition 8 (Shape optimization setting).
Let Ω ⊂ D be a domain of class C 2 ; then the signed distance function d Ω is of class C 2 on the open set U := D\Σ. Hence β := ∇d Ω is a vector field of class C 1 on U ; the associated flow map η : W → D\Σ is a diffeomorphism of class C 1 , whose expression reads: ∀(y, s) ∈ W, η(y, s) = y + s∇d Ω (y), (3.9) where W is the set defined by (2.3). The inverse flow mapping η −1 : U → W is given by The divergence of the vector field β = ∇d Ω and the Jacobian |Dη| of the flow map η are respectively given by (1 + sκ i (y)). (3.12)

Computing curvatures and detecting the skeleton for direct integration along the rays
Before going to the numerical aspects of our variational method, we clarify important practical details which are required in the implementation of the direct integration along characteristics involved in the calculation of (2.6). In Section 3.2.1, we discuss the delicate issue of detecting the triangles in the computational mesh of D crossing the skeleton Σ of Ω when traveling along the rays (note that this step is not required by our variational method). Then, we detail in Section 3.2.2 the method we used to compute the curvature κ (there is only one curvature κ := κ 1 in 2-d) required in the line integral formula (1.3). Note that these steps serve only for comparison purposes with our variational method. These are fairly classical numerical issues which could otherwise be addressed with more sophisticated techniques, see e.g. [55,27]. For our present numerical applications, the hold-all domain D is equipped with a simplicial mesh T featuring a discretization of the domain Ω as a submesh. However, the only information we use about Ω is an accurate approximation d h of the signed distance function d Ω as an element of the space V h of Lagrange P 1 finite element, where h is the maximum mesh element size. Multiple methods are available to achieve the latter calculation [57,63]; here we rely on the numerical algorithm from the previous work [23].

Detection of the skeleton Σ of Ω and identification of normal rays
The numerical detection of Σ in the course of the identification of the set ray(y) for some given point y ∈ ∂Ω is achieved by assessing the following criterion (independently of the dimension), holding at the continuous level (see [26]): ∀y ∈ ∂Ω, ∀z 0 , z 1 ∈ ray(y), d Ω (z 1 ) − d Ω (z 0 ) (z 1 − z 0 ) · ∇d Ω (y) = 1. In our implementation, when computing the ray emerging at some point y ∈ ∂Ω (which is detected by the fact that d h (y) = 0), we travel the triangles in the mesh T in the normal direction n = ∇d h (y), and we stop the calculation of the ray in the triangle T ∈ T where the entering and exiting points z 0 and z 1 satisfy: where tolRay a small tolerance (set to 0.3 in our implementation). This provides meanwhile an approximate location of the skeleton Σ, up to a tolerance of the order of the mesh size. Our criterion (3.13) differs from that used in the related works [6,50]. In there, the authors detect Σ by looking at changes in the monotonicity of the signed distance function d h along the ray, i.e. they rely on the following property of the (continuous) signed distance function d Ω : (3.14) Our personal experiment with the above criterion suggests that it may sometimes fail to detect the skeleton Σ accurately, because such a change in monotonicity may simply not occur when the ray is supposed to cross Σ in the neighborhood of center of curvatures (see Remark 13). Our criterion (3.13) may also fail depending on the chosen tolerance parameter, but it offered visible improvements ( Figure 3) in our academic test-cases. Note that when integrating along the ray, the last triangle, where the skeleton is hit, is included in the integration.

Estimating the curvature κ of a 2-d subdomain based on its signed distance function
In this part, we detail our method for the numerical approximation, in 2-d, of the principal curvature κ (there is only one in 2-d) from the knowledge of a P 1 discretization d h of the signed distance function d Ω at the nodes of the mesh T . We essentially rely on the fact that in 2-d, κ is given by the trace of ∆d Ω on the boundary ∂Ω (in view of (3.11). In 3-d, the estimation of ∆d Ω would not be sufficient to evaluate the values of both principal curvatures κ 1 and κ 2 : these could e.g. be computed from the eigenvalues of the Hessian matrix ∇ 2 d Ω .
Our first step towards estimating ∆d Ω consists in calculating a P 1 interpolation g h ∈ V h × V h of the piecewise constant gradient ∇d h by solving the following variational problem: The approximation of the divergence div(∇d Ω ) is then calculated as the (piecewise constant) divergence of the reconstructed field g h . Unfortunately, this procedure generally produces a very noisy approximation characterized by a lot of spurious oscillations when the mesh resolution increases (see Figure 4). In order to overcome this difficulty, we calculate a regularization κ h of this noisy estimation with a Laplace kernel, namely we solve: where γ h > 0 is a regularization length scale (equal to 3h max where h max is the maximum edge length in the mesh). This procedure yields satisfying results in practice (even with shapes Ω characterized by  discontinuous curvatures, up to some over smoothing near the discontinuities), although we do not have a proof of convergence of the approximation κ h towards the exact function ∆d Ω .
Let us illustrate our method by considering the example of an ellipse Ω inside a square-shaped hold-all domain D: let D and Ω be defined by where a = 1.5 and b = 1. The skeleton of Ω is explicitly known in this case: Σ = {(x 1 , 0) | |x 1 | < a − b 2 /a} and the curvature κ of ∂Ω at a point y = (y 1 , y 2 ) is given by κ(y) = ab/γ 3 with γ = b 2 a 2 y 2 1 + a 2 b 2 y 2 2 . The difference between the exact curvature κ(y) of ∂Ω and its reconstruction κ h (at the boundary nodes discretizing ∂Ω) using both procedures (3.15) and (3.16) is represented in Figure 4

Admissible numerical weights built upon the signed distance function
In this section, we discuss the numerical resolution of the variational problem (2.5) in the shape optimization context (2.7) and (2.8) (see also Proposition 8). The numerical setting is the same as that of the previous Section 3. 2: the hold-all domain D is equipped with a simplicial mesh T (i.e. composed of triangles in 2-d, or tetrahedra in 3-d, although our method would work with other kinds of meshes) featuring a discretization of the domain Ω as a submesh, and we rely on the Lagrange P 1 finite element method for the discretization of (2.5) on account of its robustness and ease of implementation. In other terms, the approximation u h to the solution u ∈ V ω of (2.5) is sought in the space V h of continuous piecewise linear functions on each simplex of T .
We start in Section 3.3.1 by introducing some motivations for the use in (2.5) of weights ω vanishing on the skeleton, and we provide (in Lemma 6 below) a formula for analytical and admissible weights satisfying approximately this property. We then discuss the issue of computing numerically these weights in Section 3.3.2. Finally, we perform a few numerical experiments in Section 3.3.3 where we show that weights vanishing on the skeleton make P 1 finite elements able to capture discontinuous test functions across the skeleton, which allows to confirm numerically the latter motivations.

Motivations for weights vanishing on the skeleton
As already mentioned, our final target is to calculate the trace (2.6) on Γ = ∂Ω of u. In the continuous setting of Section 2, this trace does not depend on the choice of the weight ω as long as it fulfills (H1) to (H3), so that in principle, any such weight could be used. However, when U is a cracked domain (typically, the skeleton Σ is a crack in the working domain, i.e., U = D\Σ) and the crack is not explicitly discretized in the mesh T , then the most simple choice ω = 1 (which is an admissible weight on account of Lemma 6 below) might not work well in practice. Indeed and as we shall illustrate below in Section 3.3.3, test functions of (2.5) which belong to V ω with ω = 1 are in general discontinuous across Σ and not well captured by P 1 finite elements. In many shape optimization applications, discretizing the skeleton Σ of the current domain Ω at every iteration of the optimization process or resorting to discontinuous finite elements (in order to be able to set ω = 1) is very inconvenient, for instance if working with fixed structured meshes as it is performed in many applications built on the level set method [4,61].
Therefore, we shall be interested in determining weights ω adapted to our commitment to use the space V h of Lagrange P 1 finite elements for the resolution of (2.5) without the need for an accurate discretization of the skeleton Σ (alternative approaches could be to use discontinuous finite elements close to the skeleton, or to have a zero weight on the degrees of freedom corresponding to modes close to the skeleton and to remove the null space in the corresponding linear system, but they seemed to be more complicated to implement, at least to us). The weight ω should be chosen in such a way that arbitrary functions v ∈ V ω are well approximated (in the norm || · || Vω ) by functions v h ∈ V h , as is reflected by the classical Céa's lemma (see e.g. [32]): for a constant C > 0 (which possibly depends on ω). The space V h of Lagrange P 1 elements is a conformal finite element space in the sense that the inclusion V h ⊂ V ω always holds (because functions of V h are smoother than those of V ω ), however V h may be "too small" to guarantee a correct approximation of discontinuous solutions u ∈ V ω in the sense of (3.18). Heuristically, and without looking for a very precise statement, these considerations call for the choice of a weight ω almost vanishing on Σ, so that the approximation error ||u − v h || Vω in (3.18) attributes a lesser weight to a neighborhood of Σ where u is expected to be discontinuous while the functions v h ∈ V h are continuous.
We now provide explicit candidates for weights ω which fulfill the conditions (H1) to (H3) while taking small values near the skeleton, that we are going to use in our practical implementations. Lemma 6. For any real numbers q ≥ 0, r ≥ 0, the weight satisfies the conditions (H1) to (H3) (this includes in particular the constant weight ω = 1/2 for q = r = 0).
• When it comes to solving (2.5) by relying on Proposition 5 with some data f ∈ L ∞ (U ), it is useful to observe that f belongs to L 2 ω −1 (U ) as soon as the weight ω satisfies ω −1 ∈ L 1 (U ) (see Remark 11), which is the case if it is of the form (3.19) with r < 2. In the following numerical experiments, we shall see however that using values for r which are larger than 2 still provides good results in practice: in general, taking higher values of q and r yields a faster decay of ω near the skeleton.
• Taking r > 0 ensures that the weight (3.19) will vanish at points x ∈ Σ that are centers of curvatures (for which ∆d Ω blows up). Taking q > 0 allows to make sure that ω = 1/2 on ∂Ω whatever the value of r, and to accentuate the decay of ω near the skeleton. Not that in general, ω will not vanish on points x ∈ Σ that are not centers of curvatures. However, it still takes very small values on Σ and it is convenient to use in the implementation. There is no unicity of weights appropriate for the numerical computation and variants can easily be imagined. • For the most general setting where U is an arbitrary open set, the methodology of this section extends naturally by considering weights ω which vanish on the cracked parts of ∂U that are not explicitly meshed.

Numerical computations of the weight based on the Laplacian of the signed distance function
For our numerical applications below, we shall use the weights of Lemma 6 in the definition of our variational formulation (2.5). This requires the computation of the Laplacian ∆d Ω on the triangulated mesh T based on the P 1 estimation of the signed distance function d Ω . For this purpose, we use the same regularization method outlined in Section 3.2.2 for the computation of the numerical curvature. Importantly, from Proposition 5, the variational formulation (2.5) is rather insensitive to the choice of the weight ω, and as a result the estimation of ∆d Ω does not need to be very accurate as long as it takes large values near the skeleton (as we shall illustrate below in Section 3.3.3). In contrast, the estimation of the principal curvatures κ i for the direct method would need to be accurate. From a numerical standpoint, the formula (3.19) featuring ∆d Ω at the denominator is convenient to obtain numerically vanishing weights near the skeleton (even if (3.19) truely vanish for r > 0 at centers of curvatures). Indeed, ∇d Ω is discontinuous across the skeleton, which should reflect in high numerical values of ∆d Ω when computing numerically the divergence div(∇d Ω ) with the method of Section 3.2.2.

3.3.3.
Assessing the choice of the weight ω when using a P 1 discretization: generating numerical test functions constant along rays In this section, we perform several numerical experiments about the influence of the choice of the weight function ω in the resolution of the variational problem (2.5) using the Lagrange P 1 finite element method. With this perspective in mind, and in the 2-d numerical setting described in Sections 3.2 and 3.3, we consider the issue of generating numerical functions v ∈ V ω which are constant along the normal rays to the shape Ω. Namely, we solve the boundary-value problem for given data v 0 ∈ L 2 (∂Ω).
Using the variational setting of Section 2, we show that it is possible to obtain the solution v to (3.24) by solving a variational problem of the same nature of (2.5).

Proposition 9.
Let ω ∈ C 0 (D\Σ, R * + ) be a weight satisfying (H1) to (H3). There exists a unique solution v ∈ V ω to the following variational problem: The solution v is independent of the choice of ω as long as (H1) to (H3) are satisfied, and it is given by the formula: Proof. To see that (3.24) and (3.25) are equivalent, it is sufficient to take v ∈ V ω constant along rays in (3.25) as in the proof of Proposition 5, which yields v = v 0 on ∂Ω, and then ∇d Ω · ∇v = 0.
The formulation (3.25) is to be compared to the so-called Galerkin Least Square formulation and SUPG methods for advection-reaction problems [32], with the difference that usual assumptions of uniformly bounded divergence (which do not hold in our applications) are replaced with the regularity assumptions of Section 2.1.
Remark 15. The problem of building constant functions along normal rays of the form (3.24) may be solved on unstructured meshes owing to variants of the fast marching algorithm; see e.g. [19].
We now verify that a good approximation of the solution v to (3.24) (or more precisely (3.25)), which is in particular discontinuous across the skeleton Σ of Ω, can be obtained either by truncating manually the skeleton from the computational mesh T and taking the weight ω = 1, or by selecting a weight ω taking "small" values near Σ. For the latter experiment, we use the weight ω = 1/(1 + |d Ω | 3.5 |∆d Ω | 2 ) (see also Remark 14 about the "small" values of ω near Σ). Of course, removing the skeleton from the mesh is not a straightforward task in full generality but it is performed here for the sake of comparison.
Let us consider again the ellipse example of (3.17). We consider two different computational meshes T and T . The former is a triangular mesh of D, and the latter T is a triangular mesh of D\Σ (i.e. the skeleton Σ has been manually removed). In both T and T , the considered shape Ω is explicitly discretized as a submesh; see Figure 5 for an illustration. The variational problem (3.25) is numerically solved for a boundary datum v 0 ∈ L 2 (∂Ω) given by (see Figure 6d): ∀(y 1 , y 2 ) ∈ ∂Ω, v 0 (y 1 , y 2 ) = cos(3y 1 ) 2 + 20y 2 , and the computed finite element solution is plotted on Figure 6 in the following three situations.
• The mesh T is used, in which D and Ω ⊂ D are meshed explicitly, and where a thin layer around Σ has been manually removed (see Figure 5c). The solution v to (3.25) is computed using the constant weight ω = 1 and the result is displayed on Figure 6a. As expected, the fact that Σ is absent from T allows the numerical solution v to have very different values on either sides of Σ. • The mesh T of D (where Σ has not been removed) is used (see Figure 5b), and v is computed with the constant weight ω = 1; the result is represented on Figure 6b. The formulation (3.25) proves numerically stable with the choice ω = 1, but it tends to over smoothen the sharp discontinuities of v near the skeleton Σ, which results in a loss of accuracy for the extension problem (3.24). • The mesh T of D is used again, but the solution v to (3.25) is now computed by using the weight ω = 2/(1 + |d Ω | 3.5 |∆d Ω | 2 ); the result is represented on Figure 6c. size), where sharp variations are observed, as expected. The numerical procedure in this case seems therefore to achieve the same order of accuracy than in the experiment using the truncated mesh T .

Numerical comparisons between the variational method and direct integration along rays
We now investigate the numerical evaluation of the function u ∈ L 2 (∂Ω) in (2.6) for several 2-d academic configurations as far as D, Ω and the function f are concerned. In particular, we compare the evaluation of u obtained by direct integration along rays (i.e. implementing directly the formula (2.6)) to that obtained by solving the variational formulation (2.5) on meshes T (resp. T ) of D in which Ω is explicitly discretized and the skeleton Σ of Ω is not removed (resp. is removed).

A domain with trivial skeleton: the case of a circle
We first consider the case where Ω is a disk enclosed in a larger disk D: Figure 7. In this case, the skeleton Σ is reduced to the point 0. The considered function f is: which belongs to the finite element space V h , and is therefore amenable to an exact integration when the travel procedure along rays of Section 3.2.1 is used. In this situation, the sought function u, given by (2.6), is known analytically; a calculation in polar coordinates indeed yields: ∀(y 1 , y 2 ) ∈ ∂Ω, u(y 1 , y 2 ) = 8 3 y 2 . (3.28) Comparisons are displayed on Figure 8 between the version of u obtained after direct integration along rays, and the numerical solutions of (2.5) obtained for various choices of weight functions ω on the mesh T of D (where Σ has not been removed). We observe a good match between the variational and the direct method. As expected, the solutions computed thanks to our variational method are less accurate when the constant weight ω = 1 is chosen. A significant increase in accuracy is achieved by selecting a weight ω vanishing at the center of the circle, i.e. ω = 1/(1 + |d Ω | 3.5 |∆d Ω | 3.5 ). Note that this weight does not fulfill the condition ω ∈ L 1 ω −1 (U ) (see Remark 12), but works well in numerical practice.

A C 2 domain with non trivial skeleton
We take on the example (3.17) where Ω is an ellipse inside a square-shaped hold-all domain D, reusing both meshes T and T depicted on Figure 5. To evaluate the influence of the skeleton on our numerical method, we compute the quantity (2.5) for a function f which is supported in Ω only (see Figure 9): The presence of the constant 1 in (3.29) avoids the "simplification" of having f vanishing on Σ. After an explicit calculation based on (2.6), the exact solution is given by: The numerical trace u (which is extended along rays using (3.25) to ease the visualization) obtained by solving the variational formulation (2.5) on T or T , and for various choices of weights, is represented on Figure 10. On Figure 11, we compare for each of these strategies the boundary values of the numerical trace u to that computed with (1.3) by direct integration along rays, and to the analytical expression (3.30).
We note that the solutions obtained by integrating along rays are generally characterized by small spurious oscillations. We explain this error by the fact that our criterion (3.13) detects the skeleton up to an error proportional to the mesh size (the actual location of the skeleton is not estimated within the last travelled triangle; see Figure 12). Such a procedure could be improved by using more accurate skeleton detection methods; see for instance the works [9,10] in the field of computer graphics. We verify that the amplitude of these oscillations decreases with the size of the mesh, as shown on Figure 11d.  The triangle paths are colored in transparent blue (the darker, the more often visited); depicting that some mesh triangles might not be crossed by all the rays or some triangles may be visited more often than others.

A C 1 domain, with discontinuous curvature
We now consider the case where Ω is a stadium, i.e. the reunion of a rectangle and two half-disks. Define (see Figure 13) and 32) The domain Ω is not of class C 2 , and the curvature κ of the boundary ∂Ω is discontinuous at the points where x 1 = ±0.5. Hence, this example does not fall into the admissible setting of Proposition 8, and there is, in principle, no guarantee that our variational method based on (2.5) should still work.
In this example, the skeleton Σ of Ω is explicitly given by Σ = {(x 1 , 0) | |x 1 | ≤ 0.5}. We calculate the function u in (2.6) associated to the datum function f defined by (3.29); the analytical solution u for this problem is given by In particular, u is ill-defined at the points x ∈ ∂Ω where the curvature is discontinuous.  Both methods for the numerical evaluation of u (i.e. the direct integration along rays and our variational method) are applied, and the results are represented on Figure 14. Even though this example does not fit into the admissible setting of Section 2, the results indicate that our numerical method still works, up to over-smoothing inaccuracies of the computed solution u near the discontinuities of the curvature κ.

A test case involving a Lipschitz domain Ω with angular corners
In this last example, we consider the following situation: together with the datum function f : as illustrated on Figure 15. Again, the theoretical framework of Section 2 does not apply to the present situation because: • The vector field ∇d Ω is not of class C 1 on U and is not U -filling: the reunion of all the rays emerging from points y ∈ ∂Ω does not cover the whole set U = D\Σ; • The flow η of this vector field is not of class C 1 , because the curvature κ of the boundary (and even the normal vector) is discontinuous at the corners. Numerically, |∆d Ω | blows up in the whole area filled by the skeleton and the set of points that are not covered by normal rays. Still, the quantity u given by (2.6) can be calculated analytically wherever it makes sense. A few elementary, albeit technical computations yield: 1+t 2 t and φ is a primitive function of λ → cos(6t(1 + λ/ √ 1 + t 2 )) 2 (1 + λ/(1 + t 2 ) 3/2 ). A numerical approximation to this function u is computed by using either a direct integration along rays or our variational method based on (2.5), on a single mesh T where the skeleton Σ is not removed, and using three possible choices of the weight ω. The results are represented on Figure 16. In all cases, the observed numerical inaccuracies, characterized by very high values, are concentrated near the angular corners of Ω, while the method yields satisfying accuracy on the remaining smooth parts of ∂Ω. Again, we observe the poor accuracy associated to the choice of a constant weight ω = 1. For this example, a uniformly small weight ω = 1e −10 seems to yield a better accuracy than the previous choice ω = 2/(1 + |d Ω | 3.5 |∆d Ω | 3.5 ).

Applications to maximum and minimum thickness constraints in shape optimization
In the shape optimization setting which is the original motivation for our work, we show how the general method proposed in Section 2 and the inferred numerical methods in Section 3 allow to efficiently implement geometric constraints, namely maximum and minimum thickness constraints. Here no comparison is made between our new variational approach and the previous method (using direct integration along rays) to evaluate the shape derivatives of the related shape functionals. The optimized shapes and topologies resulting from the variational method are very similar to those obtained in the previous works [6,50]. There is no clear gain in computational time but there is a very substantial simplification of the implementation (which would be tremendous in 3-d).

Shape optimization of linearly elastic structures
We consider shapes, that is, smooth bounded domains Ω enclosed in a fixed, bounded and Lipschitz hold-all domain D ⊂ R n . Any such shape Ω is clamped on a part Γ D of its boundary, and traction loads g ∈ L 2 (Γ N ) are applied on a disjoint region of ∂Ω; the complement Γ := ∂Ω\(Γ D ∪ Γ N ) is traction-free and body forces are omitted for simplicity; Γ is the only region of the boundary ∂Ω which is subject to optimization. In this situation, the displacement u Ω of the shape is characterized as the unique solution in H 1 (Ω, R n ) to the linearized elasticity system: where e(u) := 1 2 (∇u + ∇u T ) is the strain tensor associated to the displacement u and A is the Hooke's law, defined for any symmetric n × n matrix by Ae(u) = 2µe(u) + λTr(e(u))I, involving the Lamé coefficients λ, µ which characterize the physical properties of the constituent material.
In this context, we consider structural optimization problems of the form min Ω⊂D J(Ω), s.t. P (Ω) ≤ 0, (4.2) where J(Ω) is a performance criterion, which will typically be the volume Vol(Ω) or the compliance C(Ω) of shapes: Vol(Ω) = involving the signed distance function d Ω to Ω and a given smooth function j : R → R.
In order to achieve the numerical resolution of (4.2), we rely on Hadamard's method of boundary variations (see for instance [1,4,40,51,58]), whereby variations of a given shape Ω are considered under the form Ω θ := (Id + θ)(Ω), θ ∈ W 1,∞ (D, R n ), ||θ|| W 1,∞ (D,R n ) < 1. (4.5) Definition 8. A function F (Ω) of the domain is shape differentiable at a shape Ω if the underlying mapping θ → F (Ω θ ), from W 1,∞ (D, R n ) into R, is Fréchet differentiable at θ = 0. The corresponding derivative θ → F (Ω)(θ) is the shape derivative of F at Ω and the following expansion holds in the vicinity of θ = 0: Remark 16. In our setting, the set of admissible deformation fields θ used in (4.5) is restricted to smooth vector fields, which vanish on the non optimizable subset Γ N and Γ D of the boundaries of shapes.
The shapes derivatives of the volume and compliance functions defined in (4.3) are classically given by (see for instance [7]): where we recall that dy stands for the surface measure on Γ. As far as the shape derivative of the geometric constraint P (Ω) in (4.4) is concerned, the following (equivalent) surface and volumetric expressions were obtained in our previous work [2]: where the "Eulerian" derivative d Ω (θ) of the signed distance function d Ω is defined on U = D\Σ by: ∀x ∈ D\Σ, d Ω (θ)(x) = −θ(p ∂Ω (x)) · n(p ∂Ω (x)), (4.8) and the scalar function u : Γ → R reads: 1≤i≤n−1 (1 + κ i (y)d Ω (s))ds, as follows from an application of the coarea formula with the help of the material recalled in Section 3.1.
Using the conclusions of Section 2, the function u in (4.7) may be conveniently evaluated by solving the variational problem: for a suitable weight ω satisfying (H1) to (H3), as discussed in Section 3. Note how easily (4.7) is retrieved by taking the test function v = −d Ω (θ) in (4.9). Since θ · n ∈ L ∞ (∂Ω), the derivative given by (4.8) is indeed an admissible test function d Ω (θ) ∈ V ω , for it belongs to L ∞ (D \ Σ) and is constant along normal rays.
In our numerical implementation, we set ω = 1/(1 + 100|d Ω ∆d Ω | 3.5 ) in order to solve (4.9) with P 1 finite elements, the constant 100 being selected to increase the slope of the weight near the skeleton.
Let us now outline briefly the numerical implementation of the above shape optimization framework. All the involved finite element computations are performed within the FreeFem++ environment [39]-whether they are related to the resolution of the physical system (4.1) or to that of our variational problem (4.9). The resolution of the constrained optimization program (4.2) is carried out by the discretization of a specific gradient flow in the spirit of [56], see the details of our algorithm in our recent work [36]. When it comes to representing shapes and their evolution in the course of the iterative resolution of (4.2), the level set based mesh evolution method of [3,35] is used, as a convenient combination of the 'classical' level set method for shape and topology optimization [7,61] and the mmg open-source mesh modification algorithm [20].

Shape optimization under a maximum thickness constraint
We first consider structural optimization problems featuring a maximum thickness constraint. Following the work in [6,50], a shape Ω ⊂ D is said to have maximum thickness lower than d max > 0 provided: (4.10) Loosely speaking, this amounts to saying that the skeleton Σ of Ω lies at a distance of at most d max from the boundary ∂Ω. Following closely [6,50], the pointwise constraint (4.10) is relaxed into a single integral constraint formulated in terms of the following penalty functional P MaxT (Ω): and h is a regularized Heaviside function centered at d max /2: The parameter α f in (4.12) tunes the level of regularization; in our context, it is set to α f = 4hmax 5dmax , where h max is the maximum edge length of the computational mesh of D.
A simple calculation yields the shape derivative of P MaxT (Ω): (4.13) where the term involving d Ω (θ) is then computed using the variational formulation (4.9).
In the forthcoming examples of Sections 4.2.1 and 4.2.2, we solve the optimization problem of minimizing the volume Vol(Ω) of the structure while imposing that the compliance C(Ω) do not exceed a given threshold g max , as well as the maximum thickness constraint (4.11), namely: (4.14)

Optimization of the shape of a two-dimensional arch
Our first test-case reproduces that considered in §8.1.1 of [6], whose setting is displayed on Figure 17a: in a square-shaped domain D, the considered shapes are clamped on both their bottom left and bottom right corners, and a vertical load g = (0, −20) is applied at the middle of their bottom side. Starting from an initial shape arbitrarily perforated with several holes, the optimization problem (4.14) is solved with and without including the maximum thickness constraint, using the numerical values g max =7.00 and d max =0.12 when they are relevant.
The resulting optimized shapes in both cases are displayed on Figs. 17b and 17c. Several intermediate shapes as well as the convergence histories of the computation are displayed on Figs. 18 to 20. The obtained shapes are quite analogous to those obtained by [6,50] where the calculation of the shape derivative of P MaxT (Ω) relied on a direct numerical integration along rays.

Optimization of a two-dimensional MBB-Beam
We consider now the classical MBB-Beam test-case depicted on Figure 21a: in a box D with dimensions 3 × 1, a material shape Ω is constrained to no horizontal motion on the left boundary, and to no vertical motion on the bottom right corner. A vertical load g = (0, −10) is applied on the top left corner, and the optimization problem (4.14) is considered again, with the numerical values for the thresholds g max =30.00 and d max =0. 16. The resulting optimized shapes with and without including the maximum thickness constraint in (4.14) are represented on Figs. 21b and 21c and the convergence histories of the computation are shown on Figs. 22 to 23.

Shape optimization examples under a minimum thickness constraint
We now turn to the implementation of a minimum thickness constraint thanks to our variational method. Following [6,50], we say that a shape Ω has minimum thickness greater than d min if ∀y ∈ ∂Ω, ζ − (y) < −d min /2. In other words, the boundary ∂Ω is at a minimum distance d min /2 of the part Σ ∩ Ω of the skeleton located inside the shape. Enforcing a minimum thickness as a hard constraint (i.e. rather than a penalty term in the objective function) is by no means a straightforward task, because: (1) Our definition (4.15) of minimum thickness involve the distance to the skeleton ζ − , which is not differentiable with respect to the shape, (2) It is not clear how to formulate (4.15) by mean of a penalty functional such as (4.11) to penalize localizations on the shape that do not meet the thickness requirement, (3) Even if we were able to enforce the constraint at each iteration, such would prevent topology changes to occur naturally, which would be an issue in numerical practice.
We propose in the following a more flexible setting to enforce such a requirement in a structural optimization problem. Elaborating on ideas proposed in [50,16,17,47], the minimum thickness requirement is implemented in the objective function rather than in the constraints: we minimize a penalty functional P MinT (Ω) for the minimum thickness under constraints on the volume and compliance of shapes, i.e. we   This strategy is expected to work because our optimization algorithm is designed to satisfy violated constraints first, before then attempting to reduce the objective function while maintaining the constraints respected.
The penalty functional P MinT (Ω) we considered in our case when solving (4.16) is taken from [50] and is specially designed to have a zero derivative when the constraint (4.15) is satisfied: (4.17) The shape derivative of P MinT (Ω) is given by P MinT (Ω)(θ) = − Ω 2(d Ω max(d Ω + d min /2, 0) 2 + d 2 Ω max(d Ω + d min /2, 0))d Ω (θ)dx. (4.18) Note that an increase in perimeter entails a decrease in the value of P MinT (Ω), but this behavior is tempered in our case by the volume constraint. Variants can be considered to address such an issue, and we refer to [50] for the details.

Optimization of the shape of a 2-d cantilever beam
We first consider the classical two-dimensional cantilever benchmark example, as depicted on Figure 24: in a box D with size 2 × 1, shapes Ω are clamped on their left-hand side, and a vertical load g = (0, −10) is applied on the middle of their right-hand side. The optimization problem (4.16) is solved with the parameter values g max =70.00 and V max =0.80. The resulting optimized shapes are displayed on Figure 25 without including minimum thickness constraint (we minimize the volume Vol(Ω) subject to the compliance constraint as in (4.14)), and with the minimum        Figure 32. One observes that for the last case of Figure 29c with d min = 0.2, the minimum thickness constraint is not satisfied everywhere but a substantial improvement is visible over the first design of Figure 29a. Notably, this approach is sufficiently flexible to guide the optimization path towards shapes with very different topologies.
(a) Optimized shape without including minimum thickness constraint.