S Amstutz

C Dapogny

A Ferrer

A CONSISTENT RELAXATION OF OPTIMAL DESIGN PROBLEMS FOR COUPLING SHAPE AND TOPOLOGICAL DERIVATIVES

In this article, we introduce and analyze a general procedure for approximating a 'black and white' shape and topology optimization problem with a density optimization problem, allowing for the presence of 'grayscale' regions. Our construction relies on a regularizing operator for smearing the characteristic functions involved in the exact optimization problem, and on an interpolation scheme, which endows the intermediate density regions with fictitious material properties. Under mild hypotheses on the smoothing operator and on the interpolation scheme, we prove that the features of the approximate density optimization problem (material properties, objective function, etc.) converge to their exact counterparts as the smoothing parameter vanishes. In particular, the gradient of the approximate objective functional with respect to the density function converges to either the shape or the topological derivative of the exact objective. These results shed new light on the connections between these two different notions of sensitivities for functions of the domain, and they give rise to different numerical algorithms which are illustrated by several experiments. Contents 7.2. The case of the linear elasticity system 26 8. Description of the numerical methods 27 8.1. Practical implementation of the regularizing operator L ε 28 8.2. Two numerical optimization algorithms 29 9. Numerical illustrations 30 9.1. Numerical assessment of the convergence of the approximate shape optimization problem to its exact counterpart 31 9.2. A self-adjoint objective function for the elasticity problem 32 9.3. A non self-adjoint objective function 32 Appendix A. Sketch of the proof of Theorem 4.5 34

Most optimal design problems arising from the fields of physics or continuum mechanics inherently feature a 'black-and-white' design, that is, a 'solid' shape Ω, embedded in a fixed computational domain D, clearly separated from the 'void' (or a different phase) D \ Ω. The geometry or topology of Ω are optimized, for instance using steepest descent algorithms based on the notions of shape [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF][START_REF] Sokolowski | Introduction to Shape Optimization[END_REF] or topological derivative (see [START_REF] Novotny | Topological derivatives in shape optimization[END_REF] for an overview). A popular paradigm -termed relaxation -to circumvent some of the difficulties posed by these shape and topology optimization problems, such as the lack of existence of an optimal design or the tediousness of numerical implementation, cuts off this physical context; see for instance [START_REF] Bucur | Variational Methods in Shape Optimization Problems[END_REF] Chap. 3 or [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF] Chap.7 for general presentations. It consists in enriching the set of admissible designs to allow for the presence of intermediate 'grayscale' regions, composed of a mixture of solid and void. The original, mathematically rigorous way to carry out this change in points of view is the homogenization method, pioneered by F. Murat and L. Tartar (see [START_REF] Murat | Calculus of Variations and Homogenization[END_REF]), which has then raised a significant enthusiasm in the shape and topology optimization community: we refer for instance to [START_REF] Allaire | Shape optimization by the homogenization method[END_REF][START_REF] Bendsøe | Topology Optimization. Theory, Methods, and Applications[END_REF][START_REF] Suzuki | A homogenization method for shape and topology optimization[END_REF], and to [START_REF] Allaire | Shape optimization by the homogenization method[END_REF] for a comprehensive presentation; the conventional representation of shapes is traded for the data of a density function h : D → [0, 1], accounting for the local proportions of material (identified by h = 1) and void (corresponding to h = 0), and of a homogenized tensor, describing the local, microscopic arrangement of solid and void, and thereby the material properties of the design. Another way to achieve this relaxation is the more heuristic and quite popular Solid Isotropic Material Penalization (SIMP) method, introduced in [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF]; see [START_REF] Bendsøe | Topology Optimization. Theory, Methods, and Applications[END_REF] for numerous applications. In a nutschell, the only retained information in the SIMP setting is that of the density function h, and the local material properties of the design (i.e. the homogenized tensor) are approximated by an empirical interpolation law, mimicking physical material properties for the regions with intermediate densities: if the material properties (conductivity, linear elasticity tensor) of D \ Ω and Ω are represented by the respective quantities A 0 , A 1 , the material properties associated to an intermediate density h ∈ [0, 1] are of the form:

(1.1) Several such interpolation schemes and related properties are discussed in [START_REF] Bendsøe | Material interpolation schemes in topology optimization[END_REF] in the setting of the SIMP method for structural mechanics, and more recently in [START_REF] Borrvall | Topology optimization of fluids in stokes flow[END_REF][START_REF] Evgrafov | Topology optimization of slightly compressible fluids[END_REF] in the context of fluid mechanics. In any event, such a relaxation procedure does not alleviate all the difficulties posed by optimal design problems; see [START_REF] Sigmund | Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima[END_REF] for an overview. In the numerical setting, a common issue, related to the homogenized nature of optimal designs, is mesh dependency: using finer and finer computational meshes to increase the numerical accuracy of the Finite Element method often results in designs showing smaller and smaller features, to the point that they become impossible to manufacture. So as to impose a minimum lengthscale, many ideas have been thought off: relying on the theoretical result of [START_REF] Ambrosio | An optimal design problem with perimeter penalization[END_REF], one may impose a constraint on the perimeter of shapes [START_REF] Beckers | Topology optimization using a dual method with discrete variables[END_REF]; other investigations have proposed to enforce a constraint on the gradient or the total variation of the density function [START_REF] Haber | A new approach to variable-topology shape design using a constraint on perimeter[END_REF][START_REF] Petersson | Some convergence results in perimeter-controlled topology optimization[END_REF][START_REF] Petersson | Slope constrained topology optimization[END_REF]. Closer to the present work, another possibility consists in filtering either the sensitivities of the optimization problem [START_REF] Sigmund | On the design of compliant mechanisms using topology optimization[END_REF], or directly the density function h, as proposed in [START_REF] Bruns | Topology optimization of non-linear elastic structures and compliant mechanisms[END_REF] (see [START_REF] Bourdin | Filters in topology optimization[END_REF] for a proof of existence of optimal designs in this context). The idea underlying this last class of remedies is to replace the density function h in (1.1) with a smoothed approximation L ε h obtained with the help of a regularizing operator L ε (the small parameter ε measuring the degree of smoothing). In practice, L ε stands for a convolution with a regularizing kernel [START_REF] Bourdin | Filters in topology optimization[END_REF], or it is obtained by using the regularizing effect of an elliptic equation [START_REF] Lazarov | Filters in topology optimization based on Helmholtz-type differential equations[END_REF]. Interestingly, all the aforementioned techniques also allow to deal with checkerboards, a numerical artifact caused by the use of low-order Finite Element methods.

A 0 + (A 1 -A 0 )ζ(h),
Elaborating on the previous work [START_REF] Amstutz | Connections between topological sensitivity analysis and material interpolation schemes in topology optimization[END_REF] of the first author, and relying on techniques reminiscent of [START_REF] Allaire | Multi-phase optimization via a level set method[END_REF], the purpose of this article is to investigate in which capacity the relaxed and filtered version of an optimal design problem is consistent with the original and more intuitive 'black-and-white' shape and topology optimization problem. To be more specific, we introduce a model problem -featuring an objective criterion J(Ω) involving 'black and white' shapes Ω D -, as well as a relaxed and filtered version, in which designs are density functions h ∈ L ∞ (D, [0, 1]), which are 'filtered' by a regularization operator L ε before being interpolated by a profile ζ accounting for the material properties (1.1) of regions with intermediate densities. This leads to a relaxed, smoothed objective criterion J ε (h) depending on the density function. We then show that the features of this smoothed problem (the value of the objective functional and the optimality conditions) converge to those of the exact shape optimization problem of J(Ω) as the smoothing parameter ε vanishes. Notably, the gradient of J ε (h) (which is defined in the whole computational domain) converges to either the shape derivative of J(Ω) (when restricted to the boundary of Ω) or to its topological gradient (when restricted to inside Ω or D \ Ω).

These results seem interesting for at least two reasons. On the one hand, they prove that density optimization problems may be understood as consistent approximations of 'black-and-white', realistic shape and topology optimization problems, provided the interpolation law for the fictitious material properties is adequately chosen. On the other hand, they give an insight about a common origin for the quite different notions of shape and topological derivatives. At this point, let us mention the article [START_REF] Céa | The shape and topological optimizations connection[END_REF] where another connection between these notions was evidenced, expressing topological derivatives as 'limits of shapes derivatives'.

The remainder of this article is organized as follows. In Section 2, we present the two-phase conductivity equations, a setting under scrutiny in the main part of this work, and we recall some definitions and existing theoretical results as regards the model shape and topology optimization problem of interest in this work. In Section 3, we introduce and analyze the smoothing and relaxation of this problem obtained by considering not only 'black-and-white' shapes as designs, but also density functions. Then, in Sections 4 and 5, we investigate the convergence of this approximate problem to its exact counterpart as the smoothing parameter ε goes to 0. More precisely, the convergence of the approximate conductivity and voltage potential are proved in Section 4, and the convergence of the objective function and that of its derivative are dealt with in Section 5. In Section 6, we describe the extension of the previous material to the linearized elasticity setting. Then, in Section 7, we explain how these results may be used so that the approximate density optimization problem becomes consistent with the exact shape and topology optimization problem as ε → 0, in the sense that the derivative of the approximate objective function with respect to the density function converges to either the shape derivative or the topological derivative of the exact objective function. In Section 8, we present two numerical algorithms for dealing with the two-phase optimization problem, based upon the previous considerations, and some numerical results are discussed in Section 9.

The exact two-phase shape optimization problem in the conductivity setting

In this section, we present the model two-phase conductivity setting considered in the bulk of this article; the targeted applications in linearized elasticity, which are dealt with in Section 6, are natural (albeit sometimes technically tedious) extensions of the considerations in this context which allows for a simpler exposition.

The two-phase conductivity setting.

Let D ⊂ R d be a fixed, smooth bounded domain, and let Ω D be a subset of class C ∞ , with boundary Γ := ∂Ω. Note that a great deal of the material in this article can be adapted to the case where Γ ∩ ∂D = ∅; we ignore such generality for the sake of simplicity. The domain D is filled with two materials with different conductivities 0 < γ 0 = γ 1 < ∞, which respectively occupy the two phases Ω 0 := D \ Ω and Ω 1 := Ω induced by Ω (see Figure 1).

   -div(γ Ω ∇u Ω ) = 0 in D, u Ω = 0 on Γ D , γ 1 ∂uΩ ∂n = g on Γ N .
Introducing the restrictions u i := u Ω | Ω i , recall that (2.2) encompasses transmission conditions: both u Ω and its normal flux γ Ω ∂uΩ ∂n are continuous across the interface Γ between Ω 1 and Ω 0 :

(2.3) [u Ω ] = 0 and γ Ω ∂u Ω ∂n = 0 in the sense of traces in H 1/2 (Γ) and H -1/2 (Γ) respectively,
where n is the unit normal vector to Γ, pointing outward Ω 1 , and [α] := α 0 -α 1 is the jump of a discontinuous quantity α, taking values α 0 , α 1 in Ω 0 , Ω 1 respectively.

Our goal is to optimize the repartition of the two phases Ω 0 and Ω 1 inside D with respect to a mechanical criterion C(Ω), under a constraint on the volume Vol(Ω) = Ω dx of the phase Ω 1 . To set ideas, in our analyses, we limit ourselves with the compliance of the total domain D:

(2.4) C(Ω) = D γ Ω |∇u Ω | 2 dx = Γ N gu Ω ds,
and the volume constraint is incorporated in the problem via a penalization by a fixed Lagrange multiplier > 0. In other terms, the exact shape optimization problem reads:

(2.5) min

Ω∈U ad J(Ω), where J(Ω) = C(Ω) + Vol(Ω),
and U ad is the set of admissible shapes; so as to avoid discussions about the tedious regularity requirements for shape and topological derivative results, we impose that all admissible shapes are 'smooth':

U ad = {Ω D, Ω is of class C ∞ } .
Let us emphasize that extending the conclusions of this article to constrained optimization problems and more general objective functions poses no additionnal difficulty, but only unnecessary technicalities. See Remark 5.3 and the numerical example of Section 9.3 about this point.

Notice also that, contrary to a common custom in the literature, (2.5) is labelled as a 'shape optimization problem'. This terminology is improper since (2.5) leaves room for optimizing the geometry and the topology of Ω, and its purpose is to contrast with the density optimization problem defined in (3.5) below.

The optimization process of the shape and topology of Ω naturally brings into play the sensitivity of J(Ω) with respect to the domain: from the theoretical point of view, it allows to write necessary (local) optimality conditions; from an algorithmic point of view, it paves the way for practical minimization algorithms for J(Ω), e.g. steepest descent methods.

Differentiating J(Ω) with respect to the domain may be achieved in several ways, and we now briefly describe the two complementary ones of interest in this article.

2.2.

Hadamard's boundary variation method and shape derivative.

In order to assess the sensitivity of J(Ω) with respect to variations of the boundary of Ω, we rely on Hadamard's boundary variation method (see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF][START_REF] Murat | Sur le contrôle par un domaine géométrique[END_REF]), whereby variations of a bounded, Lipschitz domain Ω ⊂ R d are considered under the form:

Ω θ := (Id + θ)(Ω), θ ∈ W 1,∞ (R d , R d ), ||θ|| W 1,∞ (R d ,R d ) < 1.
This naturally leads to the following notion of shape differentiability:

Definition 2.1. A function F (Ω) of the domain is shape differentiable at Ω if the underlying function θ → F (Ω θ ), from W 1,∞ (R d , R d ) into R is Fréchet differentiable at θ = 0.
The shape derivative F (Ω) of F at Ω is the corresponding differential, and the following asymptotic expansion holds in the vicinity of

0 ∈ W 1,∞ (R d , R d ): (2.6) F (Ω θ ) = F (Ω) + F (Ω)(θ) + o(θ), where lim θ→0 |o(θ)| ||θ|| W 1,∞ (R d ,R d ) = 0.
In practice, the considered deformations θ are confined to a class Θ ad of admissible deformations θ, so that variations Ω θ of admissible shapes Ω ∈ U ad stay admissible. In the context of this article, one has:

Θ ad = θ ∈ C ∞ (R d , R d ), θ = 0 on ∂D .
Let us now state a classical result about the differentiation of the compliance J(Ω) defined in (2.4), with respect to the domain; see [START_REF] Ito | Variational approach to shape derivatives[END_REF][START_REF] Pantz | Sensibilité de l'équation de la chaleur aux sauts de conductivité[END_REF] for a proof.

Theorem 2.1. The function J(Ω) defined by (2.4) is shape differentiable at any (smooth) shape Ω ∈ U ad , and its shape derivative reads, for an arbitrary deformation θ ∈ Θ ad ,

(2.7) J (Ω)(θ) = Γ g S Ω θ • n ds, where g S Ω := [γ Ω ]∇ Γ u Ω • ∇ Γ u Ω - 1 γ Ω γ Ω ∂u Ω ∂n γ Ω ∂u Ω ∂n + ,
and

∇ Γ f := ∇f -(∇f • n)n denotes the tangential gradient of a regular enough function f : Γ → R.
Remark 2.1. Notice that both terms in the above expression of g S Ω do make sense owing to the transmission conditions (2.3) at the interface Γ.

Topological sensitivity analysis.

A concurrent means to perform variations of a shape Ω consists in subtracting (or adding) a small ball from (to) Ω. Let x ∈ Ω 0 ∪ Ω 1 , and let r > 0 be small enough so that the open ball B(x, r) with center x and radius r does not intersect Γ; we now consider perturbations of Ω of the form:

Ω x,r := Ω \ B(x, r) if x ∈ Ω, Ω ∪ B(x, r) if x ∈ D \ Ω.
Definition 2.2. A function F (Ω) of the domain has a topological derivative at a particular point x ∈ Ω 0 ∪Ω 1 if there exists a real number g T Ω (x) such that the following asymptotic expansion holds: (2.8) s Ω (x) = 1 when x ∈ Ω 0 and s Ω (x) = -1 when x ∈ Ω 1 .

F (Ω r,x ) = F (Ω) + s Ω (x)g T Ω (x)|B(x,
In the context of Section 2.1, the result of interest is the following [10]:

Theorem 2.2. The functional J(Ω) defined by (2.4) has a topological derivative g T Ω (x) at every admissible shape Ω ∈ U ad , and point x ∈ Ω 0 ∪ Ω 1 . Its expression reads:

g T Ω (x) = -(γ 1 -γ 0 )k i ∇u Ω (x) • ∇u Ω (x) + for x ∈ Ω i ,
where the coefficients k i read:

k 0 = 2γ0 γ0+γ1 , k 1 = 2γ1 γ0+γ1 if d = 2 and k 0 = 3γ0 2γ0+γ1 , k 1 = 3γ1 γ0+2γ1 if d = 3.

Local optimality conditions.

Relying on the variations of admissible shapes introduced in Sections 2.2 and 2.3, we are in position to define local minimizers of a function F (Ω) defined over the set U ad of admissible shapes.

Definition 2.3. A shape Ω ∈ U ad is a local minimum for F (Ω) if:
• For any point x ∈ Ω 0 ∪ Ω 1 , and for r > 0 small enough, F (Ω x,r ) ≥ F (Ω).

• For any 'small' vector field θ ∈ Θ ad , F (Ω θ ) ≥ F (Ω).

Using the Definitions 2.1 and 2.2 of shape and topological derivatives, we easily derive the necessary first-order optimality conditions associated to F (Ω).

Proposition 2.3. Let F (Ω) be a function of the domain, and let Ω ∈ U ad be a local minimizer of F (Ω). If F (Ω) has a shape derivative g S Ω at Ω and a topological derivative g T Ω (x) at every point x ∈ Ω 0 ∪ Ω 1 , then: (i) For all x ∈ Ω 1 , g T Ω (x) ≤ 0, (ii) For all x ∈ Γ, g S Ω (x) = 0, (iii) For all x ∈ Ω 0 , g T Ω (x) ≥ 0, In particular, relying on Theorem 2.1 and Theorem 2.2, these considerations immediately apply to the minimization problem (2.5) under scrutiny.

Remark 2.2. Section 2.3 and the present one have only considered spherical inclusions when it comes to the notion of topological sensitivity and the inferred optimality conditions for shape and topology optimization problems. Actually, Definition 2.2 and Theorem 2.2 may be generalized to encompass inclusions of a much more general kind (see e.g. [START_REF] Ammari | Polarization and Moment Tensors; With Applications to Inverse Problems and Effective Medium Theory[END_REF]), and one could for instance be tempted to impose optimality of shapes with respect to topological perturbations caused by elliptic inclusions. As exemplified in [START_REF] Amstutz | Connections between topological sensitivity analysis and material interpolation schemes in topology optimization[END_REF], doing so gives rise to too stringent optimality conditions, and the local minimizers for the relaxed optimization problem tend to show large areas with intermediate densities.

The approximate density optimization problem

As was originally evidenced in the pioneering work of Murat and Tartar [START_REF] Murat | Calculus of Variations and Homogenization[END_REF], one efficient way to cope with many theoretical and numerical difficulties inherent to shape optimization problems of the form (2.5) is the so-called relaxation process: the main idea is to enlarge the set of admissible designs (typically to include density functions), so that it enjoys nicer theoretical properties (in terms of convexity, or compactness notably), or lends itself to simpler or more efficient optimization algorithms. Such a relaxation may be performed rigorously, by the homogenization method [START_REF] Allaire | Shape optimization by the homogenization method[END_REF], or by simpler and more heuristic means, such as the celebrated SIMP method; see [START_REF] Bendsøe | Generating optimal topologies in structural design using a homogenization method[END_REF] or the monograph [START_REF] Bendsøe | Topology Optimization. Theory, Methods, and Applications[END_REF] and references therein.

In this spirit, we presently describe an approximate setting that reformulates the shape and topology optimization problem presented in Section 2 in terms of density functions.

Reformulation of the shape and topology optimization problem as a density optimization problem.

The proposed relaxation of (2.5) consists in replacing the characteristic function χ Ω of the optimized shape Ω involved in the definition (2.1) of the conductivity γ Ω by an interpolated, smoothed density function: for a given h ∈ L 2 (D), we define the (regularized) conductivity γ h,ε by:

(3.1) γ h,ε := γ 0 + (γ 1 -γ 0 )ζ(L ε h) ∈ L ∞ (D).
This procedure relies on two key ingredients:

• The bounded, 'regularizing' operator

(3.2) L ε : L 2 (D) → L ∞ (D),
is self-adjoint when it is seen as an operator from L 2 (D) into itself: L * ε = L ε . Several constructions are possible when it comes to L ε , with different theoretical and numerical assets. In this article, L ε will be one of the two operators L conv

(3.3) ζ, ζ and ζ ∈ L ∞ (R), -ρ ≤ ζ ≤ 1 + ρ, and
ζ(0) = 0, ζ(1) = 1.
Its definition involves a fixed real parameter ρ > 0, which is chosen sufficiently small so that min(γ 0 , γ 1 ) -ρ|γ 1γ 0 |> 0.

The potential u h,ε associated to γ h,ε is the solution in H 1 Γ D (D) to the system:

(3.4)    -div(γ h,ε ∇u h,ε ) = 0 in D, u h,ε = 0 on Γ D , γ h,ε ∂u h,ε ∂n = g on Γ N .
Our relaxed version of (2.5), hereafter referred to as the density approximation problem, is posed over the closed convex subset K := L ∞ (D, [0, 1]) of L 2 (D) composed of density functions:

(3.5) min h∈K J ε (h), where J ε (h) = C ε (h) + Vol ε (h),
the approximate compliance C ε (h) is defined by

(3.6) ∀h ∈ L 2 (D), C ε (h) = D γ h,ε |∇u h,ε | 2 dx = Γ N gu h,ε dx,
and the approximate volume reads:

∀h ∈ L 2 (D), Vol ε (h) = D L ε h dx.
When it comes to the differentiability of J ε (h), the result of interest is the following.

Proposition 3.1. The approximate compliance J ε (h) defined by (3.6) is of class C 1 on L 2 (D), and its Fréchet derivative in an arbitrary direction h ∈ L 2 (D) reads:

(3.7) J ε (h)( h) = D g h,ε h dx, where g h,ε = -(γ 1 -γ 0 )L ε ζ (L ε h)|∇u h,ε | 2 + L ε 1.
Proof. Even though results of this type are standard in optimal control theory, we outline the proof for convenience. Let us first discuss the differentiability of J ε (h): clearly, the mapping

L 2 (D) h → γ h,ε ∈ L ∞ (D) is of class C 1 .
Then, a classical use of the implicit function theorem (see e.g. [START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF], Th. 5.3.2) reveals that the mapping

L 2 (D) h → u h,ε ∈ H 1 Γ D (D) is also of class C 1 .
Therefore, h → J ε (h) has C 1 regularity. We then calculate the Fréchet derivative of J ε (h) by relying on Céa's method [START_REF] Céa | Conception optimale ou identification de formes, calcul rapide de la dérivée directionnelle de la fonction coût[END_REF], which is not merely formal in this case but perfectly rigorous since we already know that h → u h,ε and h → J ε (h) are differentiable. Define the Lagrangian L :

L 2 (D) × H 1 Γ D (D) 2 by: L(h, u, p) = D γ h,ε |∇u| 2 dx + D L ε h dx + D γ h,ε ∇u • ∇p dx - Γ N gp dx.
For a particular h ∈ L 2 (D), let us search for the critical points (u, p) ∈ H

∀ u ∈ H 1 Γ D (D), 2 D γ h,ε ∇u h,ε • ∇ u dx + D γ h,ε ∇p • ∇ u dx = 0,
whence we readily identify: p = -2u h,ε . Now, we use the fact that:

∀h ∈ L 2 (D), ∀p ∈ H 1 Γ D (D), L(h, u h,ε , p) = J ε (h)
; differentiating with respect to h in an arbitrary direction h then evaluating at p = -2u h,ε finally leads to:

J ε (h)( h) = - D ∂γ h,ε ∂h ( h)|∇u h,ε | 2 dx + D L ε h dx,
and the desired result follows after a straightforward calculation.

Remark 3.1. The choice of a self-adjoint operator L ε is a simple commodity in the analysis; the results of this section are easily adapted to the case where this hypothesis is omitted. In particular Proposition 3.1 holds, up to the change in the expression (3.7) of the approximate gradient g h,ε :

(3.8) g h,ε = -(γ 1 -γ 0 )L * ε ζ (L ε h)|∇u h,ε | 2 + L * ε 1.
On the one hand, (3.5) is a relaxation of (2.5) since the set of admissible shapes Ω ∈ U ad is extended to the larger set of density functions h ∈ K. On the other hand, (3.5) is also a smoothing of (2.5) since the values of the conductivity and the objective function of the approximate problem when the design h is the characteristic function χ Ω of an admissible shape Ω are smooth approximations of γ Ω and J(Ω) respectively.

In the following, we shall be interested in the smoothed, non relaxed version of (2.5) obtained by considering only characteristic functions χ Ω of shapes Ω ∈ U ad in (3.5). In this perspective, we denote by χ Ω,ε := L ε χ Ω , γ Ω,ε := γ χΩ,ε , and u Ω,ε := u χΩ,ε , the corresponding (smoothed) characteristic function, conductivity and potential. Likewise, we denote by C ε (Ω), Vol ε (Ω), J ε (Ω) and g Ω,ε the smoothed compliance C ε (χ Ω ), volume Vol ε (χ Ω ), objective J ε (χ Ω ) and gradient g χΩ,ε . We now define the smoothed shape optimization problem by: (3.9) min Ω∈U ad J ε (Ω).

Two examples of constructions of a regularizing operator

L ε : L 2 (D) → L ∞ (D).
In this section, we present the two instances L conv ε and L ell ε of the regularizing operator L ε : L 2 (D) → L ∞ (D) that shall be considered throughout the article. (3.10)

• η has compact support in the unit ball B(0, 1) of R d , • η is nonnegative: η ≥ 0, • η has integral 1, i.e. R d η dx = 1, • η is radial, that is, η(x) = f (|x|)
for some (smooth) function f . It is then well-known that the family of functions η ε , defined by:

(3.11) ∀x ∈ R d , η ε (x) = 1 ε d η
x ε is an approximation of the identity (or a mollifier) in the sense that: Another possible construction of a bounded operator from L 2 (D) into L ∞ (D) relies on the regularizing effects of elliptic partial differential equations: for h ∈ L 2 (D), we denote by L ell ε h the unique solution q ∈ H 1 (D) to the system:

η ε ε→0 -→ δ in D (R d ),
(3.12) -ε 2 ∆q + q = h in D, ∂q ∂n = 0 on ∂D.
By standard elliptic regularity, L ell ε h actually belongs to H 2 (D); since the spatial dimension is d = 2 or 3, Sobolev's embedding theorem (see [START_REF] Adams | Sobolev spaces[END_REF]) implies that H 2 (D) is continuously embedded in the space C(D) of continuous functions up to the boundary ∂D, and therefore in L ∞ (D). Hence, L ell ε does define a bounded operator from L 2 (D) into L ∞ (D), which is easily seen to be self-adjoint from L 2 (D) into itself. Remark 3.2. Many other choices are available as far as the regularizing equation (3.12) for L ell ε is concerned. For instance, one may define

L ell ε h = q(ε, •), where q ∈ L 2 (0, T, H 1 (D)) ∩ H 1 (0, T, (H 1 (D))
) is the unique solution to the heat equation:

   ∂q ∂t -∆q = 0 on (0, T ) × D, ∂q ∂n (t, x) = 0 for (t, x) ∈ (0, T ) × ∂D, q(t = 0, •) = h on D,
a system in which h acts as an initial condition, and T > 0 is a fixed final time.

For further purpose, we now discuss an integral representation of L ell ε . Let us first recall (see [START_REF] Kythe | Fundamental Solutions for partial differential operators, and applications[END_REF], §2.3) that the function G ε (x, y) defined by

(3.13) ∀x = y, G ε (x, y) = 1 ε d G |x -y| ε , where G(r) = 1 (2π) d 2 1 r d 2 -1 K d 2 -1 (r),
and K α is the modified Bessel function of second type with parameter α ≥ 0, is the fundamental solution of the partial differential operator q → -ε 2 ∆q + q in the whole space R d : for a given y ∈ R d , G ε (•, y) satisfies:

-ε 2 ∆ x G ε (x, y) + G ε (x, y) = δ x=y in D (R d ).
In the above identity, δ x=y ∈ D (R d ) is the Dirac distribution centered at y. We now introduce a correction R ε for G ε in such a way that N ε (x, y) := G ε (x, y) + R ε (x, y) becomes the fundamental solution of the operator q → -ε 2 ∆q + q with Neumann boundary conditions on ∂D, i.e. it satisfies:

-

ε 2 ∆ x N ε (x, y) + N ε (x, y) = δ x=y in D, ∂Nε ∂nx (x, y) = 0 on ∂D.
For a given point y / ∈ ∂D, R ε (•, y) is the unique solution in H 1 (D) to the following system:

(3.14) -ε 2 ∆ x R ε (x, y) + R ε (x, y) = 0 for x ∈ D, ∂Rε ∂nx (x, y) = -∂Gε ∂nx (x, y) for x ∈ ∂D,
which makes sense, since for y / ∈ ∂D, the function x → G ε (x, y) is smooth on a neighborhood of ∂D. Now, for h ∈ L 2 (D), extended by 0 to R d \ D, the solution q = L ell ε h to (3.12) reads:

L ell ε h(x) = R d N ε (x, y) h(y) dy.
For further use, the following lemma collects some properties of the kernels G ε , N ε and R ε .

Lemma 3.2.

(i) The correction R ε asymptotically vanishes far from the boundary ∂D as ε → 0: let V be an open neighborhood of Ω such that Ω V D; then:

(3.15) ||R ε (•, y)|| L ∞ (D) ε→0 -→ 0, uniformly in y ∈ V.
(ii) The kernel G satisfies:

(3.16) R d G(|x|) dx = 1.
Proof. (i). Let us first recall the following asymptotic estimates of the K α (see [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF], §9.6.8, §9.6.9 and §9.7.2):

(3.17) K α (r) ∼ r→0 1 2 Γ(α)( 1 2 r) -α if α > 0, K 0 (r) ∼ r→0 -log r , and K α (r) ∼ r→+∞ π 2r e -r ,
where Γ stands for the Euler Gamma function. If V is an open set such that: Ω V D, the decay property of K α (r) as r → ∞ featured in (3.17) combined with standard elliptic regularity estimates for the system (3.14) lead to (3.15).

(ii).

A switch to polar coordinates yields:

R d G(|x|) dx = s d ∞ 0 r d-1 G(r) dr = s d (2π) d/2 ∞ 0 r d/2 K d 2 -1 (r) dr,
where

s d = 2π d/2 Γ( d 2 )
is the measure of the (d -1)-dimensional sphere. Using the formula (see [START_REF] Abramowitz | Handbook of mathematical functions: with formulas, graphs, and mathematical tables[END_REF], §11.4.22):

∞ 0 r α-1 K ν (r) dr = 2 α-2 Γ( α -ν 2 )Γ( α + ν 2 ), for α > ν,
we obtain (3.16).

Remark 3.3. Formula (3.16) may be proved in the following alternative way: taking advantage of the fact that q : x → G(|x|) is the solution of -∆q + q = δ in the sense of distributions in R d , the Fourier transform q(ξ)

:= R d q(x)e -2iπx•ξ dx of q equals q(ξ) = 1 1+4π 2 |ξ| 2 . Hence, R d G(|x|) dx = q(0) = 1.

Optimality conditions for the approximate problem and orientation.

In a nutschell, at this point we have built two approximate problems for (2.5):

• The first one (3.5) is a relaxation and a smoothing of (2.5). The set of designs is enlarged from that of 'black-and-white' shapes Ω ∈ U ad -or equivalently the corresponding characteristic functions, taking only values 0 and 1 -to that of density functions h ∈ K, taking intermediate, 'grayscale' values between 0 and 1. These density functions appear in a smoothed version via the operator L ε . • The second one (3.9) is the restriction of (3.5) to the set of characteristic functions of admissible shapes, and is only a smoothing of (2.5): the optimized designs are still shapes Ω ∈ U ad , but the dependence of the minimization criterion on Ω is smoothed via L ε . These problems feature designs of different natures but both are easier to handle than (2.5) from the theoretical point of view -for instance, one proves along the line of [START_REF] Bourdin | Filters in topology optimization[END_REF] that there exists a minimizer to (3.5) -as well as from the numerical one -the smoothing effect induced by L ε is bound to overcome the mesh-dependency and checkerboards issues mentionned in the introduction. We shall see in Section 8 that they give rise to two different, interesting numerical algorithms, with competing assets and drawbacks.

For the moment, we aim at evaluating in which capacity (3.5) and (3.9) are 'close' approximations to (2.5). This comparison relies on the optimality conditions for (3.9):

Proposition 3.3. Let Ω ∈ U ad be a local minimum for J ε (Ω) in the sense of Definition 2.3; then: (i) For all x ∈ Ω 1 , g Ω,ε (x) ≤ 0, (ii) For all x ∈ Γ, g Ω,ε (x) = 0, (iii) For all x ∈ Ω 0 , g Ω,ε (x) ≥ 0, Proof. Let x ∈ Ω 1 ;
since Ω is a local minimum for J ε (Ω), one has, for r > 0 small enough:

J ε (Ω x,r ) = J ε (χ Ω,ε -χ B(x,r) ) ≥ J ε (Ω).
Since χ B(x,r) → 0 strongly in L 2 (D) as r → 0, using Proposition 3.1, we obtain that g Ω,ε (x) ≤ 0, which proves (i).

Let us now consider a point x ∈ Γ; for r > 0 small enough, one has:

J ε (Ω ∪ B(x, r)) = J ε (χ Ω + χ c Ω∩B(x,r) ) ≥ J ε (Ω).
Using again Proposition 3.1 together with the fact that χ c Ω∩B(x,r) → 0 strongly in L 2 (D) as r → 0, we obtain that g Ω,ε (x) ≥ 0. The same argument, using now

J ε (Ω ∩ c B(x, r)) = J ε (χ Ω -χ Ω∩B(x,r) ) ≥ J ε (Ω)
yields the converse inequality g Ω,ε (x) ≤ 0, which proves (ii). Eventually, (iii) is proved in the same way as (i).

In the next two Sections 4 and 5, we shall prove that, under suitable assumptions on the interpolation profile ζ involved in (3.1), for a given admissible shape Ω ∈ U ad , the objective and the optimality conditions of the approximate shape optimization problem (3.9) 'converge' to their counterparts in the context of the exact shape optimization problem (2.5). More precisely, in Section 4, we investigate the convergence of the approximate conductivities γ Ω,ε and voltage potential u Ω,ε to their exact counterparts γ Ω and u Ω for a given shape Ω ∈ U ad . In Section 5, we then study the convergence of J ε (Ω) to J(Ω) and that of the optimality conditions of Proposition 3.3 to those in Proposition 2.3.

Asymptotic behaviors of the approximate conductivity γ Ω,ε and potential u Ω,ε

In this section, we examine the convergence of the approximate conductivity γ Ω,ε and potential u Ω,ε featured in the approximate shape optimization problem (3.9) towards their respective counterparts γ Ω and u Ω .

Let us first recall a few useful facts. Since, in the setting of Section 2, Ω D is a bounded domain of class C 2 with boundary Γ, the projection mapping is well-defined and of class C 1 on some open neighborhood V of Γ; see e.g. [START_REF] Delfour | Shapes and Geometries: Metrics, Analysis, Differential Calculus, and Optimization[END_REF]. Hence, the unit normal vector n to Γ, pointing outward Ω, may be extended into a unit vector field (still denoted by n) on V via the formula:

∀x ∈ V, n(x) = n(p Γ (x)).
By the same token, an arbitrary tangential vector field τ : Γ → R d has a natural extension to V , and so do the notions of normal derivative ∂v ∂n and tangential gradient ∇ Γ v of a smooth function v : V → R.

4.1.

Asymptotic behavior of the approximate conductivity γ Ω,ε .

Our purpose in this section is to prove the following proposition:

Proposition 4.1. The following convergence results hold:

(i) χ Ω,ε → χ Ω and γ Ω,ε → γ Ω strongly in L p (D) as ε → 0, for 1 ≤ p < ∞, (ii) For any open set U Ω 0 ∪ Ω 1 , and for any n ≥ 1, γ Ω,ε → γ Ω strongly in H n (U ), (iii)
There exists an open neighborhood V of Γ such that for any tangential vector field τ : Γ → R d , ∂γΩ,ε ∂τ → ∂γΩ ∂τ = 0 strongly in L p (V ), for 1 ≤ p < ∞. These properties stem from quite different arguments depending on whether the regularizing operator In this subsection, γ Ω,ε is defined by the formula:

L ε : L 2 (D) → L ∞ (D) at play is that L conv ε defined in Section 3.2.1 or that L ell ε defined in Section 3.2.2.
γ Ω,ε = γ 0 + (γ 1 -γ 0 )ζ(χ Ω,ε ), where χ Ω,ε = η ε * χ Ω ,
and η ε is the smoothing kernel given by (3.10) (3.11). Let us start with a fairly easy lemma about χ Ω,ε : ∂χ

Ω,ε ∂n L ∞ (Γε) ≤ C ε , and ∂χ Ω,ε ∂τ L ∞ (Γε) ≤ C,
for a constant C which only depends on Ω and not on ε.

Proof. Since Ω is of class C 2 , an application of Green's formula reveals that the derivative of the characteristic function χ Ω in the sense of distributions reads:

∇χ Ω = -δ Γ n,
where δ Γ ∈ D (R d ) is the surface measure distribution of Γ, defined by:

∀ϕ ∈ D(R d ), δ Γ , ϕ = Γ ϕ ds,
where •, • is the duality pairing between D (R d ) and D(R d ). For an arbitrary test function

ϕ ∈ D(R d ), it comes: ∂χ Ω,ε ∂n , ϕ = η ε * (∇χ Ω ), ϕn = - Γ n(x) • η ε * (ϕn)(x) ds(x),
and, using elementary manipulations, we identify:

∂χ Ω,ε ∂n (x) = - Γ η ε (x -y) n(x)
• n(y) ds(y), and, likewise

∂χ Ω,ε ∂τ (x) = - Γ η ε (x -y) τ (x) • n(y) ds(y). Now, we estimate, for x ∈ Γ ε , ∂χ Ω,ε ∂n (x) = 1 ε d Γ η( x -y ε )n(x) • n(y) ds(y) ≤ C ε d Γ∩B(x,ε) ds(y), ≤ C ε ,
where the constant C only depends on the mollifier η and the properties of the surface Γ.

The corresponding estimate for ∂χΩ,ε ∂τ follows in the same way, using the additional ingredient that there exists a constant C such that for any two points x, y

∈ Γ ε , |x-y|< ε ⇒ |τ (x)•n(y)|= |τ (x)•(n(y)-n(x))|≤ Cε. Proof of Proposition 4.1 in the case L ε = L conv ε . (i).
As a well-known property of mollifiers, for any 1 ≤ p < ∞, χ Ω,ε → χ Ω , strongly in L p (D). Note that, since the χ Ω,ε are continuous and χ Ω is not, such a convergence cannot possibly hold in L ∞ (Ω). (i) is therefore an immediate consequence of the fact that ζ is of class C 1 .

(ii). This point is immediate since, for ε small enough, γ Ω,ε = γ Ω on U .

(iii). This point is a straighforward consequence of Lemma 4.2.

The case where L

ε = L ell ε .
In this section, we assume that γ Ω,ε reads:

γ Ω,ε = γ 0 + (γ 1 -γ 0 )ζ(χ Ω,ε ),
where χ Ω,ε is the unique solution in H 1 (D) to the system:

(4.2) -ε 2 ∆χ Ω,ε + χ Ω,ε = χ Ω in D, ∂χΩ,ε ∂n = 0 on ∂D.
Proof of Proposition 4.1 in the case L ε = L ell ε . (i). We start by showing that χ Ω,ε → χ Ω strongly in L 2 (D). The variational formulation for (4.2) reads:

(4.3) ∀v ∈ H 1 (D), ε 2 D ∇χ Ω,ε • ∇v dx + D χ Ω,ε v dx = D χ Ω v dx.
Hence, we obtain the a priori estimate:

(4.4) ε 2 ||∇χ Ω,ε || 2 L 2 (D) d +||χ Ω,ε || 2 L 2 (D) ≤ ||χ Ω || 2 L 2 (D) .
Therefore, for a fixed test function

v ∈ C ∞ c (R d ), D χ Ω,ε v dx ε→0 -→ D χ Ω v dx.
Owing to a standard density argument, the above convergence actually holds for arbitrary v ∈ L 2 (D), so that χ Ω,ε → χ Ω weakly in L 2 (D) as ε → 0. Now taking limits in (4.4) yields:

||χ Ω || L 2 (D) ≤ lim inf ε→0 ||χ Ω,ε || L 2 (D) ≤ ||χ Ω || L 2 (D) ,
which reveals that the convergence χ Ω,ε → χ Ω holds strongly in L 2 (D) (and thus in L p (D) for 1 ≤ p ≤ 2). On a different note, the maximum principle for (4.2) (see e.g. [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], §9.7) implies that:

||χ Ω,ε || L ∞ (D) ≤ ||χ Ω || L ∞ (D) .
Now, for an arbitrary 2 ≤ p < ∞, we simply calculate:

D |χ Ω,ε -χ Ω | p dx ≤ ||χ Ω,ε -χ Ω || p-2 L ∞ (D) ||χ Ω,ε -χ Ω || 2 L 2 (D) ; therefore χ Ω,ε → χ Ω strongly in L p (D) for 1 ≤ p < ∞, whence (i) easily follows.
(ii). It is an immediate consequence of the standard (interior) regularity theory for elliptic equations (and of the fact that the right-hand side in (4.2) is smooth on U ). See for instance [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], §9.6.

(iii). The proof is similar to that of the first point, up to a localization and local charts argument analogous to that used in the proof of Theorem 4.5 performed in Appendix A (which arises in a slightly more difficult context), and we do not repeat the argument.

4.2.

Asymptotic behavior of the approximate potential u Ω,ε .

In this section, we analyze the asymptotic properties of the approximate potential u Ω,ε and of its tangential and normal fluxes; these are quite independent of the nature of the regularizing operator L ε involved, as long as γ Ω,ε satisfies the conclusions of Proposition 4.1, which we assume throughout this section.

Let us first recall Meyer's theorem (see [START_REF] Bensoussan | Asymptotic analysis of periodic structures[END_REF], p. 38):

Theorem 4.3. Let D ⊂ R d be a smooth bounded domain whose boundary ∂D is decomposed as:

∂D = Γ D ∪ Γ N ; let 0 < m < M < ∞ be fixed. Let A ∈ L ∞ (D
) d×d be a matrix-valued function such that:

∀ξ ∈ R d , m|ξ| 2 ≤ A(x)ξ • ξ ≤ M |ξ| 2 , a.e. in x ∈ D.
Then, there exists p > 2 and a constant C, depending only on Ω, m and M such that, for any right-hand side f ∈ W -1,p (D), the unique solution u ∈ H 1 Γ D (D) to the system:

   -div(A∇u) = f in D, u = 0 on Γ D , (A∇u) • n = 0 on Γ N ,
actually belongs to W 1,p (D) and satisfies:

||u|| W 1,p (D) ≤ C||f || W -1,p (D) .
Remark 4.1. Meyer's theorem is a local result; it has numerous generalizations, notably to different types of boundary conditions, and to the case of systems (e.g. the linearized elasticity system).

As far as the convergence of the potential u Ω,ε is concerned, the result of interest is the following. Proposition 4.4. Let Ω ∈ U ad , and assume that γ Ω,ε satisfies Proposition 4.1. Then, u Ω,ε → u Ω strongly in W 1,p (D) for some p > 2.

Proof. The variational formulations associated to (2.2) and (3.4) read: for an arbitrary

v ∈ H 1 Γ D (D), D γ Ω ∇u Ω • ∇v dx = Γ N gv dx, D γ Ω,ε ∇u Ω,ε • ∇v dx = Γ N gv dx.
Subtracting and rearranging yields:

(4.5) ∀v ∈ H 1 Γ D (D), D γ Ω,ε ∇(u Ω,ε -u Ω ) • ∇v dx = D (γ Ω -γ Ω,ε )∇u Ω • ∇v dx.
Taking v = u Ω,εu Ω in the previous identity, we obtain:

D γ Ω,ε |∇(u Ω,ε -u Ω )| 2 dx = D (γ Ω -γ Ω,ε )∇u Ω • ∇(u Ω,ε -u Ω ) dx.
Now, using Meyer's theorem, there exists p > 2 such that u Ω ∈ W 1,p (D), and a constant C such that

||u Ω || W 1,p (D) ≤ C||g|| H 1/2 (Γ N ) .
Let q < ∞ be the real number defined by the relation 1 p + 1 q + 1 2 = 1. An application of Hölder's inequality yields the existence of a constant C such that:

||u Ω,ε -u Ω || H 1 (D) ≤ C||γ Ω,ε -γ Ω || L q (D) ||g|| H 1/2 (Γ N ) ,
which proves the desired strong convergence result in H 1 (D) norm. The proof that this convergence holds in W 1,p (D) for some p > 2 follows the same path, applying Meyer's Theorem 4.3 from the identity (4.5)

We now turn to study the convergence of the tangential derivative and the normal flux of u Ω,ε . The proof of this result is postponed to Appendix A. Theorem 4.5. For a given admissible shape Ω ∈ U ad , (i) Let U be an open set such that U Ω 0 or U Ω 1 . Then, for any n ≥ 1, the following convergence holds:

u Ω,ε → u Ω strongly in H n (U ). (ii) Let V be a neighborhood of Γ such that the projection p Γ given by (4.1) is well-defined on V . Then

(4.6) ∂u Ω,ε ∂τ → ∂u Ω ∂τ strongly in H 1 (V ), and γ Ω,ε ∂u Ω,ε ∂n → γ Ω ∂u Ω ∂n strongly in H 1 (V ).
Remark 4.2. A closer look to the proof of Theorem 4.5 reveals that the convergence results (4.6) actually hold strongly in W 1,p (V ) for some p > 2, but we shall not need this fact.

Convergence of the approximate objective function and its derivative

In this section, we investigate the convergence of the objective J ε (χ Ω ), and that of the derivative g Ω,ε defined in (3.7), for a given shape Ω ∈ U ad .

First convergence results.

The convergence of the approximate objective J ε (Ω) towards its exact counterpart J(Ω) straightforwardly follows from Proposition 4.4.

Proposition 5.1. Let Ω ∈ U ad ; then J ε (Ω) ε→0 -→ J(Ω).
Let us now turn to the convergence of the approximate gradient g Ω,ε inside each phase Ω 0 and Ω 1 . Proposition 5.2. Let Ω ∈ U ad . Then, for any 1 ≤ p < ∞,

g Ω,ε ε→0 -→ -(γ 1 -γ 0 )ζ (0)|∇u Ω | 2 + strongly in L p (Ω 0 ), -(γ 1 -γ 0 )ζ (1)|∇u Ω | 2 + strongly in L p (Ω 1 ).
Proof. Let x 0 ∈ Ω 0 ∪ Ω 1 , and r > 0 be small enough so that the ball B(x 0 , r) does not intersect Γ; let also n ≥ 1 be a large enough integer so that the continuous embedding H n (B(x 0 , r)) ⊂ C 1 (B(x 0 , r)) holds.

Recalling the expression (3.7) of g Ω,ε , Proposition 4.1 and Theorem 4.5 (i), one obtains, for 1 ≤ p < ∞:

g Ω,ε ε→0 ---→ -(γ 1 -γ 0 )ζ (0)|∇u Ω | 2 + strongly in L p (B(x 0 , r)) if x 0 ∈ Ω 0 , -(γ 1 -γ 0 )ζ (1)|∇u Ω | 2 + strongly in L p (B(x 0 , r)) if x 0 ∈ Ω 1 .
Now recalling the expression of the approximate conductivity γ Ω,ε and the two possible constructions L conv ε and L ell ε carried out in Section 3.2, it is easily seen that g Ω,ε is bounded in L ∞ (D) independently of ε. A use of Lebesgue's dominated convergence theorem allows to conclude. Proposition 5.2 really supplies information about the convergence of the approximate gradient g Ω,ε (x) at points x lying in the interior of Ω 0 or Ω 1 , and we shall make the connection between its limit and the topological derivative g T Ω of the exact shape optimization problem (2.5) in Section 7. However, the lack of uniform convergence of g Ω,ε in the neighborhood of Γ indicates that it has a quite different asymptotic behavior in this region, which we now investigate.

Convergence of the approximate gradient on the boundary Γ.

The main result of this section is the following: Theorem 5.3. In the setting of Sections 2 and 3, the restriction to Γ of the approximate gradient g Ω,ε at an admissible shape Ω ∈ U ad converges to the exact shape derivative g S Ω given by (2.7):

g Ω,ε ε→0 -→ g S Ω in L 1 (Γ).
Remark 5.1. The conclusion of Theorem 5.3 is quite remarkable insofar as the information contained in the approximate gradient g Ω,ε is by essence of a topological nature: loosely speaking, g Ω,ε does not distinguish the bulk and the boundary of Ω.

The proof of Theorem 5.3 mainly consists of two technical lemmas.

ψ ε = L ε (ϕ (L ε χ Ω )) .
Then, for any point x 0 ∈ Γ, one has:

(5.1) lim ε→0 ψ ε (x 0 ) = ϕ(1) -ϕ(0). Remark 5.2. Since ψ ε is bounded in L ∞ (Γ)
uniformly with respect to ε, Lebesgue's dominated convergence theorem implies that the above convergence actually holds in L p (Γ), 1 ≤ p < ∞.

Proof of Lemma 5.4. The proof relies on completely similar arguments in both situations L ε = L conv ε and L ε = L ell ε , and we deal with it in a joint fashion insofar as possible. To achieve this, for any density h ∈ L 2 (D) (extended by 0 outside D), we write the function L ε h ∈ C(D) under the form:

L ε h(x) = R d N ε (x, y)h(y) dy, with N ε (x, y) = G ε (x, y) + R ε (x, y).
In this notation, G ε (x, y) has the structure

G ε (x, y) = 1 ε d G |x -y| ε ,
where the (nonnegative) real function G : (0, ∞) → R is smooth, satisfies R d G(|x|) dx = 1 and the following decay conditions: there exists constants C > 0 and M > 0 such that:

(5.2) G(r) ≤ C|log r| if d = 2, G(r) ≤ C r d-2 if d ≥ 3, for r ≤ 1 and G(r) ≤ Ce -M r for r ≥ 1;
the correction R ε (x, y) satisfies:

(5.3) For any open set V D, ||R ε (•, y)|| L ∞ (D) ε→0 -→ 0 uniformly in y ∈ V.
As we have seen in Sections 3.2.1 and 3.2.2, both situations L ε = L conv ε and L ε = L ell ε fall into this setting with:

• In the case

L ε = L conv ε , G(r) = η(r) and R ε (x, y) = 0. • In the case L ε = L ell ε , G(r) is defined in (3.13
) and R ε (x, y) is the solution to the system (3.14). According to the above representation, one has, for x ∈ Γ:

ψ ε (x) = R d N ε (x, y) ϕ R d N ε (y, z)χ Ω (z) dz dy.
We now divide the proof of the lemma into three steps.

Step 1: Elimination of the correction induced by R ε (x, y).

Define the scalar function ψ

ε : Γ → R as:

ψ (1) ε (x) = R d G ε (x, y) ϕ R d G ε (y, z)χ Ω (z) dz dy.
We prove that, for a given x 0 ∈ Γ:

(5.4)

|ψ ε (x 0 ) -ψ (1) ε (x 0 )| ε→0 -→ 0.
To this end, using that

||R ε (x 0 , •)|| L ∞ (D) ε→0
-→ 0 together with the fact that ϕ is bounded, we see first that:

|ψ ε (x 0 ) -ψ ε (x 0 )| ε→0 -→ 0, where ψ ε (x 0 ) := R d G ε (x 0 , y) ϕ R d N ε (y, z)χ Ω (z) dz dy.
Now let r > 0 be so small that B(x 0 , r) D. We decompose:

ψ ε (x 0 ) = B(x0,r) G ε (x 0 , y) ϕ R d N ε (y, z)χ Ω (z) dz dy + R d \B(x0,r) G ε (x 0 , y) ϕ R d N ε (y, z)χ Ω (z) dz dy.
Since ϕ is bounded, the second term in the right-hand side vanishes in the limit ε → 0 as a consequence of the decay property (5.2). As for the first one, using the fact that:

||R ε (y, •)|| L ∞ (D) ε→0 -→ 0 uniformly in y ∈ B(x 0 , r),
it follows:

R d N ε (y, z)χ Ω (z) dz - R d G ε (y, z)χ Ω (z) dz ε→0 -→ 0, uniformly in y ∈ B(x 0 , r)
and since ϕ is Lipschitz and G ε (x 0 , •) is bounded in L 1 (R d ) independently of ε, we obtain:

ψ ε (x 0 ) - B(x0,r) G ε (x 0 , y) ϕ R d G ε (y, z)χ Ω (z) dz dy ε→0 -→ 0.
Eventually, using the convergence

ψ (1) ε (x 0 ) - B(x0,r) G ε (x 0 , y) ϕ R d G ε (y, z)χ Ω (z) dz dy ε→0 -→ 0, which is proved thanks to the fact ||G ε (x 0 , •)|| L ∞ (R d \B(x0,r)) ε→0
-→ 0, we eventually get the sought result (5.4).

Step 2: Reduction to the case where Γ is flat around x 0 .

Precisely, we prove that:

(5.5)

|ψ (1) ε (x 0 ) -ψ (2) ε (x 0 )| ε→0 -→ 0, where ψ (2) ε (x 0 ) := R d G ε (x 0 , y) ϕ R d G ε (y, z)χ(z) dz dy,
and χ is the characteristic function of the half-space

x ∈ R d , (x -x 0 ) • n(x 0 ) < 0 .
Grossly speaking, the proof of (5.5) consists in estimating the 'difference' between this half-space and the domain Ω 'near' the point x 0 , which is the only region in R d where G ε (x 0 , •) takes large values.

To this end, simple changes of variables in the definitions of ψ

(1) ε and ψ

(2) ε yield:

ψ (1) ε (x 0 ) = R d G(|y|)ϕ R d G(|z|) χ Ω (x 0 -εy -εz) dz dy, and
ψ (2) ε (x 0 ) = R d G(|y|)ϕ R d G(|z|) χ(x 0 -εy -εz) dz dy.
To prove (5.5), owing to the Lipschitz character of ϕ , it is therefore enough to show that:

(5.6)

R d R d G(|y|)G(|z|) |χ Ω (x 0 -εy -εz) -χ(x 0 -εy -εz)| dydz ε→0 -→ 0.
Using the decay rate (5.2) it is in turn sufficient to prove that, for any fixed M > 0,

(5.7)

I ε := B(0,M ) B(0,M ) G(|y|)G(|z|)|χ Ω (x 0 -εy -εz) -χ(x 0 -εy -εz)| dydz ε→0 -→ 0,
which we now proceed to do. Using again the decay condition (5.2) for G and Hölder's inequality, we have the rough bound, for some p > 1:

I p ε ≤ C B(0,2M )
|χ Ω (x 0εy)χ(x 0εy)| p dy.

Here and throughout the proof, the constant C represents a generic positive constant which changes from one line to the next, but is independent of ε. Now, without loss of generality, we assume that x 0 = 0 and n(x 0 ) = e d (the d th coordinate vector). Since Ω is of class C 2 , there exists a neighborhood U of x 0 and a function f :

R d-1 → R of class C 2 such that Ω ∩ U is the subgraph of f , that is: Ω ∩ U = {x ∈ U s.t. x d < f (x 1 , ..., x d-1 )} , Γ ∩ U = {x ∈ U s.t. x d = f (x 1 , ..., x d-1 )} .
Hence, we obtain:

I p ε ≤ C [-2M,2M ] d-1 |f ((x 0,1 -εy 1 , ..., x 0,d-1 -εy d-1 ))| dy 1 ...dy d-1 ,
where the right-hand side simply corresponds to the red dashed area on Figure 2. Since f is of class C 2 , f (x 0 ) = 0 and ∇f (x 0 ) = 0, it follows that:

I ε ≤ Cε 2/p .
This terminates the proof of (5.5). Step 3: Reduction to a one-dimensional situation and end of the proof of Lemma 5.4

We eventually prove the limit lim

ε→0 ψ (2)
ε (x 0 ) = ϕ(1)ϕ(0). Without loss of generality, we retain the assumption that x 0 = 0 and n(x 0 ) = e d . Recall the alternative expression for ψ

(2) ε obtained by using changes of variables, (5.8)

ψ (2) ε (x 0 ) = R d G(|y|) ϕ R d G(|z|) χ(x 0 -εy -εz) dz dy.
The key point now consists in remarking that, since χ is the characteristic function of the lower half-space, the integral featured in (5.8)

y → R d G(|z|) χ(x 0 -εy -εz) dz = {z∈R d , -y d ≤z d } G(|z|) dz
depends on y only via its d th -component. To be precise, let us introduce the function G ∈ L 1 (R) given by:

G(s) = R d-1 G(|(z 1 , ..., z d-1 , s)|) dz 1 ...dz d-1 ;
notice this definition is possible owing to (5.2). A repeated use of Fubini's theorem yields:

(5.9)

ψ (2) ε (x 0 ) = R d G(|y|) ϕ +∞ -y d G(s) ds dy, = R d G(|y|) ϕ ( G * χ)(y d ) dy, = R G(s) ϕ ( G * χ)(s) ds.
In the previous inequalities, we have introduced the usual (one-dimensional) Heaviside function χ ∈ L ∞ (R):

χ(s) = 1 if s > 0, 0 otherwise.
Now, the function G * χ has the expression:

( G * χ)(t) = {y∈R d , y d <t} G(|y|) dy = t -∞ G(s) ds, for a.e. t ∈ R.
It is therefore absolutely continuous, with derivative G, and has limits:

lim t→+∞ ( G * χ)(t) = R d G(|y|) dy = 1, and lim t→-∞ ( G * χ)(t) = 0.
It eventually follows from (5.9) that:

ψ (2) ε (x 0 ) = R d ds ϕ G * χ ds = lim t→+∞ ϕ(( G * χ)(t)) -lim t→-∞ ϕ(( G * χ)(t)), = ϕ(1) -ϕ(0),
which completes the proof.

The second ingredient in the proof of Theorem 5.3 is the following:

Lemma 5.5. Let L ε : L 2 (D) → L ∞ (D) stand for one of the two operators L conv ε and L ell ε constructed in Section 3.2. Let v ε , v ∈ H 1 (D) and f ε ∈ L ∞ (D) be such that v ε ε→0 -→ v strongly in H 1 (V ), and ||f ε || L ∞ (D) ≤ C,
for an open neighborhood V of Γ and a constant C > 0 which is independent of ε. Then the following convergence holds:

r ε := L ε (f ε v 2 ε ) -L ε (f ε )v 2 ε ε→0 -→ 0 strongly in L 1 (Γ).
Proof. We retain the notations introduced in the preamble of the proof of Lemma 5.4 in order to deal jointly with both cases L ε = L conv ε and L ε = L ell ε . At first, easy calculations using (5.2) yield, for an arbitrary M > 0:

Γ |r ε (x)| ds(x) ≤ C Γ R d G(|y|)|v 2 ε (x -εy) -v 2 ε (x)| dy ds(x) + o(1), ≤ C Γ B(0,M ) G(|y|)|v 2 ε (x -εy) -v 2 ε (x)| dy ds(x) +C Γ R d \B(0,M ) G(|y|)|v 2 ε (x -εy) -v 2 ε (x)| dy ds(x) + o(1), ≤ C Γ B(0,M ) G(|y|)|v 2 ε (x -εy) -v 2 ε (x)| dy ds(x) + Ce -αM ||v ε || 2 H 1 (V ) +o(1),
where o(1) → 0 as ε → 0. Since v ε is uniformly bounded in H 1 (D), it is enough to prove that, for an arbitrary, fixed M > 0,

I ε := Γ B(0,M ) G(|y|)|v 2 ε (x -εy) -v 2 ε (x)| dy ds(x) ε→0 -→ 0.
To do so, using polar coordinates and (5.2), we obtain:

I ε = Γ M 0 ∂B(0,1) r d-1 G(r)|v 2 ε (x -εrω) -v 2 ε (x)| ds(ω) dr ds(x), ≤ C Γ M 0 ∂B(0,1) |v 2 ε (x -εrω) -v 2 ε (x)| ds(ω) dr ds(x), ≤ C Γ M 0 ∂B(0,1) |v 2 ε (x -εrω) -v 2 (x -εrω)| ds(ω) dr ds(x), +C Γ M 0 ∂B(0,1) |v 2 (x -εrω) -v 2 (x)| ds(ω) dr ds(x) + C Γ M 0 ∂B(0,1) |v 2 ε (x) -v 2 (x)| ds(ω) dr ds(x), ≤ C Γ M 0 ∂B(0,1) |v 2 (x -εrω) -v 2 (x)| ds(ω) dr ds(x) + C||v ε -v|| H 1 (V ) .
Hence, the proof is complete if we can show that, for an arbitrary function v ∈ H 1 (V ), the integral:

Γ M 0 ∂B(0,1) |v 2 (x -εrω) -v 2 (x)| ds(ω) dr ds(x)
converges to 0 as ε → 0; but this last point easily follows from a standard density argument of smooth functions in H 1 (V ).

We are now in good position to complete the proof of Theorem 5.3.

Proof of Theorem 5.3. Let us decompose the expression (3.7) of g Ω,ε as:

g Ω,ε = -(γ 1 -γ 0 ) L ε ζ (L ε χ Ω ) ∂u Ω,ε ∂τ 2 -(γ 1 -γ 0 ) L ε ζ (L ε χ Ω ) γ 2 Ω,ε γ Ω,ε ∂u Ω,ε ∂n 2 + L ε 1.
Let us first remark that, as an easy consequence of Lemma 5.5 and Theorem 4.5:

(5.10)

g Ω,ε -g Ω,ε ε→0 -→ 0 in L 1 (Γ),
where:

(5.11)

g Ω,ε := -(γ 1 -γ 0 ) (L ε ζ (L ε χ Ω )) ∂u Ω,ε ∂τ 2 -(γ 1 -γ 0 ) L ε ζ (L ε χ Ω ) γ 2 Ω,ε γ Ω,ε ∂u Ω,ε ∂n 2 + . Now, using Lemma 5.4 with ϕ(t) = γ 0 + (γ 1 -γ 0 )ζ(t), ϕ (t) = (γ 1 -γ 0 )ζ (t), then ϕ(t) = 1 γ0+(γ1-γ0)ζ(t) , ϕ (t) = -(γ1-γ0)ζ (t) (γ0+(γ1-γ0)(1-ζ(t))) 2 , we obtain: (γ 1 -γ 0 )L ε (ζ (L ε χ Ω )) ε→0 -→ γ 1 -γ 0 a.e. on Γ, (γ 1 -γ 0 )L ε ζ (L ε χ Ω ) γ 2 Ω,ε ε→0 -→ 1 γ 0 - 1 γ 1 a.e. on Γ;
in light of (5.10), (5.11) and Theorem 4.5, this terminates the proof.

Remark 5.3. Let us comment on how the previous material extends to the case where the function C(Ω) involved in (2.5) is not the compliance (2.4).

The main difference between this new case and the previous one lies in the expressions of the shape and topological derivatives g S Ω and g T Ω : as is well-known (see e.g. [START_REF] Allaire | Conception optimale de structures[END_REF][START_REF] Henrot | Variation et optimisation de formes, une analyse géométrique[END_REF]), they now involve an adjoint state p Ω ∈ H 1 Γ D (D), solution to a system of the form (2.2), with a different right-hand side (depending on the expression of J(Ω)). For instance, Formula (2.7) for the shape derivative g S Ω becomes:

g S Ω := -[γ Ω ]∇ Γ u Ω • ∇ Γ p Ω + 1 γ Ω γ Ω ∂u Ω ∂n γ Ω ∂p Ω ∂n + ,
and a similar transformation occurs in the expression of g T Ω . The derivative of the approximate functional J ε (h), cooked from J(Ω) along the lines of Section 3, involves in turn an approximate adjoint state p h,ε ∈ H 1 Γ D (D), solution to a similar system to (3.4), with, again, a different right-hand side.

In this context, the convergence of the approximate objective J ε (Ω) and that of its derivative are studied as in the proofs of Propositions 5.1, 5.2 and Theorem 5.3, with the additional study of the convergence p Ω,ε → p Ω which is carried out in a similar way to that of Theorem 4.5.

Remark 5.4. We have hitherto assumed the 'hold-all' domain D to be smooth; this is merely for technical convenience, and the results of Sections 2, 3, 4 and 5 hold true when D is only Lipschitz regular, provided their conclusions are restricted to a fixed subset D D such that Ω D D.

Extension to the linear elasticity case

In this section, we point out how the proposed approximation process of shape and topology optimization problems extends from the setting of the conductivity equation to that of the linearized elasticity system. Both situations are similar from the mathematical point of view, and we mainly discuss the differences; so as to emphasize the parallel between them, we re-use the notations from the previous sections for their counterparts in the present context insofar as possible.

Description of the linear elasticity setting.

The geometric situation is that of Section 2: the computational domain D is divided into two complementary phases Ω 0 = D \ Ω and Ω 1 = Ω, separated by the interface Γ = ∂Ω, which are filled with two different linearly elastic materials with respective Hooke's law A 0 , A 1 . We assume that both laws are isotropic, i.e. they are of the form: ∀e ∈ S(R d ), A i e = 2µ i e + λ i tr(e)I, i = 0, 1, where I is the identity matrix of size d, and µ i , λ i are the Lamé coefficients of the i th material, which satisfy:

µ i > 0, and λ i + 2µ i /d > 0.
Recall that the bulk modulus κ of such a material equals κ = λ + 2 d µ, and that the properties of a linearly elastic material may be equivalently described in terms of its Young's modulus E and Poisson's ratio ν, defined by:

E = 4κµ κ + µ , ν = κ -µ κ + µ if d = 2, and E = 9κµ 3κ + µ , ν = 3κ -2µ 2(3κ + µ) if d = 3.
The total Hooke's law A Ω inside the structure D is discontinuous, and reads: (6.1)

A Ω = A 0 + (A 1 -A 0 )χ Ω .
The structure D is clamped on a subset Γ D ⊂ ∂D, and surface loads g ∈ H -1/2 (Γ N ) d are applied on the complementary part Γ N = ∂D \ Γ D ; we omit body forces for simplicity. The displacement of D is then the unique solution u Ω ∈ H 1 Γ D (D) d to the system:

(6.2)    -div(A Ω e(u Ω )) = 0 in D, A Ω e(u Ω )n = g on Γ N , u Ω = 0 on Γ D ,
where e(u Ω ) = 1 2 (∇u T Ω + ∇u Ω ) is the strain tensor associated to u Ω , and the corresponding stress tensor σ(u Ω ) is defined by σ(u Ω ) = A Ω e(u Ω ).

As in the case of the conductivity equation (2.2), the system (6.2) encompasses transmission conditions at the interface Γ. To express them, let us introduce, at each point of Γ, a local basis of R d obtained by gathering the unit normal vector n (pointing outward Ω 1 ) and a collection of unit tangential vectors, denoted by τ , such that (τ, n) is an orthonormal frame. For a symmetric d × d matrix M written in this basis, we introduce the notation

M = M τ τ M τ n M nτ M nn
where M τ τ stands for the (d -1) × (d -1) minor of M, M τ n is the vector of the (d -1) first components of the d-th column of M, M nτ is the row vector of the (d -1) first components of the d-th row of M, and M nn is the (d, d) entry of M.

Then it is well-known that u Ω is continuous across Γ (in the sense of traces in H 1/2 (Γ)), that the tangential components e(u Ω ) τ τ of the strain tensor, and the normal components σ(u Ω ) τ n and σ(u Ω ) nn of the stress tensor σ(u Ω ) = A Ω e(u Ω ) are continuous across Γ in the sense of traces in H -1/2 (Γ).

The (exact) shape optimization problem at stake is still of the form: (6.3) min

Ω∈U ad J(Ω), where J(Ω) = C(Ω) + Vol(Ω),
where C(Ω) is the compliance of the total structure:

(6.4) C(Ω) = D A Ω e(u Ω ) : e(u Ω ) dx = Γ N g • u Ω ds.
Here, we have denoted by : the usual Frobenius inner product over matrices. Again, more general objective functions could be considered with similar conclusions. The following result about the differentiation of J(Ω) was proved in [START_REF] Allaire | Damage evolution in brittle materials by shape and topological sensitivity analysis[END_REF]: Theorem 6.1. The objective function J(Ω) defined in (6.3) is shape differentiable at any admissible shape Ω ∈ U ad , and its shape derivative reads, for an arbitrary deformation θ ∈ Θ ad ,

J (Ω)(θ) = Γ g S Ω θ • n ds, where g S Ω = -σ(u Ω ) nn : [e(u Ω ) nn ] -2σ(u Ω ) nτ : [e(u Ω ) nτ ] + [σ(u Ω ) τ τ
] : e(u Ω ) τ τ + , where we recall that the notation [α] stands for the jump α 0α 1 of a discontinuous quantity α across Γ, taking values α 0 , α 1 inside Ω 0 and Ω 1 respectively. Remark 6.1. For further use, it is useful to give another expression to the above shape derivative. Using the transmission conditions associated to (6.2), a calculation reveals that the discontinuous components of the strain and stress tensors may be expressed in terms the continuous ones as:

[e(u Ω ) τ n ] = 1 2µ σ(u Ω ) τ n , [e(u Ω ) nn ] = 1 2µ + λ σ(u Ω ) nn - λ 2µ + λ tr(e(u Ω ) τ τ ), [σ(u Ω ) τ τ ] = [2µ] e(u Ω ) τ τ + 2µλ 2µ + λ tr(e(u Ω ) τ τ ) + λ 2µ + λ σ(u Ω ) nn δ i,j .
Using these relations, it is a simple matter to prove that:

(6.5) g S Ω = [2µ]e(u Ω ) τ τ : e(u Ω ) τ τ + 2µλ 2µ + λ tr(e(u Ω ) τ τ ) 2 - 1 2µ + λ σ(u Ω ) 2 nn + 2λ 2µ + λ σ(u Ω ) nn tr(e(u Ω ) τ τ ) - 1 µ σ(u Ω ) τ n • σ(u Ω ) τ n + .
The calculation of the topological derivative of J(Ω) is performed in [START_REF] Amstutz | Sensitivity analysis with respect to a local perturbation of the material property[END_REF] in the particular case where the Hooke's tensors A 0 , A 1 are proportional. In our precise case of interest (and in the more general case where these tensors are anisotropic), we refer to [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Giusti | Topological sensitivity analysis in heterogeneous anisotropic elasticity problem. Theoretical and computational aspects[END_REF], from which we quote the result in the case of two space dimensions. Theorem 6.2. For any two-dimensional admissible shape Ω ∈ U ad , the function J(Ω) defined in (6.3) has a topological derivative g T Ω (x) at any point x ∈ Ω 0 ∪ Ω 1 . Its expression reads: g T Ω (x) = s Ω (x) Pσ(u Ω (x)) : e(u Ω (x)), where s Ω (x) is defined by (2.8) and P is the fourth-order Pólya-Szegö polarization tensor, given by:

(6.6) ∀e ∈ S(R d ), Pe = 1 ρ 2 ρ 3 + τ 1 (1 + ρ 2 )(τ 1 -ρ 3 )e + 1 2 (ρ 1 -ρ 2 ) ρ 3 (ρ 3 -2τ 3 ) + τ 1 τ 2 ρ 1 ρ 3 + τ 2 tr(e)I ,
where (6.7)

ρ 1 = 1 + ν 1 -ν , ρ 2 = 3 -ν 1 + ν , ρ 3 = E E , τ 1 = 1 + ν 1 + ν , τ 2 = 1 -ν 1 -ν and τ 3 = ν (3ν -4) + 1 ν(3ν -4) + 1 .
and we have posed:

• If x ∈ Ω 1 , E = E 1 , E = E 0 , ν = ν 1 , and ν = ν 0 , • If x ∈ Ω 0 , E = E 0 , E = E 1 , ν = ν 0 , and ν = ν 1 .
6.2. The smoothed approximation in the linear elasticity setting.

Let us now describe the relaxed density optimization problem associated to (6.3). For an arbitrary density function h ∈ L 2 (D), we define a smoothed Hooke's tensor A h,ε via the formula: (6.8) ∀e ∈ S(R d ), A h,ε e = 2µ h,ε e + λ h,ε tr(e)I, cooked from the approximate Lamé coefficients: (6.9)

λ h,ε = λ 0 + (λ 1 -λ 0 )α(L ε h), and µ h,ε = µ 0 + (µ 1 -µ 0 )β(L ε h).
In the above formula, L ε : L 2 (D) → L ∞ (D) is either of the regularizing operators presented in Sections 4.1.1 and 4.1.2, and α, β ∈ C 2 (R) are two functions satisfying (3.3). In this context, the approximate solution u h,ε to the linear elasticity system is the unique solution in H 1 Γ D (D) d to:

(6.10)

   -div(A h,ε e(u h,ε )) = 0 in D, A h,ε e(u h,ε )n = g on Γ N , u h,ε = 0 on Γ D .
Our relaxed density optimization problem is still of the form:

(6.11) min h∈K J ε (h), where J ε (h) = C ε (h) + Vol ε (h),
and C ε (h) reads: (6.12)

C ε (h) = D A h,ε e(u h,ε ) : e(u h,ε ) dx.
The differentiation of J ε (h) is investigated in the same way as in Proposition 3.1, and we omit the proof:

Proposition 6.3. The relaxed functional J ε (h) is Fréchet differentiable on L 2 (D) and:

∀ h ∈ L 2 (Ω), J ε (h)( h) = D g ε h dx,
where:

g h,ε = -2(µ 1 -µ 0 )L ε ((β (L ε h)e(u h,ε ) : e(u h,ε ))) -(λ 1 -λ 0 )L ε α (L ε h)div(u h,ε ) 2 + L ε 1.
Last but not least, in the particular case where h = χ Ω is the characteristic function of a smooth domain Ω D, we use the notations

A Ω,ε , λ Ω,ε , µ Ω,ε , u Ω,ε , J ε (Ω) and g Ω,ε instead of A h,ε , λ h,ε , µ h,ε , u h,ε , J ε (h)
and g h,ε respectively. The density optimization problem (6.11) gives rise to a smoothed shape optimization problem of the form (3.9) in the line of Section 3.1: (6.13) min Ω∈U ad J ε (Ω).

Convergence of the approximate problem to its exact equivalent.

As in Sections 4 and 5, we examine the convergence of the approximate shape optimization problem (6.13) to the exact one (6.3).

In this context, one can prove along the line of Section 4.1 that the equivalent of Proposition 4.1 holds regarding the convergence of the approximate Lamé parameters λ Ω,ε and µ Ω,ε to their exact counterparts λ Ω and µ Ω for a given admissible shape Ω ∈ U ad .

When it comes to the asymptotic behavior of the spatial derivatives of the elastic displacement u Ω,ε , the following proposition is the exact counterpart of the convergence results of Theorem 4.5 in the context of the linear elasticity equations. We state it precisely for the sake of convenience, but the proof is omitted. Proposition 6.4. Let Ω ∈ U ad be an admissible shape. Then, (i) Let U be an open set such that U Ω 0 or U Ω 1 . Then for any n ≥ 1, one has:

u Ω,ε → u Ω strongly in H n (U ) d .
(ii) Let V be a neighborhood of Γ on which the projection mapping p Γ given by (4.1) is well-defined. Then the following convergences hold:

e τ τ (u Ω,ε ) ε→0 -→ e τ τ (u Ω ) strongly in H 1 (V ) (d-1) 2 , σ τ n (u Ω,ε ) ε→0 -→ σ τ n (u Ω ) strongly in H 1 (V ) d-1 , σ nn (u Ω,ε ) ε→0 -→ σ nn (u Ω ) strongly in H 1 (V ).
The main result of this section is now the following: Theorem 6.5. Let Ω ∈ U ad be an admissible shape. Then, (i) The approximate objective function J ε (Ω) converges to the exact one J(Ω) as ε → 0.

(ii) For any 1 ≤ p < ∞, the following 'interior' convergence of the approximate gradient holds:

g Ω,ε ε→0 -→ -2(µ 1 -µ 0 )β (0)e(u Ω ) : e(u Ω ) -(λ 1 -λ 0 )α (0)(div(u Ω )) 2 + strongly in L p (Ω 0 ), -2(µ 1 -µ 0 )β (1)e(u Ω ) : e(u Ω ) -(λ 1 -λ 0 )α (1)(div(u Ω )) 2 +
strongly in L p (Ω 1 ).

(iii) The following 'boundary' convergence of the approximate gradient holds:

g Ω,ε ε→0 -→ g S Ω in L 1 (Γ).
Proof. The proofs of (i) and (ii) are analogous to those of Propositions 5.1 and 5.2, and we focus on that of (iii). To this end, we rely on the shorthands:

µ ε ≡ µ Ω,ε , λ ε ≡ λ Ω,ε , e ε ≡ e(u Ω,ε
), and σ ε ≡ A Ω,ε e(u Ω,ε ). Proposition 6.4 reveals which are the components of the tensors e ε and σ ε that behave well in the limit ε → 0, and we express all their remaining components in terms of these ones; in particular:

e ε τ n = 1 2µ ε σ ε τ n and e ε nn = 1 2µ ε + λ ε (σ ε nn -λ ε tr(e ε τ τ )).
Now, the expression of the approximate gradient g Ω,ε reads:

g Ω,ε = -2(µ 1 -µ 0 )L ε (β (L ε χ Ω )e ε : e ε ) -(λ 1 -λ 0 )L ε α (L ε χ Ω )tr(e ε ) 2 + L ε 1,
and we expand it as:

g Ω,ε = -2(µ 1 -µ 0 )L ε (β (L ε χ Ω )e ε τ τ : e ε τ τ ) -L ε (µ1-µ0)β (LεχΩ) µ 2 ε σ ε τ n • σ ε τ n -L ε 2(µ1-µ0)β (LεχΩ) (2µε+λε) 2 (σ ε nn -λ ε tr(e ε τ τ )) 2 -(λ 1 -λ 0 )L ε α (L ε χ Ω )( 2µε 2µε+λε tr(e ε τ τ ) + 1 2µε+λε σ ε nn ) 2 + L ε 1.
which we may alternatively rewrite as, reordering terms:

g Ω,ε = -2(µ 1 -µ 0 )L ε (β (L ε χ Ω )e ε τ τ : e ε τ τ ) -L ε (µ1-µ0)β (LεχΩ) µ 2 ε σ ε τ n • σ ε τ n -L ε 2(µ1-µ0)β (LεχΩ)+(λ1-λ0)α (LεχΩ) (2µε+λε) 2 (σ ε nn ) 2 -L ε 2λ 2 ε (µ1-µ0)β (LεχΩ)+4µ 2 ε (λ1-λ0)α (LεχΩ) (2µε+λε) 2 tr(e ε τ τ ) 2 +L ε 4λε(µ1-µ0)β (LεχΩ)-4µε(λ1-λ0)α (LεχΩ) (2µε+λε) 2
tr(e ε τ τ )σ ε nn + L ε 1. Like in the proof of Theorem 5.3, one may prove that

g Ω,ε -g Ω,ε ε→0 -→ 0 in L 1 (Γ),
where g Ω,ε is defined by:

g Ω,ε = -2(µ 1 -µ 0 )L ε (β (L ε χ Ω )) e ε τ τ : e ε τ τ -L ε (µ1-µ0)β (LεχΩ) µ 2 ε σ ε τ n • σ ε τ n -L ε 2(µ1-µ0)β (LεχΩ)+(λ1-λ0)α (LεχΩ) (2µε+λε) 2 (σ ε nn ) 2 -L ε 2λ 2 ε (µ1-µ0)β (LεχΩ)+4µ 2 ε (λ1-λ0)α (LεχΩ) (2µε+λε) 2 tr(e ε τ τ ) 2 +L ε 4λε(µ1-µ0)β (LεχΩ)-4µε(λ1-λ0)α (LεχΩ) (2µε+λε) 2
tr(e ε τ τ )σ ε nn + . We now apply Lemma 5.4 repeatedly to each term in the above expression, using the shorthands:

λ(t) = λ 0 + (λ 1 -λ 0 )α(t) and µ(t) = µ 0 + (µ 1 -µ 0 )β(t).
• Using ϕ(t) = µ(t), we obtain:

lim ε→0 ((µ 1 -µ 0 ) L ε β (L ε χ Ω )) = µ 1 -µ 0 , • Using ϕ(t) = 1 µ(t) , ϕ (t) = -(µ1-µ0)β (t) µ(t) 2
we obtain:

lim ε→0 L ε (µ 1 -µ 0 )β (L ε χ Ω ) µ 2 ε = 1 µ 0 - 1 µ 1 , • Using ϕ(t) = 1 2µ(t)+λ(t) , ϕ (t) = -2µ (t)+λ (t) (2µ(t)+λ(t)) 2 , we obtain: lim ε→0 L ε 2(µ 1 -µ 0 )β (L ε χ Ω ) + (λ 1 -λ 0 )α (L ε χ Ω ) (2µ ε + λ ε ) 2 = 1 2µ 0 + λ 0 - 1 2µ 1 + λ 1 , • Using ϕ(t) = 2µ(t)λ(t) 2µ(t)+λ(t) , ϕ (t) = 4µ 2 (t)λ (t)+2λ 2 (t)µ (t) (2µ(t)+λ 2 (t)) 2 yields: lim ε→0 L ε 2λ 2 ε (µ 1 -µ 0 )β (L ε χ Ω ) + 4µ 2 ε (λ 1 -λ 0 )α (L ε χ Ω ) (2µ ε + λ ε ) 2 = 2µ 1 λ 1 2µ 1 + λ 1 - 2µ 0 λ 0 2µ 0 + λ 0 , • Using ϕ(t) = 2λ(t) 2µ(t)+λ(t) , ϕ (t) = 4µ(t)λ (t)-4µ (t)λ(t) (2µ(t)+λ(t)) 2 yields: lim ε→0 L ε 4λ ε (µ 1 -µ 0 )β (L ε χ Ω ) -4µ ε (λ 1 -λ 0 )α (L ε χ Ω ) (2µ ε + λ ε ) 2 = 2λ 0 2µ 0 + λ 0 - 2λ 1 2µ 1 + λ 1 .
In view of the expression (6.5), this ends the proof of the theorem.

Connection with the shape and topological derivatives of the exact problem

At this point, the results obtained about the asymptotic behavior of the smoothed shape optimization problem (3.9) may be summarized as follows, in the language of the two-phase conductivity equation for simplicity: for a fixed admissible shape Ω ∈ U ad ,

• The smoothed conductivity γ Ω,ε defined in (3.1)

converges to γ Ω strongly in L p (D) (1 ≤ p < ∞) as ε → 0 (Theorem 4.1
). • The smoothed potential u Ω,ε defined in (3.4) converges to u Ω in the sense that Theorem 4.5 holds.

• The smoothed value J ε (Ω) of the objective function converges to its exact counterpart J(Ω).

• The restriction to Γ of the approximate gradient g Ω,ε defined by (3.7) converges to the shape gradient g S Ω (x) of the exact objective J(Ω), defined by (2.7) (see Theorem 5.3):

g Ω,ε ε→0 -→ g S Ω , in L p (Γ), 1 ≤ p < ∞.
• The restriction to Ω 0 ∪ Ω 1 of g Ω,ε behaves as follows (see Proposition 5.2):

g Ω,ε ε→0 ---→ -(γ 1 -γ 0 )ζ (0)|∇u Ω | 2 + strongly in L p (Ω 0 ), -(γ 1 -γ 0 )ζ (1)|∇u Ω | 2 + strongly in L p (Ω 1 ),
and we shall soon connect these expressions to that of Theorem 2.2 for the topological gradient g T Ω of J(Ω).

As we are about to see, Proposition 5.2 and Theorem 6.5 make it possible to identify the limit of the approximate gradient g Ω,ε outside the interface Γ with the topological gradient g T Ω of the original cost J(Ω). This can be achieved through appropriate choices of the interpolation functions ζ, α and β, as we now discuss.

The case of the conductivity equation.

In this context, using Theorem 2.2 and Proposition 5.2, one easily observes that g Ω,ε → g T Ω in L p (Ω 0 ∪ Ω 1 ) as ε → 0, provided the interpolation profile ζ satisfies:

(7.1) ζ(0) = 0, ζ(1) = 1, ζ (0) = k 0 and ζ (1) = k 1 .
Multiple possibilities are therefore available as for the function ζ:

• The simplest choice is a third-order polynomial function ζ p (t) of the form:

ζ p (t) = a 3 t 3 + a 2 t 2 + a 1 t + a 0 ,
where a 0 , a 1 , a 2 , a 3 are uniquely determined by the relations (7.1).

• An alternative possibility is inspired by the homogenization theory: a classical result (see [START_REF] Allaire | Conception optimale de structures[END_REF] or [START_REF] Allaire | Shape optimization by the homogenization method[END_REF], Th. 2.2.3.1) states that a d × d tensor with eigenvalues λ 1 , ..., λ d , can be realized as a homogenized conductivity tensor obtained by mixing the two materials with conductivities γ 0 and γ 1 in proportions t and 1t (for t ∈ (0, 1)) if and only if it fulfills the so-called Hashin-Shtrikman bounds. Assuming for a moment that γ 0 < γ 1 to ease notations, these bounds read:

(7.2) λ - t ≤ λ i ≤ λ + t , where λ - t = t γ 0 + 1 -t γ 1 -1 , λ + t = tγ 0 + (1 -t)γ 1 , and (7.3)
d i=1 1 λ i -γ 0 ≤ 1 λ - t -γ 0 + d i=1 d -1 λ + t -γ 0
, and

d i=1 1 γ 1 -λ i ≤ 1 γ 1 -λ - t + d i=1 d -1 γ 1 -λ + t .
From (7.2) and ( 7.3), we easily infer a necessary and sufficient condition for a scalar conductivity (corresponding to an isotropic conductivity tensor λ 1 = ... = λ d ) to be achievable as a mixture of the two materials with conductivities γ 0 and γ 0 . Since (7.2) and ( 7.3) behave as rational functions of the volume fraction t, a natural idea consists in defining ζ(t) = ζ r (t) under the form:

ζ r (t) = a 2 t 2 + a 1 t + a 0 bt + 1 ,
where the coefficients a 0 , a 1 , a 2 and b are again uniquely determined by the relations (7.1). Interestingly, in two space dimensions, in the limit γ 0 → 0, a simple calculation reveals that the interpolating functions ζ p (t) and ζ r (t) produced by the above procedures amount to the same second-order polynomial ζ p (t) = ζ r (t) = t 2 . In Figure 3, the common values of the polynomial and rational interpolation functions ζ p (t) and ζ r (t) are represented in two space dimensions, in the (formal limit) case γ 0 = 0 and γ 1 = 1.

The case of the linear elasticity system.

So as to directly observe and check the positive definiteness of the interpolated Hooke's tensor (6.8) (6.9), we rely on the description of this tensor in terms of the shear and bulk moduli µ and κ respectively.

Assuming an interpolation of the form (6.9) between the material properties of both phases Ω 0 , Ω 1 , the bulk modulus κ h,ε of the approximate Hooke's tensor A h,ε reads:

κ h,ε = κ 0 + (κ 1 -κ 0 ) α(L ε h), where α(t) := 2 d β(t) + α(t).
Then, Theorem 6.2 imposes explicit conditions on the values and derivatives of α(t) and β(t) (thus of α(t)) at t = 0 and 1 so that the convergence g Ω,ε → g T Ω holds in L p (Ω 0 ∪ Ω 1 ), for 1 ≤ p < ∞. Like in the conductivity case, polynomial functions α p (t), β p (t) and rational functions α r (t), β r (t) are possible when it comes to constructing explicitly such interpolation functions.

Recall from [START_REF] Amstutz | Connections between topological sensitivity analysis and material interpolation schemes in topology optimization[END_REF] that in two space dimensions, and in the particular case where ν 0 = ν 1 = 1/3, elementary albeit tedious calculations reveal that both polynomial interpolation profiles α p (t) for λ and β p (t) for µ coincide. Hence, (6.8) reduces to:

A h,ε = A 0 + (A 1 -A 0 )β p (L ε h);
besides, as E 0 → 0, β p (t) ≈ t 3 , which is incidently the most commonly used power law for interpolationg Hooke's tensors in the context of the SIMP method.

An illustration of both interpolation procedures is presented In Figure 4. Notice that the interpolation functions shown here satisfy the Hashin-Shtrikman bounds. Depending on the elastic coefficients of the interpolated materials, this is not always the case using polynomial functions, as already remarked in the context of the SIMP method [START_REF] Bendsøe | Topology Optimization. Theory, Methods, and Applications[END_REF], in contrast to the suggestion of our numerical investigations with rational functions.

E 0 = 0.001, E 1 = 1, ν 0 = ν 1 = 1/3.
In both figures, the Hashin-Shtrikman bounds are depicted in blue, and the polynomial (resp. rational) interpolation functions correspond to the yellow (resp. brown) curves. The red segments illustrate the conditions imposed on their derivatives at t = 0 and t = 1.

Description of the numerical methods

In this section, we describe the two main ingredients in the device of numerical algorithms from the theoretical investigations in the previous sections. We first examine several possibilities offered by the numerical framework as regards the regularizing operator L ε : L 2 (D) → L ∞ (D). For simplicity, our discussion takes place in the two-dimensional setting, but its extension to three space dimensions is straightforward.

Throughout this section, the domain D is equipped with a conforming, triangular mesh T composed of K triangles T k , k = 1, ..., K, and J vertices p j , j = 1, ..., J. Let us recall the classical Finite Element spaces:

• V 0 ⊂ L 2 (D)
is the finite-dimensional space of Lagrange P 0 Finite Element functions on T , i.e. of constant functions in restriction to each triangle T k ∈ T . A basis of V 0 is composed of the functions N 0 k , k = 1, ..., K, where

N 0 k ≡ 1 on T k and N k ≡ 0 on T k , k = k . • V 1 ⊂ L 2 (D)
is the finite-dimensional space of Lagrange P 1 Finite Element functions, i.e. of affine functions in restriction to each triangle T k ∈ T . A basis of V 1 is composed of the functions N 1 j , j = 1, ..., J, where N 1 j is the unique element in V 1 such that N 1 j (p j ) = 1 if j = j and 0 otherwise. We now propose three discrete versions of the operator L ε , mapping an arbitrary density h ∈ L 2 (D) to a piecewise constant function L ε h ∈ V 0 . The first two ones are quite heuristic; they are inspired by the convolution operator L conv ε described in Section 3.2.1 and take advantage of the Finite Element framework. The third one is a direct discretization of the operator L ell ε of Section 3.2.2.

8.1.1. The ' P 0 kernel' operator P 0 .

Our first operator P 0 : L 2 (D) → V 0 is the L 2 (D) orthogonal projection from L 2 (D) to the subspace V 0 :

∀h ∈ L 2 (D), P 0 h = k 1 |T k | (N 0 k , h) L 2 (D) N 0 k ,
i.e. for each triangle T k ∈ T , P 0 h| T k is the average of h over T k . P 0 is easily seen to be a self-adjoint operator from L 2 (D) into itself. In this case, the mesh size -i.e. the maximum length of an edge in the mesh -plays the role of the smoothing parameter ε.

When P 0 is used as the regularizing operator L ε , the gradient g h,ε of the optimized functional J ε (h), calculated via (3.7), naturally belongs to the space V 0 , and it is therefore constant in restriction to each triangle of the mesh. This may cause difficulties when the numerical algorithm used to solve the optimization problem (3.5) or (3.9) requires that this information be supplied at the vertices of T (this is the case of the level set algorithm presented in Section 8.2.2 below). Then, a gradient defined at nodes can be recovered by using e.g. a Clément Finite Element interpolate; see e.g. [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. 8.1.2. The ' P 1 kernel' operator P 1 .

Let us now define the operator P 1 : L 2 (D) → V 0 by:

∀h ∈ L 2 (D), P 1 h = K k=1 (Q k , h) L 2 (D) N 0 k ,
where for any triangle T k , with vertices p j1 , p j2 , p j3 , the ' P 1 kernel' Q k ∈ V 1 is defined by (see Figure 5):

(8.1) Q k = 1 3 i=1 (N 1 ji , 1) L 2 (D) 3 i=1 N 1 ji .
This operator is not self-adjoint from L 2 (D) into itself; nevertheless, as observed in Remark 3.1, the considerations of Section 3.1 are easily adapted to this case. The gradient g h,ε of J ε (h), which is then calculated as (3.8), naturally belongs to V 1 . Therefore, P 1 smoothes the density h in such a way that the values of P 1 h inside one triangle T k of the mesh not only involve the values of h on T k but also those on the neighboring triangles. Our last operator L ell T : L 2 (D) → V 0 is defined as follows: for h ∈ L 2 (D), L ell T h = P 0 q T ,h , where q T ,h ∈ V 1 is the numerical solution of (3.12) given by the P 1 Lagrange Finite Element method on the mesh T ; more precisely:

q T ,h = J i=1   J j=1 R ij (N 1 j , h) L 2 (D)   N 1 i , where R = (ε 2 K + M ) -1 ,
and K, M stand respectively for the J × J stiffness and mass matrices associated to the resolution of (3.12) by means of the Lagrange P 1 Finite Element method:

K ij = D ∇N 1 j • ∇N 1 i dx, and M ij = D N 1 j N 1 i dx.
L ell T is not self-adjoint as an operator from L 2 (D) into itself, but the mapping h → q T ,h is. Again, the corresponding discrete gradient g h,ε given by (3.7) is naturally a V 1 function.

The smoothing effect induced by L ell T is controlled at the same time by the parameter ε and the mesh size. It is natural to demand that the mapping h → q T ,h satisfy the discrete maximum principle, so that regularized densities L ell T h take values between 0 and 1. This requirement prevents ε from being 'very small'; our experience suggests a choice of ε of the order of the mesh size.

On a different note, notice that, contrary to the previous operators P 0 and P 1 , the dependence of L ell T h on h is non local, i.e. even if h has small support, L ell T h is non trivial on the whole computational domain D (which reflects the strong maximum principle satisfied by the continuous equation (3.12)). Therefore, densites or characteristic functions smoothed using L ell T are expected to show large grayscale areas.

8.2. Two numerical optimization algorithms.

We now present two numerical algorithms inspired by the previous considerations. The first one is a simple projected gradient algorithm devoted to the resolution of the density optimization problem (3.5), and the second one relies on the level set method for solving the approximate shape optimization problem (3.9).

These algorithms are motivated by the conclusions of Proposition 5.2 and Theorem 5.3, whereby performing variations of a shape Ω with respect to the approximate gradient g Ω,ε defined in (3.7) should closely approximate at the same time a shape optimization algorithm based on the shape derivative, and a topology optimization algorithm based on the topological gradient, provided the interpolation profiles ζ, α and β are chosen according to the discussion in Section 7.

For simplicity, the discussion unfolds in the case of the two-phase conductivity equation of Section 2.

8.2.1.

A projected gradient algorithm for solving the density optimization problem.

Our first numerical algorithm is a simple projected gradient algorithm for the resolution of the density optimization problem (3.5), which is very close in essence to the standard filtered SIMP algorithm. Starting from an initial density h 0 ∈ L ∞ (D, [0, 1]), the algorithm produces a sequence h n ∈ L ∞ (D, [0, 1]), n = 0, ..., passing from h n to h n+1 by following the steepest descent direction g hn,ε supplied by (3.7), then thresholding in such a way that h n+1 ∈ L ∞ (D, [0, 1]) (hence the terminology 'projected gradient'), i.e.: h n+1 (x) = max(0, min(1, h n (x)τ n g hn,ε (x))), for a sufficiently small pseudo-time step τ n . The process ends when the relative difference ||h n+1 -h n || L 2 (D) /||h n || L 2 (D) between two consecutive values of the optimized density is smaller than a fixed, user-defined threshold r pr f . Notice that the resulting optimized design from this algorithm is a density function, which may a priori contain regions with intermediate densities, even if they are penalized by the interpolation profile ζ. A popular method to eliminate possible such regions appearing in the final design consists in projecting the latter onto the set of functions taking only values 0 and 1 as a post-processing in order to 'interprete' it as a true 'black and white' shape. 8.2.2. An algorithm using the level set method.

Our second algorithm is devoted to the numerical resolution of the approximate shape optimization problem (3.9); hence, it presents the important feature that it always involves a real, 'black and white' shape (and not a density function). This algorithm is similar to the level set based topology optimization algorithm introduced in [START_REF] Amstutz | A new algorithm for topology optimization using a level-set method[END_REF], then analyzed in [START_REF] Amstutz | Analysis of a level set method for topology optimization[END_REF].

According to the general philosophy of the level set method (see [START_REF] Osher | Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[END_REF] for the seminal article and [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF] for an overview), every considered shape Ω is viewed as the negative subdomain of an associated 'level set' function

φ : D → R, i.e.:    φ(x) < 0 if x ∈ Ω, φ(x) = 0 if x ∈ Γ, φ(x) > 0 if x ∈ D \ Ω,
Relying on the necessary optimality condition for Problem (3.9) expressed in Proposition 3.3, if a shape Ω ∈ U ad is optimal with respect to J ε (Ω) (in the sense of Definition 2.3), then g Ω,ε is an associated level set function for Ω. Noting in addition that if ψ is one level set function for Ω, then so is cψ for any constant c > 0, we are led to search for one function ψ such that:

(8.2) ||ψ|| L 2 (D) = 1, and ψ = 1 ||g Ω,ε || L 2 (D)
g Ω,ε , with Ω = {ψ < 0} .

Applying a fixed-point algorithm with relaxation to (8.2), elementary calculations produce the sequences of level set functions ψ n and associated shapes Ω n := {x ∈ D, ψ n (x) < 0} defined by: (8.3)

ψ n+1 = 1 sin a n (sin((1 -τ n )a n )ψ n + sin(τ n a n ) g n ) ,
where

g n = 1 ||gn|| L 2 (D)
g n is the normalized version of the gradient g n := g Ωn,ε , a n ∈ [0, π] is the angle a n = arccos((ψ n , g n ) L 2 (D) ), and τ n plays the role of a time step.

The process described by (8.3) is iterated until the relative difference ||χ Ωn+1,εχ Ωn,ε || L 2 (D) /||χ Ωn,ε || L 2 (D) becomes smaller than a user-defined threshold r ls f . Remark 8.1. There are traditionnally two ways for evaluating numerically the gradient information g h,ε (or g Ω,ε ) see e.g. [START_REF] Quarteroni | Numerical Models for Differential Problems[END_REF], §16.13. The 'optimize-then-discretize' strategy simply consists in discretizing the continuous expression (3.7) for g h,ε , e.g. by using Finite Elements. In the applications of this article, we rely on the competing 'discretize-then-optimize' paradigm: density functions h are first discretized -in our case, they are decomposed on a basis of the Finite Element space V 0 (see Section 8) -, then g h,ε is calculated as the (finite-dimensional) gradient of the resulting finite-dimensional discretization of J ε (h).

Numerical illustrations

In this section, we present several numerical illustrations of the previous material, limiting ourselves to the physically more relevant context of the linearized elasticity equations of Section 6. In all the considered examples, the computational domain D is a polygon (see Remark 5.4 about the validity of the theoretical results of this paper in this context). It is equipped with a triangular mesh, and the elastic coefficients of the materials in Ω 0 and Ω 1 are E 0 = 10 -3 , E 1 = 1, ν 0 = ν 1 = 1/3. We systematically rely on the rational interpolation scheme for the interpolation profiles ζ, α and β (see Section 7).

At first, Section 9.1 takes place in a model situation, and we appraise the main conclusions of Theorem 6.5 about the convergence of the approximate shape optimization problem (6.13) to its exact counterpart (6.3) as the smoothing parameter ε vanishes.

Then, in Sections 9.2 and 9.3, we apply the constructions and algorithms discussed in Sections 7 and 8 to two classical examples in shape and topology optimization. In both cases, the initial design is the whole computational domain, i.e. h 0 ≡ 1 (resp. ψ 0 ≡ -1/ |D|) on D when the projected gradient algorithm of Section 8.2.1 (resp. the level set algorithm of Section 8.2.2) is used. The parameter ε is chosen of the order of the mesh size, and as far as the stopping criterion is concerned, the common value 10 -3 is used for the thresholds r pr f and r ls f . 9.1. Numerical assessment of the convergence of the approximate shape optimization problem to its exact counterpart.

A rectangular domain D of size 2 × 1 is divided into two regions Ω 1 = Ω and Ω 0 = D \ Ω 1 corresponding to its lower and upper halves respectively; D is clamped at its left-and right-hand sides, and a unit vertical load g = (0, -1) is applied on a region Γ N located at the middle of its lower side; see Figure 6. The domain D is discretized by a regular triangular mesh T , obtained from a Cartesian grid of size h by splitting each voxel into four triangles. In this setting, the interface Γ between Ω 0 and Ω 1 is explicitely represented in T .

We first calculate a close approximation to the solution u Ω to (6.2) by using a very fine mesh of D with size h = 1 320 (containing 819, 200 triangles). We thence infer precise calculations of the compliance C(Ω), of the shape gradient g S Ω and of the topological derivative g T Ω of the exact shape optimization problem (6.3) in this context, according to Formula (6.4) and Theorems 6.1 and 6.2. These quantities serve as reference values in the experiments of this section; they are referred to as the 'exact' compliance, shape gradient and topological gradient, and, abusing notations, they are still denoted by C(Ω), g S Ω and g T Ω , respectively. We then calculate the approximate compliance C ε (Ω) given by (6.12) and the derivative g Ω,ε of the approximate shape optimization problem (see Proposition 6.3) by using the discrete version of the regularizing operator L ε based on either the 'P 1 kernel' operator of Section 8.1.2 or the L ell T operator of Section 8.1.3, on successively refined meshes, with sizes h = 1 20 , 1 40 , 1 80 and 1 160 (corresponding to 3, 200, 12, 800, 51, 000 and 204, 800 triangles respectively), and for the values ε = h, h 2 , or h 10 in the latter case. In each situation, we evaluate:

• The value of the approximate compliance C ε (Ω),

• The relative error

e T Ω,ε := ω |g T Ω -g Ω,ε | 2 dx ω |g T Ω | 2 dx
between g Ω,ε and g T Ω , on the subset ω :=

x = (x 1 , x 2 ) ∈ Ω, x 2 < 1 2 -1 20
of Ω; i.e. ω is obtained from Ω by removing one layer of elements of the coarser mesh.

• The relative error

e S Ω,ε := Γ |g S Ω -g Ω,ε | 2 dx Γ |g S Ω | 2 dx
between g Ω,ε and g S Ω . The results of these experiments are reported on Figure 7, and are in good agreement with the conclusions of Theorem 6.5: as the mesh size h (and thereby the regularization parameter) decreases, the approximate compliance C ε (Ω) tends to its exact counterpart, and the approximate derivative g Ω,ε converges to the topological derivative g T Ω of the exact shape optimization problem when it is restricted to the subset ω ⊂ Ω, and to the shape gradient g S Ω when it is restricted to Γ.

A self-adjoint objective function for the elasticity problem.

A rectangular cantilever beam D of size 2 × 1, equipped with a regular triangular mesh T made of 19600 elements, is clamped at its left-hand side, and a unit vertical load g = (0, -1) is applied at the middle of its right-hand side, as depicted on Figure 8. The objective function of interest is a weighted sum of the compliance of shapes and their volume:

J(Ω) = D A Ω e(u Ω ) : e(u Ω ) dx + Vol(Ω),
where A Ω is defined in (6.1), and the elastic displacement u Ω is the solution to (6.2). The Lagrange multiplier is set to = 1.

We apply the projected gradient algorithm of Section 8.2.1 for solving the density optimization problem (6.11) in this context, and the level set algorithm of Section 8.2.2 for solving the associated approximate shape optimization problem (6.13). For the sake of comparison, the three different instances of the regularization operators L ε described in Section 8.1 are used, corresponding to as many different ways to calculate the cost function J(Ω) and the associated gradient. The resulting optimal topologies are shown in Figures 9 and10; the corresponding convergence histories are reported in Figure 11.

A non self-adjoint objective function.

In order to appraise the robustness of the proposed algorithms, we eventually turn to the study of a non self-adjoint objective function; as pointed out in Remark 5.3, the main conclusions of this article remain valid in this case.

Let us consider the optimal design of a gripping mechanism, as shown in Figure 12; the computational domain D is discretized with a triangular mesh T containing 19600 triangles. The structure is fixed at two nodes on its symmetry axis; the region Γ N where loads are applied is decomposed as Γ N = Γ 1 N ∪ Γ 2 N and: • Vertical loads g 1 = (0, ±10) are applied on Γ 1 N , which represent the pressure exerted by the user, • Vertical loads g 2 = (0, ±1) are applied on Γ 2 N , which account for the reaction force applied by the object undergoing the action of the gripping mechanism. The considered objective function J(Ω) reads:

J(Ω) = Γ N k • u Ω ds + Vol(Ω),
where the vector k is defined by:

• k = (0, -1) on the upper side of Γ 1 N and k = (0, 1) on the lower side of Γ 1 N ,

Projected gradient

Level-set Then (see [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF], Prop. 9.3):

L ell T J ε (h) = 2.4 J ε (Ω) = 2.394 P 1 J ε (h) = 2.294 J ε (Ω) = 2.282 P 0 J ε (h) = 2.247 J ε (Ω) = 2.216
Proposition A.1. Let ϕ ∈ L p (D); the following assertions are equivalent:

(1) ϕ belongs to W 1,p (D).

(2) There exists a constant C > 0 such that: Let us pass to the proof of Theorem 4.5 properly speaking. Throughout this appendix, a shape Ω ∈ U ad is fixed; L ε stands for either one of the operators L conv ε or L ell ε constructed in Section 3.2, and we rely on the shorthands: γ ε ≡ γ Ω,ε , and u ε ≡ u Ω,ε . Also, C consistently denotes a positive constant, that may change from one line to the other, but is in any event independent of ε and the parameter h (to be introduced).

Proof of (i):

Without loss of generality, we assume that U Ω 0 , and we introduce two other subsets V, W of D such that U V W Ω 0 . Let χ be a smooth cutoff function such that χ ≡ 1 on U , and χ ≡ 0 on V . Our aim is to prove that: (A.6)

D γ ε |∇D h w ε | 2 dx = - D D h γ ε τ h (∇w ε ) • ∇D h w ε dx - D τ h (γ ε -γ Ω )∇D h v Ω • ∇D h w ε dx - D D h (γ ε -γ Ω )∇v Ω • ∇D h w ε dx + D (f ε -f Ω )D -h D h w ε dx + D D h (h ε -h Ω ) • ∇D h w ε dx.
To achieve the last identity, we have used the 'discrete integration by parts': Therefore, (A.2) holds for m = 2. The case m > 2 is obtained by iterating the previous argument.

Proof of (ii):

Since Γ is compact, it is enough to prove that the estimates (4.6) hold in an open neighborhood U of any arbitrary point x 0 ∈ Γ; namely, we introduce two other subset V, W of D such that U V W D; let χ be a smooth cutoff function, which equals 1 on U et 0 on D \ V ; we aim to prove that: where, again v ε = χu ε and v Ω = χu Ω . Without loss of generality, we assume that Γ is flat in U , that is: Ω ∩ U = {x ∈ U, x d < 0} is a piece of the lower half-space. The general case is recovered from this one by a standard argument of local charts for the smooth boundary Γ. To prove (A.8), it is therefore enough to consider the case τ = e i , i = 1, ..., d -1, the i th -vector of the canonical basis (e 1 , ..., e d ) of R d , which we now do.

Introducing f ε , f Ω ∈ L 2 (D), h ε , h Ω ∈ H 1 (D) d as in (A.5), and w ε = v εv Ω , we now follow the exact same trail as that leading to (A.6) in the proof of (i), using only vectors h of the form h = he i , with small enough h > 0:

(A.10) D γ ε |∇D h w ε | 2 dx = - D D h γ ε τ h (∇w ε ) • ∇D h w ε dx - D τ h (γ ε -γ Ω )∇D h v Ω • ∇D h w ε dx - D D h (γ ε -γ Ω )∇v Ω • ∇D h w ε dx + D (f ε -f Ω )D -h D h w ε dx + D D h (h ε -h Ω ) • ∇D h w ε dx
We now estimate each of the five terms in the right-hand side of (A.10):

• Using Hölder's inequality with 1 p + 1 q + 1 2 = 1, the first term is controlled as: is bounded since e i is a tangential direction to Γ.

D D h γ ε ∇w ε • ∇D h w ε dx ≤ C||∇v ε -∇v Ω || L q (W ) d ∂γ ε ∂x i L p (W ) ||∇D h w ε || L 2 (W ) d ,
• The second term is controlled as:

D τ h (γ ε -γ Ω )∇D h v Ω • ∇D h w ε dx ≤ C||γ ε -γ Ω || L p (D) ∂(∇v Ω ) ∂x i L q (W ) d ||∇D h w ε || L 2 (W ) d ,
where we have used the regularity theory for v Ω and the fact that e i is tangential to Γ, and where • One has:

||γ ε -γ Ω || L p (D)
D (f ε -f Ω )D -h D h w ε dx ≤ ||f ε -f Ω || L 2 (W ) ||D -h D h w ε || L 2 (W ) , ≤ ||f ε -f Ω || L 2 (W ) ||∇D h w ε || L 2 (W ) d ,
where we have used (A.1), and it is easily seen, using Proposition 4.4, that ||f ε -

f Ω || L 2 (W ) ε→0 -→ 0. • Likewise, D D h (h ε -h Ω ) • ∇D h w ε dx ≤ ∂h ε ∂x i - ∂h Ω ∂x i L 2 (W ) ||∇D h w ε || L 2 (W ) d
and using the fact that e i is a tangential direction to Γ, it comes that ∂hε ∂xi -∂hΩ ∂xi L 2 (W ) ε→0 -→ 0.

Putting things together, using Proposition A.1, we obtain that:

∂v ε ∂τ ε→0 -→ ∂v Ω ∂τ in H 1 (W ).
Remark that, using Meyer's Theorem 4.3 from the identity (A.10) together with the previous estimates for its right-hand side shows that (A.8) actually holds in W 1,p (W ) for some p > 2.

The only thing left to prove is then (A.9), where we recall the simplifying hypothesis n = e d . Actually, (A.8) combined with Proposition 4.1 and the above remark already prove that, for i = 1, ..., d -1,

∂ ∂x i γ Ω,ε ∂v ε ∂x d = γ Ω,ε ∂ 2 v ε ∂x i ∂x d = γ Ω,ε ∂ ∂x d ∂v ε ∂x i ε→0 -→ γ Ω ∂ 2 v Ω ∂x i ∂x d in L 2 (W ).
Hence, the only thing left to prove is that:

∂ ∂x d γ Ω,ε ∂v ε ∂x d ε→0 -→ ∂ ∂x d γ Ω,ε ∂v Ω ∂x d ;
this last convergence is obtained by using the original equations (2.2) and (3.4), and notably the facts that:

-div(γ Ω,ε ∇v ε ) = -div(γ Ω ∇v Ω ) = 0 on U, together with the previous convergences in (A.8) and Proposition 4.1:

∂ ∂x d γ Ω,ε ∂v ε ∂x d = - d-1 i=1 ∂ ∂x i γ Ω,ε ∂v ε ∂x i ε→0 -→ - d-1 i=1 ∂ ∂x i γ Ω ∂v Ω ∂x i = ∂ ∂x d γ Ω ∂v Ω ∂x d .
This completes the proof.

  where ζ : R → R satisfies ζ(0) = 0 and ζ(1) = 1.

Figure 1 .

 1 Figure 1. Setting of the two-phase conductivity problem. Denoting by χ Ω ∈ L ∞ (D, {0, 1}) the characteristic function of Ω, the conductivity γ Ω inside D reads: (2.1) γ Ω = γ 0 + (γ 1γ 0 )χ Ω . The domain D is insulated on a part Γ D ⊂ ∂D of its boundary, and an incoming flux g ∈ H -1/2 (Γ N ) is imposed on the complementary subset Γ N = ∂D \ Γ D . Neglecting sources, the voltage potential u Ω in D belongs to the Hilbert space H 1 Γ D (D) := u ∈ H 1 (D), u = 0 on Γ D , and is the unique solution to the conductivity equation:

  r)|+o(|B(x, r)|), where lim r→0 o(|B(x, r)|) |B(x, r)| = 0, and:

ε

  and L ell ε whose constructions are detailed in Section 3.2 below. • The 'interpolation profile' ζ ∈ C 2 (R) makes the connection between the intermediate, 'grayscale' values of the (smoothed) density L ε h and the material properties in the domain. It has the properties:

3. 2 . 1 ..

 21 Regularization via convolution with a smoothing kernel: definition of L conv ε The perhaps most intuitive construction of a bounded operator from L 2 (D) into L ∞ (D) relies on a function η ∈ C ∞ c (R d ) with the properties:

3 . 2 . 2 .

 322 where δ ∈ D (R d ) stands for the Dirac distribution. Owing to the properties of convolution, the operationh → L conv ε h := η ε * h,where h is extended by 0 outside D, accounts for a linear, bounded operator from L 2 (D) into L ∞ (D) -it actually maps L 1 (D) into C ∞ (D) -and it is clearly self-adjoint from L 2 (D) into itself. Regularization based on an elliptic partial differential equation: definition of L ell ε .

( 4 . 1 )

 41 p Γ : x → the unique y ∈ Γ such that |x -y|= d(x, Γ) := min y∈Γ |x -y|

4. 1 . 1 .

 11 The case where L ε = L conv ε .

Lemma 4 . 2 .

 42 The smoothed characteristic function χ Ω,ε belongs to D(R d ) and takes the values 1 on {x ∈ Ω, d(x, Γ) > ε} and 0 on x ∈ D \ Ω, d(x, Γ) > ε . Moreover, the following estimates hold in the tube Γ ε := {x ∈ D, d(x, Γ) < ε}:

Lemma 5 . 4 .

 54 Let ϕ : R → R be a function of class C 2 such that ϕ, ϕ and ϕ ∈ L ∞ (R). Let L ε : L 2 (D) → L ∞ (D)stand either for the operator L conv ε constructed in Section 3.2.1 or for that L ell ε of Section 3.2.2. Let the function ψ ε : Γ → R be defined as:

Figure 2 .

 2 Figure 2. Illustration of the second step in the proof of Lemma 5.4.

Figure 3 .

 3 Figure 3. Interpolation process in the conductivity setting, in the case γ 0 = 0 and γ 1 = 1. The harmonic and arithmetic means λ - t and λ + t defined in (7.2) are depicted in black. The red segments account for the conditions (7.1) on ζ. The Hashin-Shtrikman bounds for the common eigenvalue of an isotropic conductivity tensor (i.e. obtained by imposing λ 1 = λ 2 in (7.3)) are represented by the blue curves, and the common values of ζ p (t) and ζ r (t) in this case are depicted in yellow.

Figure 4 .

 4 Figure 4. Interpolation functions (left) β for the shear modulus µ and (right) α for the bulk modulus κ. The coefficients of the materials filling Ω 0 and Ω 1 read: E 0 = 0.001, E 1 = 1, ν 0 = ν 1 = 1/3. In both figures, the Hashin-Shtrikman bounds are depicted in blue, and the polynomial (resp. rational) interpolation functions correspond to the yellow (resp. brown) curves. The red segments illustrate the conditions imposed on their derivatives at t = 0 and t = 1.

8. 1 .

 1 Practical implementation of the regularizing operator L ε .

Figure 5 . 8 . 1 . 3 .

 5813 Figure 5. Graph of the (non normalized) function Q k defined in (8.1).

Figure 6 .

 6 Figure 6. Setting of the test case used in Section 9.1 for the evaluation of the convergence of the approximate shape optimization problem to its exact counterpart.

Figure 7 .Figure 8 .

 78 Figure 7. (Top, left) Convergence of the approximate compliance C ε (Ω) to its exact counterpart C(Ω) with respect to the mesh size; (top, right) relative error e TΩ,ε over the approximation of the topological derivative for various mesh sizes; (bottom) relative error e S Ω,ε over the approximation of the shape gradient for various mesh sizes.

Figure 10 .

 10 Figure 10. Optimal shapes obtained in the cantilever test-case of Section 9.2, (left column) using the projected gradient algorithm of Section 8.2.1 for the density optimization problem (3.5) and (right column) using the level set algorithm of Section 8.2.2 for the smoothed shape optimization problem (3.9). The exact density h (left) and characteristic function χ Ω (right) are represented.

Figure 11 .Figure 12 .

 1112 Figure 11. Convergence histories for the projected gradient and level-set algorithm with P 0 , P 1 and L ell ε regularization operators for the cantilever test case of Section 9.2.

( 3 )

 3 ≤ C||ψ|| L p (D) , for any v ∈ C ∞ c (D), i = 1, ..., d. There exists a constant C > 0 such that, for any open subset V D, lim sup h→0 ||D h ϕ|| L p (V ) ≤ C. Moreover, the smallest constant satisfying points (2) and (3) is C = ||∇ϕ|| L p (D) d .

-Figure 15 .

 15 Figure 15. Deformed configurations for the optimal gripping mechanisms in the setting of Section 9.3.

D

  w(D h z) dx = -D (D -h w)z dx,

Figure 16 .

 16 Figure 16. Convergence histories for the projected gradient and level-set algorithm with P 0 , P 1 and L ell ε regularization operators for the gripping mechanism test case of Section 9.3.

  in H 1 (W ), for any tangential vector field τ : Γ → R d , and(A.9) γ Ω,ε ∂v ε ∂n ε→0 -→ γ Ω ∂v Ω ∂n in H 1 (W ),

  where we have used Proposition A.1. Owing to Proposition 4.4, we know that ||∇v ε -∇v Ω || L q (D) d → 0 as ε → 0, while, from Proposition 4.1, ∂γε ∂xi L p (D)

ε→0-→ 0 (

 0 see again Proposition 4.1).• The third term in the right-hand side of (A.10) is estimated as:D D h (γ εγ Ω )∇v Ω • ∇D h w ε dx ≤ C||∇v Ω || L q (W ) ∂γ ε ∂x i -∂γ Ω ∂x i L p (W ) ||∇D h w ε || L 2 (W ) d , in which ∂γε ∂xi -∂γΩ ∂xi L 2 (W )ε→0-→ 0 because of Proposition 4.1.

  The approximate density optimization problem 3.1. Reformulation of the shape and topology optimization problem as a density optimization problem 3.2. Two examples of constructions of a regularizing operator L ε : L 2 (D) → L ∞ (D) 3.3. Optimality conditions for the approximate problem and orientation 4. Asymptotic behaviors of the approximate conductivity γ Ω,ε and potential u Ω,ε 4.1. Asymptotic behavior of the approximate conductivity γ Ω,ε 4.2. Asymptotic behavior of the approximate potential u Ω,ε 5. Convergence of the approximate objective function and its derivative 5.1. First convergence results 5.2. Convergence of the approximate gradient on the boundary Γ 6. Extension to the linear elasticity case 6.1. Description of the linear elasticity setting 6.2. The smoothed approximation in the linear elasticity setting 6.3. Convergence of the approximate problem to its exact equivalent 7. Connection with the shape and topological derivatives of the exact problem 7.1. The case of the conductivity equation 1. Introduction

	1. Introduction
	2. The exact two-phase shape optimization problem in the conductivity setting
	2.1. The two-phase conductivity setting
	2.2. Hadamard's boundary variation method and shape derivative
	2.3. Topological sensitivity analysis
	2.4. Local optimality conditions
	3.

Acknowledgements. C.D. was partially supported by the ANR OptiForm. A.F. has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n. 320815 (ERC Advanced Grant Project "Advanced tools for computational design of engineering materials" COMP-DES-MAT).

Projected gradient

Level-set

J ε (Ω) = 2.394 • k = (0, 2) on the upper side of Γ 2 N and k = (0, -2) on the lower side of Γ 2 N , so that it is expected that the elastic displacement of the resulting design shows a pinching of the jaws Γ 2 N without inducing an excessive displacement of the region Γ 1 N where the user applies forces. Again, the projected gradient algorithm and the level set algorithm are applied with the three different instances of the regularization operators L ε discussed in Section 8.1. The resulting optimal topologies are shown in Figures 13 and14, and the corresponding deformed configurations of the grip are displayed in Figure 15. The convergence histories are reported in Figure 16.

Appendix A. Sketch of the proof of Theorem 4.5

The proof is a little technical, but rests on classical arguments, and notably Nirenberg's technique of difference quotients [START_REF] Nirenberg | Remarks on strongly elliptic partial differential equations[END_REF], as exposed, e.g. in [START_REF] Brezis | Functional Analysis, Sobolev Spaces and Partial Differential Equations[END_REF] (see also the discussion in [START_REF] Dapogny | Shape optimization, level set methods on unstructured meshes and mesh evolution[END_REF], §4.8). We therefore limit ourselves to sketches of proofs.

Let us start with a useful result about the characterization of W 1,p spaces. Let 1 < p ≤ ∞, and 1 ≤ p < ∞ be such that 1 p + 1 p = 1. For a function ϕ ∈ L p (D), a subset V D, and a vector h ∈ R d with |h|< dist(V, ∂D), we define the difference quotient D h ϕ ∈ L p (V ) as:

x ∈ V, where τ h ϕ(x) := ϕ(x + h).

Projected gradient

Level set