Gabriel Delgado
email: gabriel.delgado@irt-systemx.fr

Marc Bonnet
email: mbonnet@ensta.fr

The topological derivative of stress-based cost functionals in anisotropic elasticity

The topological derivative of cost functionals J that depend on the stress (through the displacement gradient, assuming a linearly elastic material behavior) is considered in a quite general 3D setting where both the background and the inhomogeneity may have arbitrary anisotropic elastic properties. The topological derivative DJ(z) of J quantifies the asymptotic behavior of J induced by the nucleation in the background elastic medium of a small anisotropic inhomogeneity of characteristic radius a at a specified location z. The fact that the strain perturbation inside an elastic inhomogeneity remains finite for arbitrarily small a makes the small-inhomogeneity asymptotics of stress-based cost functionals quite different than that of the more usual displacement-based functionals.

The asymptotic perturbation of J is shown to be of order O(a 3 ) for a wide class of stress-based cost functionals having smooth densities. The topological derivative of J, i.e. the coefficient of the O(a 3 ) perturbation, is established, and computational procedures then discussed. The resulting small-inhomogeneity expansion of J is mathematically justified (i.e. its remainder is proved to be of order o(a 3 )). Several 2D and 3D numerical examples are presented, in particular demonstrating the proposed formulation of DJ on cases involving anisotropic elasticity and non-quadratic cost functionals.

Introduction

The topological derivative DJ(z) quantifies the perturbation induced to a cost functional J by the virtual creation of an object (in this work, an elastic inhomogeneity) occupying a region B a (z) with prescribed center z inside the solid and vanishingly small characteristic radius a. In structural optimization, the field DJ(z) helps directing the algorithm towards optimal topologies by indicating where creating new holes is most profitable, see e.g. [START_REF] Amstutz | A new algorithm for topology optimization using a level-set method[END_REF][START_REF] Samet | The topological asymptotic for the Helmholtz equation[END_REF] and also [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level-set method[END_REF] in conjunction with the shape derivative. Moreover, computational evidence [START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF][START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF][START_REF] Masmoudi | The topological asymptotic expansion for the Maxwell equations and some applications[END_REF] as well as more recent theoretical and experimental studies [START_REF] Ammari | Localization, stability, and resolution of topological derivative based imaging functionals in elasticity[END_REF][START_REF] Bellis | Acoustic inverse scattering using topological derivative of far-field measurements-based L 2 cost functionals[END_REF][START_REF] Tokmashev | Experimental validation of the topological sensitivity approach to elastic-wave imaging[END_REF] show that the topological derivative is also effective for flaw identification problems.

Most of the studies so far devoted to the topological derivative and its applications consider cost functionals that depend on the primary variable, namely the displacement field in the solid mechanics context of this work. Such formulations exploit the analysis of the asymptotic behavior as a → 0 of the perturbation induced to the displacement by the virtual creation of B a (z), for which an abundant literature is available [START_REF] Ammari | Polarization and moment tensors with applications to inverse problems and effective medium theory[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Cedio-Fengya | Identification of conductivity imperfections of small diameter by boundary measurements. Continuous dependence and computational reconstruction[END_REF][START_REF] Il'in | Matching of asymptotic expansions of solutions of boundary value problems[END_REF][START_REF] Maz'ya | Asymptotic theory of elliptic boundary value problems under a singular perturbation of the domains[END_REF]. An important component of small-inhomogeneity asymptotics of displacement fields is the elastic moment tensor (EMT), whose definition and properties are studied in e.g. [START_REF] Ammari | Complete asymptotic expansions of solutions of the system of elastostatics in the presence of an inclusion of small diameter and detection of an inclusion[END_REF] for isotropic materials and [START_REF] Beretta | Small volume asymptotics for anisotropic elastic inclusions[END_REF][START_REF] Nazarov | Polarization matrices in anisotropic heterogeneous elasticity[END_REF] for anisotropic materials.

However, some optimization problems involve stress-based cost functionals, which is equivalent (for the present context of linearly elastic solids) to considering functionals that depend on strains or displacement gradients. Examples include topology optimization of composite structures with materials constrained by yield criteria, and flaw identification using full-field kinematical measurements, as both types of problems may be cast as minimization problems involving stress-based cost functionals. The asymptotic behavior of such cost functionals is quite different, and more involved, than in the previous case, due to the fact that the strain perturbation inside an elastic inhomogeneity has a finite, nonzero limit as a → 0, while the asymptotics of displacement-based functionals rests upon the fact that the magnitude of the displacement perturbation vanishes as a → 0.

So far, only few works have investigated the small-inhomogeneity asymptotics of stressbased functionals. The 2D isotropic case is addressed for specific stress-based functionals (elastic energy, von Mises and Drucker-Prager yield criteria) in [START_REF] Amstutz | Topological derivative-based topology optimization of structures subject to Drucker-Prager stress constraints[END_REF][START_REF] Amstutz | Topological optimization of structures subject to Von Mises stress constraints[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF]. Moreover, results for general quadratic stress-dependent functionals are given in [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF] within a 2D and 3D anisotropic framework. The special case of cost functionals that are directly linked to the anisotropic elastic potential energy in 2D or 3D is considered in [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF]. The main purpose of this article is to formulate and justify the topological derivative for stress-dependent functionals, in a quite general three-dimensional framework where the functionals are defined in terms of domain integrals of arbitrary sufficiently smooth densities while both the background material and the small trial inhomogeneity are allowed to have anisotropic elastic properties. The results given herein for topological derivatives thus contain expressions given in [START_REF] Amstutz | Topological derivative-based topology optimization of structures subject to Drucker-Prager stress constraints[END_REF][START_REF] Amstutz | Topological optimization of structures subject to Von Mises stress constraints[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF] as special cases. They are illustrated by computational experiments on 2D and 3D examples involving anisotropic elasticity and stress-based quadratic or non-quadratic cost functionals, and inspired by topology optimization or flaw identification.

The article is organised as follows. Section 2 recalls the concept of topological derivative, introduces notation and collects the main facts on the elastic transmission problem and elastic moment tensor. The main result on topological derivative is stated and established in Section 3. The numerical evaluation of DJ is addressed in Section 4, and Section 5 is then devoted to the presentation and discussion of computational experiments.

Elastic transmission problem and cost functional 2.1 Notation, elastic transmission problem

Consider an elastic body occupying a smooth bounded domain Ω ⊂ R 3 . The anisotropic elastic properties of the background material (against which the effect of small inhomogeneities will be considered), assumed to be homogeneous, are characterized by the fourth-order elasticity tensor C. The boundary ∂Ω is split according to ∂Ω = Γ D ∪ Γ N (where Γ D ∩ Γ N = ∅ and |Γ D | = 0), so that a given force density g ∈ L 2 (Γ N ; R 3 ) is applied on Γ N while a given displacement ū ∈ H 1/2 (Γ D ; R 3 ) is prescribed on Γ D . Additionally, a body force density f assumed (for reasons given later) to have C 0,β (Ω) Hölder continuity for some β > 0 is applied to Ω.

The background solution, i.e. the displacement field arising in the reference solid due to the prescribed excitations (f , g, ū), is defined as the solution to div(C :

ε[u]) + f = 0 in Ω, (C : ε[u])•n = g on Γ N , u = ū on Γ D ( 1
)
where n is the unit outward normal to Ω and ε[w] denotes the linearized strain tensor associated with a given displacement w, defined by

ε[w] = 1 2 (∇w + ∇w T ). (2)
In [START_REF] Allaire | Structural optimization using topological and shape sensitivity via a level-set method[END_REF] and hereinafter, symbols ' • ' and ' : ' denote single and double inner products, e.g. (C : ε) ij = C ijk ε k , with Einstein's convention of summation over repeated indices implicitly used throughout. Alternatively, the background displacement is governed by the weak formulation

Find u ∈ W ( ū), u, w C Ω = F (w), ∀w ∈ W 0 , (3)
where u, w C D denotes the bilinear elastic energy form associated to given domain D ⊂ R 3 and elasticity tensor C, i.e.:

u, w C D := D ε[u] : C : ε[w] dV = D ∇u : C : ∇w dV (4)
(with the second equality holding by virtue of the well-known minor symmetries of C), the linear form F associated to the loading is defined by

F (w) = Ω f •w dV + Γ N g•w dS, (5)
and having introduced, for given ū ∈ H 1/2 (Γ D ; R 3 ), the spaces W ( ū) and W 0 of displacement fields that are kinematically admissible with respect to arbitrary and homogeneous prescribed Dirichlet data, respectively, i.e.:

W ( ū) := v ∈ H 1 (Ω; R 3 ), v = ū on Γ D , W 0 := W (0). (6)
The C 0,β (Ω; R 3 ) regularity assumption on f , which is stronger than the more-usual assumption f ∈ L 2 (Ω; R 3 ), ensures (e.g. from the properties of elastic volume potentials, see [START_REF]Three-dimensional problems of the mathematical theory of elasticity and thermoelasticity[END_REF], Thm. 10.4) that u is in C 2,β (D; R 3 ) for any subset D Ω. It is made in order to later permit Taylor expansions of displacements or strains about selected internal points.

Transmission problem for a small trial inhomogeneity

Now, consider a single small elastic inhomogeneity located at z ∈ Ω, of characteristic linear size a, occupying the domain

B a = z + aB,
where B ⊂ R 3 is a bounded smooth domain and a is small enough to have Ba Ω. The inhomogeneity has anisotropic properties characterized by the elasticity tensor C , so that the elastic properties of the whole solid are defined by the tensor-valued field C a given by

C a = (1 -χ(B a ))C + χ(B a )C = C + χ(B a )∆C, (7)
χ(D) being the characteristic function of the domain D and ∆C := C -C denoting the elastic tensor perturbation.

The displacement field u a ∈ W ( ū) arising in the solid containing the small inhomogeneity due to the prescribed excitations (f , g, ū) solves the transmission problem

div(C a : ε[u a ]) + f = 0 in Ω, (C : ε[u a ])•n = g on Γ N , u = ūa on Γ D . (8)
Formulation [START_REF] Amstutz | Topological optimization of structures subject to Von Mises stress constraints[END_REF] implicitly enforces, by virtue of its distributional interpretation, the perfectbonding relations

u a | + = u a | -and (C : ε[u a ])•n| + = (C : ε[u a ])•n| -on ∂B a
, where the ± subscripts indicate limiting values from outside and inside B a , respectively, and n is the unit outward normal vector to ∂B a .

The transmission problem (8) can alternatively be formulated in terms of the displacement perturbation v a := u a -u rather than the total displacement u a . Subtracting (3) from the corresponding weak formulation of (8) yields the following weak formulation for v a :

Find v a ∈ W 0 , v a , w Ca Ω = -u, w ∆C Ba , ∀w ∈ W 0 . (9)
Free-space transmission problem (FSTP). The auxiliary problem of a perfectly-bonded inhomogeneity (B, C ) embedded in an infinite elastic medium Ω ∞ = R 3 and subjected to a uniform remote stress equal to the background stress at z will play an important role in the sequel. The FSTP thus consists in finding the displacement field u B such that

div(C B : ε[u B ]) = 0 in R 3 , u B (ξ) -u ∞ (ξ) = O(|ξ| -2 ), |ξ| → ∞, (10)
where the background displacement u ∞ is defined by u ∞ (ξ) = ∇u(z) • ξ and with C B := C + χ(B)∆C. The FSTP [START_REF] Bay | Digital volume correlation: threedimensional strain mapping using X-ray tomography[END_REF] can be recast into the following weak formulation for the displacement perturbation

v B := u B -u ∞ : Find v B ∈ W ∞ , v B , w C B R 3 = -u ∞ , w ∆C B , ∀w ∈ W ∞ , (11)
with the function space W ∞ defined by

W ∞ = w ∈ L 2 loc (R 3 ; R 3 ), ∇w ∈ L 2 (R 3 ; R 3×3 ) . When B is an ellipsoid, i.e. B = {x ∈ R 3 , (x 1 /a 1 ) 2 + (x 2 /a 2 ) 2 + (x 3 /a 3 ) 2 ≤ 1} , a 1 ,
a 2 , a 3 > 0 (with the axes of the Cartesian frame (Ox 1 x 2 x 3 ) assumed, without loss of generality, to be aligned with the ellipsoid principal axes), the FSTP ( 12) is analytically solved in Eshelby's celebrated paper [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion and related problems[END_REF]. In that case, v B is found to have a constant strain and gradient inside B, with

∇v B (x) = S int : ε (z), ε (z) = -(C + ∆C : S int ) -1 : ∆C : ∇u(z) (x ∈ B), (12)
S int being the interior Eshelby tensor of B, which relates a constant eigenstrain ε in B to the constant displacement gradient in B. We note that the above definition of S int somewhat differs from the usual Eshelby tensor S, which is such that ε

[v B ](x) = S : ε (z) (x ∈ B)
instead. Using e.g. Eq. (17.14) of [START_REF] Mura | Micromechanics of Defects in Solids[END_REF], one has

S int ijmn = 1 4π C k mn S 2 N ik (ξ(ŝ))ξ (ŝ)ξ j (ŝ) dS(ŝ), ( 13
)
where S 2 is the unit sphere, ξ(ŝ) is defined for ŝ ∈ S 2 by

ξ = (ŝ 1 /a 1 , ŝ2 /a 2 , ŝ3 /a 3 ), (14)
and the tensor-valued function ξ → N (ξ) is given by N (ξ) = K -1 (ξ), where K ik (ξ) := C ijk ξ j ξ l is the Christoffel acoustic tensor. A similarly modified version of the exterior Eshelby tensor will be used in Section 4.1.

Elastic moment tensor. The elastic moment tensor (EMT) [START_REF] Ammari | Complete asymptotic expansions of solutions of the system of elastostatics in the presence of an inclusion of small diameter and detection of an inclusion[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Nazarov | Polarization matrices in anisotropic heterogeneous elasticity[END_REF] will be seen to play an important role in the small-inhomogeneity asymptotics of cost functionals. The EMT A is the fourth-order tensor defined for any value of the constant tensor ∇u(z) ∈ R 3×3 by

A : ∇u(z) = B ∆C : ∇u B dV = B ∆C : (∇u(z) + ∇v B ) dV, ( 15
)
where v B is the solution of the FSTP [START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF]. The EMT has the same symmetry properties as the elasticity tensor C: for any pair of second-order tensors E, E ∈ R 3×3 , A satisfies

E : A : E = E : A : E (major symmetry), (16a) E : A : E = E : A : E T and E : A : E = E T : A : E (minor symmetries). ( 16b
)
The EMT for an ellipsoidal inhomogeneity is readily found, using ( 12) into [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF], to have the closed form expression [START_REF] Ammari | Complete asymptotic expansions of solutions of the system of elastostatics in the presence of an inclusion of small diameter and detection of an inclusion[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Nazarov | Polarization matrices in anisotropic heterogeneous elasticity[END_REF][START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF]]

A = |B| C : (C + ∆C : S int ) -1 : ∆C. ( 17
)
Asymptotic behavior of v a [START_REF] Allaire | Damage and fracture evolution in brittle materials by shape optimization methods[END_REF][START_REF] Beretta | Small volume asymptotics for anisotropic elastic inclusions[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF]. The inner approximation ṽa of v a is given by ṽa

(x) = av B x -z a , x ∈ B a ( 18
)
where v B solves the FSTP [START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF]. Moreover, for any cut

-off function θ ∈ C ∞ c (Ω) such that θ ≡ 1 in a neighborhood D of z, let δ a ∈ H 1 (Ω; R 3 ) be defined by v a = θ ṽa + δ a . (19)
Then, there exists a constant C > 0 independent of a such that the estimates

(a) ∇ṽ a L 2 (Ω) ≤ Ca 3/2 , (b) ṽa L 2 (Ω) ≤ Ca 5/2 , (c) δ a H 1 (Ω) ≤ Ca 5/2 (20)
and (a)

∇v a L 2 (Ω) ≤ Ca 3/2 , (b) v a L 2 (Ω) ≤ Ca 5/2 . ( 21
)
hold [START_REF] Allaire | Damage and fracture evolution in brittle materials by shape optimization methods[END_REF][START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF]. Estimate (20c) also relies on the previously-made assumption that f has C 0,β regularity for some β > 0 in a neighbourhood of z. Moreover, the known O(| x| -3 ) far-field behavior of ∇v B [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Mura | Micromechanics of Defects in Solids[END_REF] implies, by rescaling, that

∇ṽ a L ∞ (Ω\D) ≤ Ca 3 . ( 22
)
Remark 1. As mentioned in [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF], the above estimates, established assuming C to be constant (homogeneous background material), are expected to also hold for heterogeneous elastic properties that are smooth in a fixed neighbourhood of z (with the EMT then defined in terms of C(z)). A numerical experiment involving a piecewise constant C for non-destructive testing is shown in Section 5.

Cost functional

Cost functionals of the form

J(C a ) = J a (u a , ∇u a ) with J a (u, d) = Ω ψ a x, u(x), d(x) dV (x) (23)
are considered, where the density ψ a : Ω × R 3 × R 3×3 → R is defined by

ψ a = (1 -χ(B a ))ψ + χ(B a )ψ = ψ + χ(B a )∆ψ, (24)
with functions ψ and ψ (and hence also ∆ψ := ψψ) assumed to be twice differentiable in all their arguments. Moreover, all second-order derivatives of ψ and ψ are assumed to have C 0,γ (Ω × R 3 × R 3×3 ) Hölder regularity for some γ > 0 with respect to all their arguments. We will denote by x ∈ Ω, u ∈ R 3 , d ∈ R 3×3 the generic arguments of a density ψ(x, u, d). Then, ∂ x ψ, ∂ u ψ, ∂ d ψ will denote the partial derivatives with respect of the corresponding arguments, with higher-order partial derivatives denoted similarly, e.g. ∂ 2 ud ψ (with similar notations for ψ ). The assumed Hölder regularity of ψ can thus be expressed as the existence of a finite constant K such that, for every combination of indices i, j = x, u, d , sup

x 0 ∈Ω u 0 ∈R 3 d 0 ∈R 3×3 ∂ 2 ij ψ(x 0 , u 0 , d 0 ) + sup (x 0 ,x 1 )∈Ω (u 0 ,u 1 )∈R 3 (d 0 ,d 1 )∈R 3×3 ∂ 2 ij ψ(x 0 , u 0 , d 0 ) -∂ 2 ij ψ(x 1 , u 1 , d 1 ) |x 0 -x 1 | γ + |u 0 -u 1 | γ + |d 0 -d 1 | γ ≤ K,
with a similar inequality holding for ψ . We note that the above assumption implies that ∂ 2 ij ψ and ∂ 2 ij ψ are bounded over Ω × R 3 × R 3×3 . The partial directional derivatives ∂ u J a and ∂ d J a of J a with respect to its first and second arguments are defined, for later use, by

∂ u J a (u, d), w = Ω\Ba ∂ u ψ(x, u, d)•w dV + Ba ∂ u ψ (x, u, d)•w dV, w ∈ W 0 , ∂ d J a (u, d), h = Ω\Ba ∂ d ψ(x, u, d) : h dV + Ba ∂ d ψ (x, u, d) : h dV, h ∈ L 2 (Ω; R 3×3 ). ( 25
)
Remark 2. The assumed Hölder and boundedness conditions on the cost functional densities might seem restrictive, but are satisfied by a number of cost functions often used in applications (compliance, energy based functionals, least squares misfit functionals, yield functions). Examples include (i) all quadratic functions of d, and (ii) the useful penalization function

ψ(d) = Ψ n (q(d)), Ψ n (t) := 1 + t n 1/n -1 (26)
(where q is any quadratic function of d), introduced in [START_REF] Amstutz | Topological derivative-based topology optimization of structures subject to Drucker-Prager stress constraints[END_REF] to approximate pointwise yield functions with thresholds, whose derivatives Ψ n and Ψ n are bounded over R + (one finds e.g. that 0 < |Ψ n | ≤ 2 1/n (n -1) for any t > 0, n > 1). Case (ii) is typically used with q(∇u) chosen as the appropriate squared (von Mises, Hill-Tsai, Drucker-Prager. . . ) yield function, normalized so that the corresponding criterion reads q(∇u) ≤ 1.

Topological derivative

In this section, we state our main result (Theorem 1) and give its proof in Section 3.1.

The generalization to piecewise-regular densities ψ is then addressed in Section 3.2. Finally, particular instances of Theorem 1 are discussed in Section 3.3.

Definition 1 (topological derivative). Assume that J(C a ) can be expanded in the form

J(C a ) = J(C) + δ(a)DJ(z) + o(δ(a)) (27)
where δ(a) is assumed to vanish as a → 0 and characterizes the small-inhomogeneity asymptotic behavior of J(C a ). Then, the coefficient DJ(z), which also depends a priori on the shape B and the moduli C, C , is called the topological derivative of J at z ∈ Ω.

Remark 3. Terminology for the concept of topological derivative varies, with "gradient" or "sensitivity" sometimes used instead of "derivative".

Theorem 1. Assume a three-dimensional setting as laid out in Section 2. Any cost functional J of the form [START_REF] Gay | Composite materials: design and applications[END_REF] and fulfilling the assumptions made in Section 2.3 admits an expansion of the form [START_REF] Il'in | Matching of asymptotic expansions of solutions of boundary value problems[END_REF], with δ(a) = a 3 and the topological derivative DJ(z) of J at z ∈ Ω given by DJ(z) = |B|∆ψ(z, u(z), ∇u(z)) -∇p(z) : A : ∇u(z)

+ ∂ d (∆ψ)(z, u(z), ∇u(z)) : B ∇v B ( x) dV ( x) + R 3 \B G(z, ∇v B ( x)) dV ( x) + B G (z, ∇v B ( x)) dV ( x). ( 28
)
The functions G and G : R 3 × R 3×3 → R are defined, for a given background solution u, by

G(z, d) := ψ(z, u(z), ∇u(z) + d) -ψ(z, u(z), ∇u(z)) -∂ d ψ(z, u(z), ∇u(z)) : d (29a) G (z, d) := ψ (z, u(z), ∇u(z) + d) -ψ (z, u(z), ∇u(z)) -∂ d ψ (z, u(z), ∇u(z)) : d, (29b)
and p ∈ W 0 is the adjoint state, defined as the solution of the weak formulation p, w

C Ω = ∂ u J 0 (u, ∇u), w + ∂ d J 0 (u, ∇u), ∇w ∀w ∈ W 0 , (30)
with ∂ u J 0 and ∂ u J 0 as defined by [START_REF] Hinton | Fully stressed topological design of structures using an evolutionary procedure[END_REF].

When the densities ψ, ψ are linear or quadratic in their third argument (i.e. when ∂ 2 dd ψ and ∂ 2 dd ψ are independent on d), letting D(z) := ∂ 2 dd ψ(z, u(z)) and D (z) := ∂ 2 dd ψ (z, u(z)), the last two terms in [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF] are given by the more explicit expression

1 2 R 3 \B ∇v B ( x) : D(z) : ∇v B ( x) dV ( x) + 1 2 B ∇v B ( x) : D (z) : ∇v B ( x) dV ( x)
Moreover, under two-dimensional plane-strain conditions (where only in-plane displacements are nonzero), the result (28) still holds (with Ω ⊂ R 2 , B ⊂ R 2 and the next-to-last integral now taken over R 2 \ B), while δ(a) = a 2 in expansion [START_REF] Il'in | Matching of asymptotic expansions of solutions of boundary value problems[END_REF].

Remark 4. The result (28) for DJ involves the gradient ∇v B of the FSTP solution, rather than just its strain ε[v B ], which explains the definition (10) chosen here for the Eshelby tensor.

Proof of Theorem 1

The detailed proof to follow concentrates on the 3D case, its adaptation to the 2D case being then outlined in a comment. The proof consists in finding the leading contribution to the difference J(C a ) -J(C) = J a (u a , ∇u a ) -J 0 (u, ∇u) as a → 0. To this end, we write

J a (u a , ∇u a ) -J 0 (u, ∇u) = J a (u, ∇u) -J 0 (u, ∇u) + J a (u a , ∇u a ) -J a (u, ∇u a ) -∂ u J a (u, ∇u), v a + J a (u, ∇u a ) -J a (u, ∇u) -∂ d J a (u, ∇u), ∇v a + ∂ u J a (u, ∇u), v a + ∂ d J a (u, ∇u), ∇v a , (31)
with ∂ u J a and ∂ u J a as defined by [START_REF] Hinton | Fully stressed topological design of structures using an evolutionary procedure[END_REF], and separately evaluate the leading contribution of each bracketed term in the right-hand side of (31); this is done in the following Lemmas 2 to 5. Using the results of the lemmas in the above decomposition then directly establishes both the expansion [START_REF] Il'in | Matching of asymptotic expansions of solutions of boundary value problems[END_REF] and the expression (28) of DJ(z) stated in Theorem 1. Lemma 2. Let ∆ψ be defined as in [START_REF] Gradshteyn | Tables of integrals, series and products[END_REF]. One has

J a (u, ∇u) -J 0 (u, ∇u) = a 3 |B|∆ψ(z, u(z), ∇u(z)) + o(a 3
). Proof. By interior regularity for u and the assumed smoothness of ψ, ψ , ∆ψ(x, u(x), ∇u(x)) is continuous at x = z. Therefore, the Lemma follows easily from:

J a (u, ∇u) -J 0 (u, ∇u) = Ω ψ a (x, u, ∇u) -ψ(x, u, ∇u) dV (x) = Ba ∆ψ(x, u, ∇u) dV (x) = a 3 |B|∆ψ(z, u(z), ∇u(z)) + o(a 3 ),
where the last step exploits the fact that the volume of B a is |B a | = a 3 |B|. Lemma 3. Let the displacement perturbation v a solve problem [START_REF] Baumgartner | SKO (soft kill option): the biological way to find an optimum structure topology[END_REF]. One has

J a (u a , ∇u a ) -J a (u, ∇u a ) -∂ u J a (u, ∇u), v a = o(a 3 ), with ∂ u J a as defined by (25)
Proof. The proof is based on Taylor expansions. A first-order expansion of ψ with respect to its first argument first yields

ψ(x, u a , ∇u a ) -ψ(x, u, ∇u a ) = ∂ u ψ(x, u, ∇u a )•v a + 1 2 v a •∂ 2 uu ψ(x, u + δ u v a , ∇u a )•v a for some δ u (x) ∈ [0, 1]
. Moreover, a zeroth-order Taylor expansion of ∂ u ψ with respect to its second argument gives

∂ u ψ(x, u, ∇u a )•v a -∂ u ψ(x, u, ∇u)•v a = ∇v a : ∂ 2 du ψ(x, u, ∇u + δ d ∇v a )•v a for some δ d (x) ∈ [0, 1].
Both expansions are valid due to the assumed regularity of ψ. Similar expansions also hold for the density ψ , for some δ u (x), δ d (x) ∈ [0, 1]. Combining all of these expansions, one finds

J a (u a , ∇u a ) -J a (u, ∇u a ) = ∂ u J a (u, ∇u), v a + R a (32)
with the remainder R a given by

R a = Ω\Ba 1 2 v a •∂ 2 uu ψ(x, u + δ u v a , ∇u a )•v a + ∇v a : ∂ 2 du ψ(x, u, ∇u + δ d ∇v a )•v a dV + Ba 1 2 v a •∂ 2 uu ψ (x, u + δ u v a , ∇u a )•v a + ∇v a : ∂ 2 du ψ (x, u, ∇u + δ d ∇v a )•v a dV
Next, thanks to the boundedness of the second-order partial derivatives of ψ and ψ , there exists a constant C > 0 such that [START_REF]Freeyams[END_REF] imply that there exists a constant C > 0 such that

R a ≤ C( v a 2 L 2 (Ω) + ∇v a L 2 (Ω) v a L 2 (Ω) ) Finally, estimates
R a ≤ Ca 4 = o(a 3 ).
Using the above estimate in (32) completes the proof. Lemma 4. One has

J a (u, ∇u a ) -J a (u, ∇u) -∂ d J a (u, ∇u), ∇v a = a 3 R 3 \B G(z, ∇v B ( x)) dV ( x) + B G (z, ∇v B ( x)) dV ( x) + o(a 3 ) (33)
where ∂ d J a is defined by [START_REF] Hinton | Fully stressed topological design of structures using an evolutionary procedure[END_REF] and the functions G and G : R 3 ×R 3×3 → R are defined by (29a,b).

Proof. The combination to be estimated is first recast in the form

J a (u, ∇u a ) -J a (u, ∇u) -∂ d J a (u, ∇u), ∇v a = Ω\Ba G(x, ∇v a ) dV (x) + Ba G (x, ∇v a ) dV (x) = J a (u, ∇u a ) -J a (u, ∇u + ∇ṽ a ) -∂ d J a (u, ∇u + ∇ṽ a ), ∇v a + ∂ d J a (u, ∇u + ∇ṽ a ), ∇v a -∂ d J a (u, ∇u), ∇v a + J a (u, ∇u + ∇ṽ a ) -J a (u, ∇u) -∂ d J a (u, ∇u), ∇ṽ a (34a)
having used functions G, G defined by (29a,b) for the last equality, and with va defined (using [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] in the second equality below) by

va = v a -ṽa = δ a + (θ -1)ṽ a . (34b)
Since ∇v a = ∇δ a + (θ -1)∇ṽ a + ∇θ⊗ ṽa , we note for later use that estimates [START_REF]Freefem++[END_REF] and ( 22) (together with the fact that the support of (θ -1)∇ṽ a is Ω \ D) imply

∇v a L 2 (Ω) = O(a 5/2 ). ( 34c
)
We first focus on contributions of integrals over Ω \ B a , i.e. of the density G, to (34a). To begin, each bracketed combination in (34a) is reformulated by exploiting Taylor expansions (of first-order with respect to its third argument, with integral remainder) of ψ, to obtain

Ω\Ba G(x, ∇v a ) dV (x) = Ω\Ba ∇v a (x) : D 2 x, u + ṽa , ∇v a (x) : ∇v a (x) dV (x) + Ω\Ba ∇v a (x) : D 1 x, ∇v a (x) : ∇ṽ a (x) dV (x) + Ω\Ba ∇ṽ a (x) : D 3 x, ∇ṽ a (x) : ∇ṽ a (x) dV (x), (34d)
with

D 1 (y, d) = 1 0 ∂ dd ψ y, u(y), ∇u(y) + td) dt, D 2 (y, w, d) = 1 0 ∂ dd ψ y, u(y), ∇w(y) + td (1 -t) dt, D 3 (y, d) = D 2 (y, u, d). (34e)
Both D 1 x, ∇v a (x) and D 2 x, u(x)+ṽ a (x), ∇v a (x) are bounded over Ω, due to the boundedness of the second-order partial derivatives of ψ. This remark is exploited by applying the Cauchy-Schwarz inequality to the first two integrals I 1 and I 2 of the right-hand side of (34d) and invoking estimates (20a) and (34c), to obtain

I 1 = O(a 5 ) = o(a 3 ) and I 2 = O(a 4 ) = o(a 3 ).
Using these estimates, applying the change of variables x = (x -z)/a to the third integral of the right-hand side of (34d) (whereby dV (x) = a 3 dV ( x)), and recalling definition [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] of ṽa , we obtain

Ω\Ba G(x, ∇ṽ a ) dV (x) = a 3 ((Ω-z)/a)\B ∇v B ( x) : D 3 (z + a x, ∇v B ( x)) : ∇v B ( x) dV ( x) + o(a 3 ) = a 3 ((Ω-z)/a)\B ∇v B ( x) : D 3 (z, ∇v B ( x)) : ∇v B ( x) dV ( x) + R + o(a 3 ) = a 3 ((Ω-z)/a)\B G(z, ∇v B ( x)) dV ( x) + R + o(a 3 ) (34f) (having noted that ∇v B ( x) : D 3 (z, ∇v B ( x)) : ∇v B ( x) = G(z, ∇v B ( x)))
, where the remainder R is such that

R := a 3 ((Ω-z)/a)\B ∇v B ( x) : D 3 (z + a x, ∇v B ( x)) -D 3 (z, ∇v B ( x)) : ∇v B ( x) dV ( x) ≤ Ca 3+γ ((Ω-z)/a)\B |∇v B ( x)| 2 | x| γ dV ( x) ≤ Ca 3+γ R 3 \B |∇v B ( x)| 2 | x| γ dV ( x), (34g)
by virtue of the inequality

|D 3 (z + a x, d) -D 3 (z, d)| ≤ Ca γ | x| γ (34h)
stemming from the assumed C 0,γ Hölder regularity of ∂ dd ψ and the known C 2 interior regularity of u in Ω, which implies that there exists τ , τ ∈ [0, 1] such that

u(x) -u(z) = a∇u(z + τ a x)• x, ∇u(x) -∇u(z) = a∇ 2 u(z + τ a x)• x. ( 34i
)
The known O(| x| -3 ) far-field behavior of ∇v B ( x) [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF][START_REF] Mura | Micromechanics of Defects in Solids[END_REF] implies that the last integral in (34g) over the unbounded domain R 3 \B is finite for any γ < 3, and hence that R = O(a 3+γ ) = o(a 3 ). Finally, taking the limit ((Ω-z)/a)\B → R 3 \B in (34f) (which is legitimate by the dominated convergence theorem since |G(z, ∇v B ( x))| is integrable over R 3 \ B), the desired asymptotic form of (34d) is obtained:

Ω\Ba G(x, ∇v a ) dV (x) = a 3 R 3 \B G(z, ∇v B ( x)) dV ( x) + o(a 3 ). ( 34j
)
The second integral in the right-hand side of (34a) can be estimated following similar arguments. A representation similar to (34d) holds, with integrals taken over B a and ψ replaced by ψ in (34e). Noting in addition that now va = δ a since θ = 1 in B a , the first two integrals in the right-hand side of the counterpart of (34d) are easily established to be of order o(a 3 ) using estimates [START_REF]Freefem++[END_REF]. Using again the change of variables x = (x -z)/a in the remaining integral, one then finds

Ba G (x, ∇v a ) dV (x) = a 3 B G (z, ∇v B ( x)) dV ( x) + R + o(a 3 ), (34k)
where, exploiting through (34h) the assumed Hölder regularity of ψ , the remainder R is such that

R := a 3 B ∇v B ( x) : D(z + a x, ∇v B ( x)) -D(z, ∇v B ( x)) : ∇v B ( x) dV ( x) ≤ Ca 3+γ B |∇v B ( x)| 2 | x| γ dV ( x) = O(a 3+γ ).
The desired asymptotic form of (34k) is therefore obtained:

Ba G (x, ∇v a ) dV (x) = a 3 B G (z, ∇v B ( x)) dV ( x) + o(a 3 ) (34l)
The lemma finally follows from using expansions (34j) and (34l) in (34a).

Finally, the leading contribution to the last bracketed combination of ( 31) is given in the following lemma in terms of an adjoint solution. Lemma 5. Let the adjoint solution p ∈ W 0 be defined by the weak formulation p, w

C Ω = ∂ u J 0 (u, ∇u), w + ∂ d J 0 (u, ∇u), ∇w ∀w ∈ W 0 , (35)
with ∂ u J 0 and ∂ d J 0 as defined by [START_REF] Hinton | Fully stressed topological design of structures using an evolutionary procedure[END_REF]. One has

∂ u J a (u, ∇u), v a + ∂ d J a (u, ∇u), ∇v a = a 3 |B| -∇p(z) : A : ∇u(z) + ∂ d (∆ψ)(z, u(z), ∇u(z)) : B ∇v B dV ( x) + o(a 3 ) (36)
Proof. Setting ∆ψ := ψψ, one has

∂ u J a (u, ∇u), v a + ∂ d J a (u, ∇u), ∇v a = ∂ u J 0 (u, ∇u), v a + ∂ d J 0 (u, ∇u), ∇v a + ∂ u ∆J(u, ∇u), v a + ∂ d ∆J(u, ∇u), ∇v a (37a)
Invoking the definition [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF] of the adjoint solution, the identity

v a , w C Ω = -u a , w
∆C Ba , ∀w ∈ W 0 verified by the transmission problem, one finds

∂ u J 0 (u, ∇u), v a + ∂ d J 0 (u, ∇u), ∇v a = p, v a C Ω = -u a , p ∆C Ba = -a 3 |B|∇p(z) : A : ∇u(z) + o(a 3 ), (37b)
where the last equality holds by virtue of [15, Lemma 2] and p having C 2 interior regularity (see Section 3.2). Next, using decomposition [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] of v a and the fact that ∆J has B a as its geometrical support (implying in particular that θ = 1 in B a ), one has

∂ u ∆J(u, ∇u), v a + ∂ d ∆J(u, ∇u), ∇v a = ∂ u ∆J(u, ∇u), (ṽ a + δ a ) + ∂ d ∆J(u, ∇u), ∇(ṽ a + δ a ) (37c)
The partial derivatives of ∆ψ(x, u(x), ∇u(x)) being bounded by virtue of the assumptions made on ψ, ψ and the C 2 interior regularity of u, there exists a constant C > 0 such that

∂ u ∆J(u, ∇u), (ṽ a + δ a ) + ∂ d ∆J(u, ∇u), ∇δ a ≤ a 3/2 ( ṽa L 2 (Ω) + δ a H 1 (Ω) )) ≤ Ca 4 , (37d)
the last inequality stemming from estimates [START_REF]Freefem++[END_REF].

The term ∂ d ∆J(u, ∇u), ∇ṽ a remains to be estimated. By the mean value theorem applied to ∂ d (∆ψ)(x, u(x), ∇u(x)), there exists t(x) ∈ [0, 1] such that

∂ d (∆ψ)(x, u(x), ∇u(x)) = ∂ d (∆ψ)(z, u(z), ∇u(z)) + ∂ 2 xd ∆ψ(x t , u t , ∇u t )•(x -z) + ∂ 2 ud ∆ψ(x t , u t , ∇u t )•[u(x) -u(z)] + ∂ 2 dd ∆ψ(x t , u t , ∇u t ) : [∇u(x) -∇u(z)
] where x t , u t , ∇u t are defined by

x t := z + t(x -z), u t := u(z) + t[u(x) -u(z)], ∇u t := ∇u(z) + t[∇u(x) -∇u(z)]. (37e)
Introducing xz = a x and expansions (34i), stemming from the C 2 interior regularity of u, in the above definitions, one obtains for

x ∈ B a ∂ d (∆ψ)(x, u(x), ∇u(x)) = ∂ d (∆ψ)(z, u(z), ∇u(z)) + O(a)
which in turn implies

∂ d ∆J(u, ∇u), ∇ṽ a = a 3 ∂ d (∆ψ)(z, u(z), ∇u(z)) : B ∇v B dV + o(a 3 ) (37f)
The lemma finally follows by substituting (37d) and (37f) into the right-hand side of (37c) and then using the resulting estimate together with (37b) in (37a).

The two-dimensional case. The proof for the two-dimensional plane-strain case is identical, except for the fact that estimates (20) to ( 22) must be replaced by their following two-dimensional counterparts:

(a) ∇ṽ a L 2 (Ω;R 2 ) ≤ Ca, (b) ṽa L 2 (Ω;R 2 ) ≤ Ca 2 | log a|, (c) δ a H 1 (Ω;R 2 ) ≤ Ca 2 , ( 38a
) (a) ∇v a L 2 (Ω;R 2 ) ≤ Ca, (b) v a L 2 (Ω;R 2 ) ≤ Ca 2 . (38b) ∇ṽ a L ∞ (Ω\D) ≤ Ca 2 , (38c)
which can be established e.g. by adapting to the two-dimensional case the proofs given in [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF].

Case of piecewise-regular cost functional densities

Now we extend the previous results to the topological derivative of functionals J ω defined by an integral over a portion ω ⊂ Ω of the elastic body Ω, of the form

J ω (C a ) = J a (u a , ∇u a ; ω), with J a (u, d; ω) = ω ψ a (x, u, ∇u) dV (x),
where the trial inhomogeneity B a is assumed to satisfy either B a ω or B a (Ω\ ω) (the case where B a ∩ ∂ω = ∅ not being considered). If B a ω, the previous analysis remains valid, with DJ still given by ( 28) and the only change concerning the adjoint solution, which now satisfies the weak formulation

Ω ε[p] : C : ε[q]dV = Ω χ ω ∂ d ψ : ∇q + ∂ u ψ • q dV, ∀q ∈ W 0 (Ω). ( 39
)
On the other hand, if B a (Ω \ ω), ψ a -ψ = 0 in ω and the cut-off function θ in decomposition [START_REF] Feijóo | A new method in inverse scattering based on the topological derivative[END_REF] can be chosen, for any sufficiently small a, such that θ = 0 in ω. This choice implies that v a = δ a in ω for any z ∈ Ω\ω, and hence, by estimate (20c), that

v a H 1 (ω) ≤ Ca 5/2 .
Consequently, retracing the proof of Theorem 1, contributions to DJ arising from ∂ dd ψ and ∂ dd ψ in Lemma 4 are o(a 3 ), and DJ is simply given, in terms of the solutions u of (3) and p of (39), by DJ(z) = -∇p(z) : A : ∇u(z).

Regularity of the adjoint solution. As it was previously seen, the point-wise evaluation of DJ at some z ∈ Ω requires the background displacement u and the adjoint solution p to have some local regularity at z, namely u, p ∈ C 2,α (D; R 3 ) for some neighborhood D Ω of z and α ∈ (0, 1). The needed regularity for u follows directly from the regularity of the body force density f ∈ C 0,α (Ω; R 3 ). When the cost function J depends only on u, the adjoint state

p solves Find p ∈ W 0 , p, w C Ω = ∂ u J 0 (u), w , ∀w ∈ W 0 . Then, if ∂ u ψ(x, u(x)) ∈ C 0,α (Ω; R 3
), p fulfills automatically the required interior regularity in Ω. On the contrary, the case when the cost functional depends on ∇u is slightly more delicate. In such a case p solves

Find p ∈ W 0 , p, w C Ω = ∂ d J 0 (∇u), ∇w , ∀w ∈ W 0 Thus if ∂ d ψ(x, ∇u(x)) ∈ C 1,α (Ω; R 3×3
), p fulfills the required regularity. We remark that the components of the tensor (∇u) ij ∈ C 1,α (Ω) so we can take any function ψ(•, d) ∈ C 1 (R 3×3 ) for this purpose. Finally when the cost function J is defined in a sub-domain ω of Ω,

∂ d J 0 (∇u), ∇w = ω ∂ d ψ(x, ∇u) : ∇w dV = Ω χ ω ∂ d ψ(x, ∇u(x)) : ∇w dV, the adjoint state p solves      -div(Cε(p)) = -div(∂ d ψ(x, ∇u)) in D, p = 0 on Γ D , C : ε(p)•n = 0 on Γ N ,
where Γ N ∩ ω = ∅ is assumed. We can easily check that div(

χ ω ∂ d ψ(x, ∇u(x))) ∈ H -1 (Ω; R 3 ) if ∂ d ψ(x, ∇u(x)) ∈ C 1,α (Ω; R 3×3 ).
Then by the Lax-Milgram theorem p ∈ W 0 (Ω). To achieve the required C 2,α local regularity of p at z ∈ Ω, we only need to select z / ∈ ∂ω. Indeed, if z ∈ (Ω \ ω), there is a smooth neighborhood D ⊂ Ω \ ω of z such that div(C : ε(p)) = 0 in D. Therefore (e.g. [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Theorem 4.16]), there exists a neighborhood D D of z where p ∈ C 2,α (D ; R 3 ).

Otherwise if z ∈ ω, there exists a smooth neighborhood D ω of z where p solves the problem -div(Cε(q)) = -div(∂ d ψ(x, ∇u)) in D, q = p on ∂D.

The adjoint state p can be decomposed as p = p 1 + p 2 , where p 1 solves the problem -div(Cε(q)) = 0 in D, q = p on ∂D and p 2 solves -div(Cε(q)) = -div(∂ d ψ(x, ∇u)) in D, q = 0 on ∂D.

Therefore, there exists a neighborhood D D of z where p 1 ∈ C 2,α (D ; R 3 ). Moreover

p 2 ∈ C 2,α (D; R 3 ) if ∂ d ψ(x, ∇u(x)) ∈ C 1,α (D; R 3×3
), thanks to the interior regularity of p 2 . Hence p ∈ C 2,α (D ; R 3 ).

Particular cases

A few particular instances of the general result given by Theorem 1, and their connections to previously-available results, are now discussed.

Displacement-based functional. In this case, ∂ d ψ = ∂ d ψ = 0, and hence G = G = 0. Only the first two terms in the expression (28) of DJ then remain; moreover the second term in the right-hand side of the adjoint problem [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF] vanishes. As a result, formula (28) reduces to known results for displacement-based functionals, e.g. [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF]Prop. 4.2] if ∆ψ(z, u(z)) = 0.

Quadratic stress-based functional. This case is such that

∂ 2 dd ψ = D(z), ∂ 2 dd ψ = D (z)
, where D and D are symmetric fourth-order tensor fields. It is studied in [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF], where D and D are constant, and otherwise arbitrary, tensors. Expression [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF] of DJ for this case is indeed found, after adjusting for notational differences, to coincide with [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF]Theorem 3.1]. Reference [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF] also gives a number of useful explicit formulas on Eshelby's solution and its use in evaluating DJ for quadratic stress-based functional, including one for the evaluation of the last two integrals in [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF] when B is the unit sphere and the tensor D is isotropic.

Quadratic energy-like functional. This case, examined in [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF], is a particular instance of the previous case with D = C and D = C , and thus corresponds to functionals of the form

J a (∇u) = u -u 0 , u -u 0 Ca Ω (u 0 ∈ H 1 (Ω; R 3 ) given).
Apply our main result [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF] to this case yields

DJ(z) = |B|ε[u -u 0 ](z) : ∆C : ε[u -u 0 ](z) -ε[p](z) : A : ε[u](z) + |B|ε[u -2u 0 ](z) : ∆C : B ε[v B ] dV = |B|ε[u 0 ](z) : ∆C : ε[u 0 ](z) + ε[u -2u 0 -p](z) : A : ε[u](z) (41)
where the last expression matches the one established by another approach in [START_REF] Bonnet | The topological derivative in anisotropic elasticity[END_REF] (wherein the adjoint solution q is linked to p through q = u 0 -u + p). We have used identities

R 3 ε[v B ] : C a : ε[v B ] dV = -|B|ε[u](z) : ∆C : B ε[v B ] dV
(obtained from [START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF] with w = v B ) in the first equality of (41), and

|B|ε[u -2u 0 ] : ∆C : B ε[v B ] dV = ε[u -2u 0 ] : A -|B|∆C : ε[u](z)
(resulting from left multiplication of ( 15) by ε[u -2u 0 ]) in the second equality of (41).

Drucker-Prager penalty functional. The penalty functional considered in [START_REF] Amstutz | Topological derivative-based topology optimization of structures subject to Drucker-Prager stress constraints[END_REF] is based on the following assumptions: (a) two-dimensional isotropic elasticity, plane-strain conditions, identical Poisson ratio in background and inhomogeneity materials, (b) circular trial inhomogeneities (i.e. B taken as the unit disk), (c) ψ(d) = Ψ n (α 2 DP (C : d)), with α 2 DP denoting the yield function associated with the Drucker-Prager criterion, and a similar definition for ψ in terms of the inhomogeneity material. Using these assumptions in equation ( 28), we recover equation ( 44) of [START_REF] Amstutz | Topological derivative-based topology optimization of structures subject to Drucker-Prager stress constraints[END_REF] as a special case.

Numerical evaluation of the topological derivative

The evaluation of the topological derivative (28) requires numerical procedures, even in the simplest cases (isotropic elasticity, spherical shape for B), in particular because of the integral over the unbounded region R 3 \B. The details of such procedure depend upon whether B is an ellipsoid, or has some other shape, as ∇v B is constant inside B in the former case but not necessarily in the latter [START_REF] Eshelby | The determination of the elastic field of an ellipsoidal inclusion and related problems[END_REF]. In the sequel, we concentrate on the ellipsoidal case, which is sufficient for most applications.

3D case, ellipsoidal trial inclusion

Let B denote an ellipsoid with principal axes a 1 , a 2 , a 3 , as in Section 2.2. In that case, ∇v B ( x) inside B is given by [START_REF] Bellis | Acoustic inverse scattering using topological derivative of far-field measurements-based L 2 cost functionals[END_REF], and the two integrals over B appearing in [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF] are easily evaluated since their densities are constant.

The integral over R 3 \ B in (28) requires a more elaborate procedure, as the integration domain is unbounded and v B is spatially-dependent (and in fact decays at infinity). First x is transformed to y according to

y = (x 1 /a 1 , x2 /a 2 , x3 /a 3 ) (|y| ≥ 1). ( 42
)
It is then natural to set y = |y| ŷ, with ŷ ∈ S 2 (where S 2 is the unit sphere) and 1 ≤ |y| ≤+∞.

Finally the transformation t = |y| -1 (with 0 < t ≤ 1) is applied to the radial variable, so that (42) represents any x ∈ R 3 \ B in the form x( ŷ, t). This results in

R 3 \B G(z, ∇v B ( x)) dV ( x) = 1 a 1 a 2 a 3 1 0 1 t 4 S 2 G(z, ∇v B ( x( ŷ, t))) dS( ŷ) dt, (43)
which is to be evaluated using numerical quadrature (using e.g. a Gaussian rule for t and a Lebedev rule on S 2 [START_REF] Lebedev | A quadrature formula for the sphere of the 131st algebraic order of accuracy[END_REF]).

The exterior Eshelby tensor and its numerical evaluation. Equation ( 43) in turn requires the evaluation of ∇v B ( x( ŷ, t)) for any given quadrature point ( ŷ, t). Having assumed an ellipsoidal shape for B, ∇v B is given outside of B by (Eq. (18.6) of [START_REF] Mura | Micromechanics of Defects in Solids[END_REF])

∇v B ( x) = S ext ( x) : ε , ( 44
)
where ε is again defined by ( 12) and S ext ( x) is the exterior (spatially varying) Eshelby tensor, given by where ξ = ξ(ŝ) and N (ξ) are defined as in [START_REF] Beretta | Small volume asymptotics for anisotropic elastic inclusions[END_REF], S ( ŷ, t) is the portion of S 2 defined by S ( ŷ, t) := {ŝ ∈ S 2 : 0 ≤ ŝ• ŷ ≤ t}, L + ( ŷ, t) the circular contour on S 2 defined by L + ( ŷ, t) := {ŝ ∈ S 2 : ŝ• ŷ = t} (Fig. 1), while dS(ŝ) and dφ(ŝ) respectively denote the solid angle differential on S 2 and the polar angle differential on the circle L + ( ŷ, t). We note that definition (44) of S ext ( x) again differs somewhat from the symmetrized version used in e.g. [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF].

S ext ijmn ( x( ŷ, t)) = 1 2π C k mn S ( ŷ,t) N ik (ξ)ξ ξ j dS(ŝ) - L + ( ŷ,t) tN ik (ξ)ξ ξ j dφ(ŝ) , (45)
It is useful to recast (45) in a form more suitable for numerical To this end, let S ( ŷ, t) be represented in terms of coordinates (z, φ) spanning the fixed domain

Q := (z, φ) ∈ [0, 1] × [0, 2π[ by S ( ŷ, t) = ŝ ∈ S 2 | ŝ = (1 -t 2 z 2 ) 1/2 ˆ (φ) + tz ŷ
with ˆ (φ) spanning the unit circle S 1 := {ŝ• ŷ = 0} (this representation stems from parameterizing ŝ using angular spherical coordinates (θ, φ) and setting z := t -1 cos θ). L + ( ŷ, t) is then the subset of S ( ŷ, t) such that z = 1. The (z, φ) parametrization implies that dS(ŝ) = t dz dφ, and also induces a corresponding representation ξ(tz, φ) of ξ as defined by [START_REF] Bhargava | Elliptic inclusion in orthotropic medium[END_REF]. Inserting it in [START_REF] Tokmashev | Experimental validation of the topological sensitivity approach to elastic-wave imaging[END_REF] and rearranging the resulting expression, the exterior Eshelby tensor S ext is given, for a given evaluation point x( ŷ, t), by an integral over the fixed domain Q:

S ext ijmn ( x( ŷ, t)) = t 2π C k mn 1 0 dz 2π 0 Σ k ij (tz, φ) -Σ k ij (t, φ) dφ (46) with Σ k ij (tz, φ) = N ik (ξ(tz, φ))ξ (tz, φ)ξ j (tz, φ).
For given x( ŷ, t), S ext ( x( ŷ, t)) can then be evaluated by means of standard numerical quadrature in Q (using e.g. product rules that are Gaussian in z and uniform in φ). In addition, a Taylor expansion in

t about t = 0 yields Σ k ij (tz, φ) -Σ k ij (t, φ) = (z -1)∂ 1 Σ k ij (0, φ) + O(t 2 ),
where ∂ 1 denotes the partial derivative w.r.t. the first argument. It is moreover straightforward to show that ∂ 1 Σ k ij (0, φ) is for any (k ij) a polynomial in (cos φ, sin φ) involving only odd-degree terms, and hence that its integral over φ ∈ [0, 2π] vanishes. Consequently, the Taylor expansion of (46) about t = 0 is S ext ijmn ( x( ŷ, t)) = O(t 3 ), implying that the integral in t of ( 43) is well defined. This remark is consistent with the otherwise known O(| x| -3 ) behavior of ∇v B ( x) as | x| → +∞.

Remark 5. For a general functional with non-quadratic dependence on stress (i.e. such that ∂ 2 d dψ and ∂ 2 d dψ depend on d), the computation of (43) by numerical quadrature must be done anew for each evaluation point z of DJ(z). This makes the computation of the topological derivative field DJ potentially expensive.

On the other hand, quadratic stress-based functionals entail much less computational work, since one has (with D as in Theorem 1)

R 3 \B G(z, ∇v B ( x)) dV ( x) = ε (z) : 1 2 R 3 \B S ext ( x) : D(z) : S ext ( x) dV ( x) : ε (z),
allowing the computation of the whole field DJ(z) using just one numerical quadrature in the ( ŷ, t) variables (if D is independent of z) or up to 45 such quadratures (for general D(z), for which the major symmetry D ijk = D k ij necessarily holds).

Isotropic background material. Letting the material be in that case characterized using its shear modulus G and Poisson's ratio ν, the Christoffel tensor N has a simple closed-form expression

N ik (ξ) = 1 G|ξ| 2 δ ik - 1 2(1 -ν) ξi ⊗ ξk ( 47
)
with ξ := ξ/|ξ|. A straightforward calculation then shows that C k mn Σ k ij = C k mn Σ k ji , implying that S ext ( x) : ε is symmetric, i.e. coincides with the value obtained from the definition of the exterior Eshelby tensor used in [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF]. The exterior Eshelby tensor for a unit ball is then found (by analytical evaluation of ( 46) using ( 47), or by completing derivations presented in [START_REF] Mura | Micromechanics of Defects in Solids[END_REF]Sec. 11]) to have the following closed-form expression (wherein x := | x|):

2(1 -ν)S ext ijk ( x) = 7 x9 - 5 x7 x i x j x k x + 1 x5 - 1 x7 δ ij x k x + 1 -2ν x5 - 1 x7 δ k x i x j + ν x5 - 1 x7 (δ ik x j x + δ jk x i x + δ i x j x k + δ j x i x k ) + 1 5x 5 - 1 -2ν 3x 3 δ ij δ k + 1 5x 5 + 1 -2ν 3x 3 (δ ik δ j + δ jk δ i ). ( 48
)
Moreover, S ext ( x) : E evaluated for an arbitrary tensor E ∈ R 3×3 sym using the above formula coincides with equation ( 23) of [START_REF] Schneider | The topological gradient in anisotropic elasticity with an eye towards lightweight design[END_REF].

2D plane strain case, elliptical trial inclusion.

The exterior Eshelby tensor for an elliptic inclusion

B = x ∈ R 2 , (x 1 /a 1 ) 2 + (x 2 /a 2 ) 2 ≤ 1 , a 1 , a 2 > 0,
embedded in a 2D infinite anisotropic medium, under plane strain conditions, can also be expressed (following a direct derivation given in Appendix A) as an integral in a form suitable for numerical quadrature.

S ext ( x) = 1 π π/2 -π/2 α(θ) ⊗ N α(θ) ⊗ α(θ) dθ : C - 2 π π/2 -π/2 α(θ(w)) ⊗ N α(θ(w)) ⊗ α(θ(w)) dw : C (49)
having set x = y a 1 cos γ, a 2 sin γ , α(θ) = a -1 1 cos(θ+γ), a -1 2 sin(θ+γ) , while the function θ(w) in the second integral is defined implicitly by sin θ = 1 -y -2 sin w.

Isotropic background material, circular inclusion. The exterior Eshelby tensor for a circular inclusion of unit radius is then found, for example by analytical evaluation of (49) using [START_REF] Yang | Inhomogeneous material property assignment and orientation definition of transverse isotropy of femur[END_REF], to have the following closed-form expression (with all indices ranging in {1, 2}):

2(1 -ν)S ext ijk ( x) = 6 x8 - 4 x6 x i x j x k x + 1 x4 - 1 x6 δ ij x k x + 1 -2ν x4 - 1 x6 δ k x i x j + ν x4 - 1 x6 (δ ik x j x + δ jk x i x + δ i x j x k + δ j x i x k ) + 1 4x 4 - 1 -2ν 2x 2 δ ij δ k + 1 4x 4 + 1 -2ν 2x 2 (δ ik δ j + δ jk δ i ), ( 50
)
5 Numerical examples

The finite element analysis for each of the following 2D and 3D test cases was carried out with the software Freefem++ [START_REF]Freefem++[END_REF]. The finite elements for the displacement and adjoint state where chosen as Lagrange P 1 elements on a triangular and tetrahedral mesh, respectively. The mesh of the surface of a human femur was obtained from the mesh database of the GAMMA project [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF], while its inner tetrahedral mesh and the surface mesh adaptation were generated thanks to TetGen [43] and FreeYams [START_REF]Freeyams[END_REF], respectively. The plot of the 2D, 3D functions and meshes was done with Medit [START_REF] Medit | [END_REF].

The term in the topological derivative [START_REF]INRIA. 3D research meshes database, Gamma project[END_REF] involving an integral on R 3 \ B was evaluated using (43) and two quadrature rules, namely a 4-point Gauss-Legendre rule on t ∈ [0, 1] and a 26-point Lebedev rule [START_REF] Lebedev | A quadrature formula for the sphere of the 131st algebraic order of accuracy[END_REF] on ŷ ∈ S 2 (Sec. 5.2), while two Gauss-Legendre quadratures were used for the 2D example of Sec. 5.1. The known analytical expression [START_REF] Yuan | Application of topological sensitivity toward tissue elasticity imaging using magnetic resonance data[END_REF] of the 3D isotropic exterior Eshelby tensor S ext for the unit sphere was directly applied (Secs. 5.1 and 5.2) to avoid additional numerical quadrature work. The examples of Sec. 5.3 rely instead on the special cases ( 40) and ( 41) of DJ, avoiding both the numerical quadrature in (43) and the recourse to S ext .

Since the topological derivative depends on the derivatives of u and p, the finite element representation of DJ(z) is P 0 (piecewise constant). To facilitate graphical post-processing, a regularized version DJ r of DJ was computed, by applying a standard regularization procedure consisting in solving the variational problem

Ω DJ r w + ∇DJ r •∇w dV = Ω DJ w dV, ∀w ∈ H 1 (Ω).
The parameter controls the diffusion and regularization of DJ, at a slight expense of accuracy. For this study was set to = 10 -6 .

Sensitivity analysis of the Hill-Tsai failure criterion

Composite materials are quite popular in industry thanks to their low weight, high fatigue resistance and good endurance against corrosion. The elastic constitutive relation for such materials, restricted for two-dimensional problems to the in-plane components of the stress tensor, may be conveniently written in the principal orthotropy directions, and using the Voigt notation, in the form

   σ σ t τ t    = 1 1 -ν t ν t   E ν t E 0 ν t E t E t 0 0 0 2G t      ε ε t ε t   
where the stresses are denoted as σ in the fiber direction, σ t in the direction transverse to the fibers, and τ t for the shear stresses. This model involves four independent elastic constants, since one must have

ν t E = ν t E t .
The Hill-Tsai failure criterion, suitable for orthotropic materials, is an adaptation of the von Mises yield criterion whereby the material strength depends on the direction, according to the orientation of the fiber reinforcement. This criterion can be written as [START_REF] Gay | Composite materials: design and applications[END_REF]

α 2 (σ) := σ 2 σ2 + σ 2 t σ2 t - σ σ t σ2 + τ 2 t τ 2 t < 1, (51)
where σ , σt and τ t denote known rupture strengths. Unlike their metal counterparts, composite structures are incapable to conduct away the extreme electrical currents and electromagnetic fields generated by lighting strikes. Hence the need for protection of composite structures has prompted the development of specialized lighting strike protection materials [START_REF] Gardner | Lightning strike protection for composite structures[END_REF]. An example of such material features metallic parts of conductive material added into the laminate. This numerical example examines, by means of the topological derivative, the sensitivity of the Hill-Tsai criterion (51) to the addition of small metallic inhomogeneities to the orthotropic material, so as to determine (and avoid) those locations for which this material addition make the laminate most vulnerable to failure.

Let the domain Ω = {(x, y) ∈ (0, 2) × (0, 1)} ⊂ R 2 be occupied by a rectangular carbon/epoxy fiber ply (i.e. a composite membrane, with fiber direction e = (e x + e y )/ √ 2 and 60% fiber volume fraction), clamped to a circular support. An uniform horizontal tensile traction g = 10 -5 E e x is applied on its right side (Fig. 2).

The elastic parameters of the composite ply are E = 135 GPa, E t = 10 GPa, G t = 5 GPa, ν t = 0.3, while the ultimate tensile failure strengths involved in the criterion (51) are σ = 1500 MPa, σt = 50 MPa, τ t = 70 MPa. The metallic inclusions are considered circular and made of aluminum, whose isotropic elastic characteristics are E = 72 GPa, ν = 0.34, while the von Mises yield strength is σ = 20MPa.

Let the densities ψ, ψ entering the definition (23) of the cost functional J(C a ) be given, in terms of the penalization function Ψ n introduced in (26), by

ψ(d) = Ψ n α 2 (C : d) , ψ (d) = Ψ n α 2 (C : d) (52)
(with n ≥ 1). The function α for the composite membrane is given by (51). The corresponding function α for the aluminum is also taken of the form (51) with σ = σt = σ and τ t = σ/ √ 3, as this choice reduces the Hill-Tsai criterion to the plane-stress von Mises criterion for isotropic materials. With this choice, J(C a ) is always nonnegative; moreover, in the limit n → ∞, J(C a ) = 0 unless the threshold ( 51) is violated at some location. The value n = 5 is used in the numerical experiments to follow. The topological derivative DJ is given by ( 28), with

G(z, d) = ψ(∇u(z) + d) -ψ(∇u(z)) -∂ d ψ(∇u(z)) : d, (53)
and with G similarly defined in terms of ψ . The derivative ∂ d ψ is found to be given by

∂ d ψ(∇u) : d = 1 + α 2n (C : ∇u) (1-n)/n α 2(n-1) (C : ∇u) ∂ d α 2 (C : ∇u) : d (54)
where, since α 2 (σ) is a symmetric quadratic form in σ,

∂ d α 2 (C : ∇u) : d = 1 2 α 2 C : (∇u(z) + d ) -1 2 α 2 C : (∇u(z) -d ) .
Numerical assessment of the topological derivative evaluation. To provide a validation of expression (28) of DJ, the Hill-Tsai cost functional J(C a ) defined by ( 23) and ( 52) is numerically evaluated (with n =1, 2 or 4) for a circular inclusion of finite radius a, centered at (x, y) = (1, 2/3). The mesh (Fig. is plotted against a in Fig. 4. A numerical test of correctness of the evaluation of DJ(z) then consists in checking that e(a) = o(a) for small a. This desired trend is achieved in Fig. 4.

Topological derivative distribution in the composite membrane. Figure 5 shows the distributions of values of the Hill-Tsai criterion (51) and its topological derivative. The metallic inclusion should not be placed in zones where DJ takes higher values.

Sensitivity of the von Mises criterion for an isotropic 3D beam

Consider a beam occupying the domain Ω = { (x, y, z) ∈ (-0.1, 0.1)×(-0.5, 0.5)×(-0.1, 0.1) }, clamped on its rear face y = -0.5, made of an isotropic elastic material (Fig. 6). A traction g = xe z is applied on the front face, so as to produce torsion around the main axis x = z = 0. The remaining faces are traction-free. The beam is meshed with ∼ 75, 500 tetrahedral elements and the isotropic elastic properties of the material are normalized and given as Young's modulus E = 1 and Poisson's ratio

ν = 0.3.
We study the variation of the von Mises criterion when a small spherical cavity is introduced in the beam. This is a typical example of sensitivity analysis for topology optimization. The von Mises yield criterion sets a threshold on the equivalent stress σ eq : α(σ) ≤ 1, with α(σ) := σ eq /σ, σ eq := 3 2 dev(σ) : dev(σ)

1/2 , (55)
where dev(σ) := σ -1 3 tr(σ)I is the deviatoric stress tensor and σ is a critical stress threshold for the material, here chosen as σ = 0.1. Under the given torsional loading, the stress state then satisfies (55) in the entire beam (see Figure 7). where Ψ n is again the penalty function [START_REF] Hwu | Two-dimensional problems of the anisotropic elastic solid with an elliptic inclusion[END_REF], with n ≥ 1 and α is defined by (55). Then if n → ∞, the value of the integral of (56) on Ω is nonzero only if (55) is violated in some part of The distribution of the values of DJ for three different cut planes is plotted in Figure 8. The generation of these three planar isosurface maps entailed, for each evaluation point z, the evaluation of (43) via numerical quadrature (see the opening part of this section) and in practice required about 1 mn of CPU time on a desktop PC. This shows the practical feasibility of generating topological derivative maps even for the non-quadratic case and under 3D conditions.

Remark 6. The similarities between the distribution of |σ dev | in Figure 7 and the topological derivative in Figure 8 are noticeable. This observation supports the idea of using |σ| as a sensitivity measure for topology optimization. In fact this property is exploited by evolutionary algorithms [START_REF] Hinton | Fully stressed topological design of structures using an evolutionary procedure[END_REF] and soft kill option algorithms [START_REF] Baumgartner | SKO (soft kill option): the biological way to find an optimum structure topology[END_REF] for lightweight design subjected to a yield criterion. In general, these algorithms search the optimal topology through a fully stressed design, by gradually removing the low stressed material w.r.t. a certain reference value.

3D anisotropic non-destructive testing

Two applications of non-destructive control in anisotropic materials are now presented, one pertaining to medical imaging and the other to composite structures. Indeed, several groups have recently investigated the topological derivative as a means for imaging hidden flaws, see for instance [START_REF] Ammari | Localization, stability, and resolution of topological derivative based imaging functionals in elasticity[END_REF][START_REF] Yuan | Application of topological sensitivity toward tissue elasticity imaging using magnetic resonance data[END_REF], and also [START_REF] Tokmashev | Experimental validation of the topological sensitivity approach to elastic-wave imaging[END_REF] for an application on experimental data. Existing investigations in this direction are based on usual displacement-based cost functionals (typically of the output least-squares type). In contrast, we examine in this example an alternative approach where the misfit to experimental data is formulated in terms of an energy (and hence strainor stress-based) cost functional. Assuming the availability of a displacement measurement u 0 over a part ω of the elastic solid Ω, using e.g. full-field kinematical data [START_REF] Bay | Digital volume correlation: threedimensional strain mapping using X-ray tomography[END_REF], the strain energy of the measurement misfit is given by

J ω (C a ) = J a (∇u a ; ω) = ω ∇(u a -u 0 ) : C a : ∇(u a -u 0 ) dV. ( 57
)
Depending on whether the trial inhomogeneity B a is located inside or outside of ω, we distinguish two cases:

1. B a ⊂ ω. This case is relevant when complete displacement measurements are available over a particular region ω inside Ω. Adapting [START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF], the topological derivative reads

DJ(z) = |B|∇u 0 (z) : ∆C : ∇u 0 (z) -∇(p + 2u 0 -u) : A : ∇u (z ∈ ω),
where the adjoint solution p satisfies the variational formulation p, w

C Ω = 2 u -u 0 , w C ω , ∀w ∈ W 0 (Ω).
2. B a ⊂ Ω\ ω. The interest of this situation is justified when we have a set of small control volumes inside a body and we want to identify the position of the anomaly outside the measured volumes. As previously seen in Section 3.2, there is no second order contribution in the topological derivative and DJ simply reads

DJ(z) = -∇p(z) : A : ∇u(z) (z ∈ Ω \ ω).
First application. It is concerned with the detection of anomalous femoral bone tissue. The local change of elastic properties in femoral bone may be provoked e.g. by cancer metastasis, traumatic or pathological fractures. Moreover, bone cell elasticity and morphology changes during the cell cycle [START_REF] Kelly | Bone cell elasticity and morphology changes during the cell cycle[END_REF], and elasticity differences between cancerous and healthy tissues of various kinds have been experimentally established [START_REF] Krouskop | Elastic moduli of breast and prostate tissues under compression[END_REF][START_REF] Suresh | Nanomedicine: elastic clues in cancer detection[END_REF]. Additionally, bone is a complex material, with a multiphasic, heterogeneous and anisotropic microstructure [START_REF] Doblaré | Modelling bone tissue fracture and healing: a review[END_REF]. In particular, femoral bone can be accurately modelled a transversely isotropic material whose principal orientations are defined based either on the trabecular structures or the harvesian system, according to whether the bone is cancellous or cortical [START_REF] Yang | Inhomogeneous material property assignment and orientation definition of transverse isotropy of femur[END_REF].

Consider the proximal part Ω of a femoral bone (Fig. 9), contained in the box (x, y, z) ∈ (0.03, 0.09) × (0.04, 0.08) × (0.01, 0.11) and meshed with ∼ 213, 600 tetrahedral elements (mesh size h = 0.001). The elastic properties of the healthy bone are by simplification assumed to be homogeneous and transversely isotropic, with normalized elastic constants given by

E x = E y = 0.5, E z = 1, ν xy = ν xz = ν yz = 0.35, G xz = G yz = 0.03.
The anomalous tissue is assumed to be a small and stiffer spherical inhomogeneity (Fig. 10) with radius 0.005 and center at (0.05, 0.06, 0.06). Its Young and shear moduli are twice those of the healthy bone, while both materials have the same Poisson ratios. The bone is clamped at z = 0.01 and loaded with a vertical force density g = -1 × 10 -3 e z applied on the top surface of the femur (z ∈ [0.10, 0.11]). The measurement region ω is defined as the vertically central zone of the femur (shown in red in Fig. 10). Simulated data is assumed to be exact for simplicity. While this constitutes a strong idealization, previous numerical experiments on flaw identification by topological derivative have shown the approach to be only moderately sensitive to measurement noise [START_REF] Bellis | A FEM-based topological sensitivity approach for fast qualitative identification of buried cavities from elastodynamic overdetermined boundary data[END_REF].

Figure 11 shows three iso-surfaces of DJ with decreasing levels η ≈ -0.19 (yellow), η ≈ -0.45 (green) and η ≈ -0.71 (blue), where an iso-surface S η with level η relative to the (negative) absolute minimum DJ min := min z∈Ω DJ(z) is defined by

S η = z ∈ Ω, DJ(z) = η|DJ min | (58)
Those iso-surfaces show the location of the anomalous tissue to be correctly identified. plies, each of them composed of a weak matrix (most often polymeric) and reinforcement fibers (carbon, glass, kevlar, etc). The main failure modes in composites are fiber rupture, matrix rupture and delamination [START_REF] Gay | Composite materials: design and applications[END_REF]. We consider here a composite cube Ω = {(x, y, z) ∈ (0, 0.2) 3 }, filled up with ∼ 23, 000 tetrahedral elements (mesh size h = 0.011), made of three stacked layers of equal thickness (Fig. 12) whose constitutive elastic properties are transversely isotropic. The normalized elastic constants for the bottom layer are x, y axes, respectively. The failure point is modeled as a spherical inhomogeneity, of centre (0.1, 0.15, 0.1) and radius 0.01, and with very low elastic moduli properties C = 10 -5 C. The considered misfit criterion is again of the form (57), this time with the measurement region ω consisting of a set of M small disconnected control volumes: ω = ∪ω j (1 ≤ j ≤ M ), with ω j ⊂ Ω and ω i ∩ ω j = ∅ (two such configurations, with M = 44 and M = 729, are shown on Fig. 13). The adjoint solution p in this case solves

Ω ∇p : C : ∇q dV = 2 M i=1 ω i ∇(u -u 0 ) : C : ∇q dV, ∀q ∈ W 0 (Ω).
Fig. 14 shows one iso-surface (in green) of DJ with negative level surrounding the minimum of DJ(z) for each control volume configuration of Figure 13 (with η ≈ -0.65 and η ≈ -0.3, respectively). As expected, the identification quality improves with the number of measurement zones. Moreover, the absolute minima of DJ are found to be DJ min ≈ -1.3×10 -3 and DJ min ≈ -6.9 × 10 -2 , respectively, which is another indication of configuration 2 being more sensitive to a small defect.

Conclusion

The topological derivative DJ of a general cost functional J depending on both the displacement and its gradient (i.e. on the stress, assuming linear elastic constitutive properties) has been established by means of a rigorous small-inhomogeneity asymptotic expansion, under arbitrary 3D anisotropic conditions. This generalized version of DJ combines the previouslyknown terms associated with the topological derivative of displacement-based functionals and new terms which arise only for the case of stress-dependent functionals. Furthermore, our result holds for the case of a trial inclusions B a appearing either inside or outside the support ω ⊂ Ω of the volume density ψ a of the functional J (i.e. only the more-complex case case where B a sits on ∂ω is left out). From a computational standpoint, one of the new terms in DJ entails carrying out a numerical integration on the unbounded region R 3 \ B, and a suitable change of variables allowing subsequent use of standard quadrature formulas was given for both the 3D and 2D cases. Moreover, the computational procedure involves both the exterior and interior Eshelby tensors, for which explicit formulas are available in a few cases, and which can otherwise be evaluated by numerical quadrature of integral representation formulas such as [START_REF] Willis | Anisotropic elastic inclusion problems[END_REF]. Numerical examples demonstrate both the feasibility and usefulness of computing the field z → DJ(z), even for the case of non-quadratic stress-based functionals for three-dimensional configurations for which substantial numerical quadrature work is necessary.

A The 2D plane-strain case

Inclusion and inhomogeneity problems in plane strain have been addressed in many references, with solution methods usually based on complex potentials, see e.g. [START_REF] Bhargava | Elliptic inclusion in orthotropic medium[END_REF][START_REF] Hwu | Two-dimensional problems of the anisotropic elastic solid with an elliptic inclusion[END_REF][START_REF] Theocaris | The inclusion problem in plane elasticity[END_REF][START_REF] Willis | Anisotropic elastic inclusion problems[END_REF]. Here, we derive the exterior Eshelby tensor S ext (x) such that the displacement u(x) produced in an infinite two-dimensional anisotropic elastic medium by the application of an eigenstrain E ∈ R 2×2 sym in an elliptical region B satisfies ∇v B [E ](x) = S ext (x) :

E x ∈ R 2 . ( 59
)
The derivation is done by means of a direct evaluation of the integral representation formula of u(x), which under the present conditions (anisotropic, plane strain) reads The above fourth-dimensional integral over (ξ, η) ∈ B × R 2 is now evaluated with the help of coordinate transformations. To this aim, B is assumed without loss of generality to be centered at the coordinate origin and with its principal directions directed along the coordinate axes. The evaluation point x ∈ R 2 is parameterized as

u(x) = i (2π) 2
x(y, γ) = y (a 1 cos(γ), a 2 sin(γ)) , γ ∈ [0, 2π], y ∈ (0, ∞),

and two mappings are introduced. First, f : (t, θ) ∈ R + \{0}×[0, 2π[ -→ (η 1 , η 2 ) ∈ R 2 \{(0, 0)} is defined by η 1 (t, θ) = tα 1 (θ) η 2 (t, θ) = tα 2 (θ) with α(θ) = (α 1 , α 2 )(θ) := a -1 1 cos(θ + γ), a -1 2 sin(θ + γ) , (62) which implies dV (η) = (a 1 a 2 ) -1 t dt dθ, Then, for given η ∈ R 2 \ {0} (i.e. for given (t, θ) ∈ R + \ {0}), g : (z 1 , z 2 ) ∈ D -→ (ξ 1 , ξ 2 ) ∈ B (where D ⊂ R 2 is the closed unit disk) is defined by

ξ 1 = a 1 ( z 1 sin θ + z 2 cos θ) ξ 2 = a 2 (-z 1 cos θ + z 2 sin θ) with -1 -z 2 2 ≤ z 1 ≤ 1 -z 2 2 , -1 ≤ z 2 ≤ 1 (63)
which implies dV (ξ) = a 1 a 2 dz 1 dz 2

Now, mappings (62) and ( 63) are substituted into the integral representation formula (60). Noting that definitions (61), ( 62) and (63) imply

α(θ)•x = y cos θ, η•(ξ -x) = t[z 2 -α(θ)•x],
one first finds, using mapping (63) and the fact that the integrand of the resulting integral over (z 1 , z 2 ) does not depend on z 1 , that B e iη•(ξ-x) dV (ξ) = a 1 a 2 e -itα(θ)•x 1 -1 2 1 -z 2 2 e itz 2 dz 2 = a 1 a 2 2π t e -itα(θ)•x J 1 (t),

where J 1 is the Bessel function of first kind and order 1, the last equality stemming from formula 3.752 of [START_REF] Gradshteyn | Tables of integrals, series and products[END_REF] together with the function z 2 -→ 1 -z 2 2 sin kz 2 being odd. Then, using the above result into (60) and applying mapping (62), one obtains i (2π) 2 The last equality above is established by using that (i) α(θ) is 2π-periodic, (ii) α(θ + π) = -α(θ) and (iii) N (-η) = N (η) (N being homogeneous of degree -2).

The inner integral in (64) in fact admits a known closed-form expression (formula 6.693(1) of [START_REF] Gradshteyn | Tables of integrals, series and products[END_REF]), which depends on the value of α(θ)•x = y cos θ: where the integral can now be evaluated by usual quadrature rules.

∞ 0 sin t α(θ)•x J 1 (t) dt t = α(θ)•x 0 ≤ α(θ)•x ≤ 1 (65a) = α(θ)•x -(α(θ)•x) 2 -1 α(θ)•x ≥ 1 (65b) If x ∈ B, (

Figure 1 :

 1 Figure 1: Parametrization of the set S ⊂ S 2 .

Figure 2 :

 2 Figure 2: Orthotropic ply Ω clamped to a circular support: configuration and mesh (∼ 3, 700 elements). A horizontal uniform load g is applied on its right side.

Figure 3 :Figure 4 :

 34 Figure 3: Orthotropic ply: mesh for the numerical validation of DJ.

  Values of the Hill-Tsai failure criterion Values of the topological derivative

Figure 5 :

 5 Figure 5: Orthotropic ply: sensitivity analysis of the Hill-Tsai failure criterion.

Figure 6 :

 6 Figure 6: Beam under torsion: geometrical configuration, with color scale of displacement modulus (undeformed and deformed views)

1. 6 × 10 - 5 3. 2 × 10 - 3 6. 4 × 10 - 3 1. 6 ×

 6105210341036 10 -3 4.8 × 10 -3

Figure 7 :

 7 Figure 7: Beam under torsion: contour plot of 2σ 2 eq /3 = |dev(σ)| 2 on the boundary (left) and inside the beam along a transversal cut (right).

Figure 8 :

 8 Figure 8: Beam under torsion: color maps of |DJ| in planes x = 0 (top left), x = 0.05 (top right) and x = 0.1 (bottom). The left boundary corresponds to the clamped face of the beam.

Figure 9 :

 9 Figure 9: Femoral bone: a uniform vertical load g = -10 -3 e z is applied on the head of the femur (simulating the body weight) while the distal horizontal section z = 0.01 is clamped (left). The right panel shows a color scale of the displacement modulus.

Figure 10 :

 10 Figure 10: Femoral bone: anomalous tissue (left, in red) and measurement region ω (right, in red).

  Three-dimensional view. (b) Transversal cut view.

Figure 11 :

 11 Figure 11: Femural bone: iso-surfaces η ≈ -0.19 (yellow), η ≈ -0.45 (green) and η ≈ -0.71 (blue) of DJ.

E x = 1 ,

 1 E y = E z = 0.05, ν xy = ν xz = ν yz = 0.35, G xz = G yz = G xy = 0.03, The middle and top layers have the same elastic constants than the lowest one, but with the horizontal principal orthotropy directions resulting from a 45 o and 90 o rotation of the 0Layered cube, and fiber orientation.

2 (

 2 b) A color scale of the displacement modulus.

Figure 12 :

 12 Figure 12: Layered cube: two uniform compression loads g 1 = -0.1e x and g 2 = 0.1e x are respectively applied on the faces of the cube x = 0.2 and x = 0. Displacements are not affine in the spatial coordinates due to the anisotropy.

( a )

 a Configuration 1 (M = 44). (b) Configuration 2 (M = 729).

Figure 13 :

 13 Figure 13: Layered cube: control volumes ω j (1 ≤ j ≤ M ).

  Configuration 1 (M = 44).

  Configuration 2 (M = 729).

Figure 14 :

 14 Figure 14: Layered cube: iso-surface S η with η ≈ -0.65 (left) and η ≈ -0.3 (right) of DJ. The grey sphere shows the correct location of the failure point.

R 2 B

 2 exp(iη•(ξ -x))N (η) ⊗ η dV (ξ) dV (η) : C : E (60)

R 2 B

 2 exp(iη•(ξ -x))N (η) ⊗ η dV (ξ) dV (η)

= 1 π π/ 2 -π/ 2 απ π/ 2 -π/ 2 α

 2222 61) implies that y ≤ 1, and hence that α(θ) • x ≤ 1 for any θ. From (60), (64) and (65a,b), u(x) is then such thatu(x) = xS int : E , ∇u(x) = S int : Ewhere S int is the plane-strain interior Eshelby tensor, given byS int (θ) ⊗ N α(θ) ⊗ α(θ) dθ : C If x ∈ R 2 \ B, (61) implies that y > 1. Let θ = arccos(1/y), so that the subset of θ ∈ [-π/2, π/2] where α(θ)•x ≥ 1 is θ ∈ [-θ, θ].In that case, using (60), (64) and (65a,b) and differentiating the resulting expression of u(x) with respect to x, one finds ∇u(x) = S ext (x) : E with the plane-strain exterior Eshelby tensor S ext (x) given byS ext (x) = S int -2 π θ θ (α(θ)•x) (α(θ)•x) 2 -1 α(θ) ⊗ N α(θ) ⊗ α(θ) dθ : CEven though this representation of the 2D general Eshelby tensor is valid, is not suited for numerical evaluation due to the term 1/ (α • x) 2 -1, which is (weakly) singular at the endpoints θ = ± θ. We recast it into a form suitable for numerical quadrature by settingsin θ = 1 -y -2 sin w, θ ∈ [-θ, θ], w ∈ [-1, 1]Then, since α(θ)•x = y cos θ, one easily finds that(α(θ)•x) dθ = y 2 -1 cos w dw, (α(θ)•x) 2 -1 = y 2 -1 cos wConsequently, S ext (x) is now expressed as S ext (x) = S int -2 (θ(w)) ⊗ N α(θ(w)) ⊗ α(θ(w)) dw : C