Marc Bonnet
email: bonnet@lms.polytechnique.fr

Andrei Constantinescu
email: andrei.constantinescu@lms.polytechnique.fr

Inverse problems in elasticity

come

Inverse problems in elasticity

Introduction

Elasticity theory describes the reversible deformation of solid bodies subjected to excitations of various physical natures: mechanical, thermal, electromagnetical etc. Such excitations, applied as distributions over the body (e.g. gravitation, Lorentz forces, thermal expansion) or over the boundary (pressure, contact forces), generate strains (i.e. local deformations) and stresses (i.e. local forces) in the material. Elasticity is a mechanical constitutive property of materials whereby (i) a one-to-one relationship between instantaneous strains and stresses on the current deformed configuration is assumed, and (ii) the material reverts to its initial state if the sollicitation history is reversed.

Almost all natural or manufactured solid materials have a deformation range within which their mechanical behaviour can be modelled by elasticity theory. For sufficiently small strains, the elastic behaviour is considered as linear, i.e. strains and stresses are assumed to be proportional to each other. A vast body of engineering experience shows that the theory of linear elasticity allows an accurate modelling of many man-made or natural objects: civil engineering structures, transportation vehicles, machines, the earth mantle (to list just a few), and provides an essential tool for analysis and design. In addition to the basic theory for threedimensional solid media, specialized approches have been developed for cases featuring two or more dissimilar length scales: composite media, slender structures (beams, plates, shells).

The theory of linearized elasticity has developed into one of the now classical areas of mathemematical physics. Equilibrium problems are governed by elliptic partial differential equations, similar to those of electrostatics but more complex in that physical quantities of interest are described by tensor fields rather than vector fields. Closed-form solutions are available only for simple geometries (usually corresponding to separable coordinate systems), so that most real-life modelling studies are based on numerical solution methods. Dynamic conditions give rise to hyperbolic partial differential equations.

The main types of inverse problems that arise in the context of linear elasticity, and more generally of the mechanics of deformable solids, are similar to those encountered in other areas of physics involving continuous media and distributed physical quantities, e.g. acoustics, electrostatics and electromagnetism. They are usually motivated by the desire or need to overcome a lack of information concerning the properties of the system (a deformable solid body or structure). Mathematical and numerical techniques for the reconstruction of buried objects of a geometrical nature, such as cracks, cavities or inclusions, are the subject of many investigations, e.g. [3, 7, 8, 15, 26, 31, 32, 44-46, 61, 76, 95, 96, 118, 120]. Mechanical waves, such as ultrasonic or Lamb waves, are also frequently used in practical non-destructive testing of structures, see e.g. [103-105, 112, 131] and the references provided therein. The identification of distributed parameters [5, 12-14, 29, 39, 48, 49, 60, 70, 81-83, 116] (e.g. elastic moduli, mass density, wave velocity) arises in connection with e.g. medical imaging of tissues [START_REF] Barbone | Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions[END_REF] or seismic exploration [START_REF] Kurpinar | Numerical solution of the inverse problem for the elasticity system for horizontally stratified media[END_REF][START_REF] Menke | Geophysical data analysis : discrete inverse theory[END_REF][START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF][START_REF] Santos | On the solution of an inverse scattering problem in seismic whiledrilling technology[END_REF][START_REF] Tarantola | Inverse problem theory[END_REF][START_REF] Weglein | Inverse scattering series and seismic exploration[END_REF]. The reconstruction of residual stresses [START_REF] Ballard | On the inversion of subsurface residual stresses from surface stress measurements[END_REF][START_REF] Gao | On the inversion of residual stresses from surface measurements[END_REF][START_REF] Lin | Uniqueness in inverse problems for an elasticity system with residual stress by a single measurement[END_REF][START_REF] Robertson | Boundary identifiability of residual stress via the Dirichlet to Neumann map[END_REF] is a related topic with important engineering implications. Models of complex engineering structures often feature local parameters that are not known with sufficient accuracy, and therefore need to be corrected by exploiting experimental information on the mechanical response of the structure. Model updating is often treated as an inverse problem [START_REF] Ben Abdallah | Inversion gaussienne appliquée à la correction paramétrique de modèles structuraux[END_REF][START_REF] Bonnet | On the estimation of the geometrical support of modelling defects using the distributed error in constitutive equation[END_REF][START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF][START_REF] Deraemaeker | Reduced bases for model updating in structural dynamics based on constitutive relation error[END_REF][START_REF] Ladev Èze | Application of a posteriori error estimation for structural model updating[END_REF][START_REF] Link | Requirements for the structure of analytical models used for parameter identification[END_REF][START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF][START_REF] Zimmermann | Structural damage detection using a minimum rank update theory[END_REF], in particular because corrections affect distributed parameters over a limited region of space which is not known beforehand. Identification of sources or inaccessible boundary values (i.e. Cauchy problems in elasticity) are also encountered [START_REF] Cimeti Ère | Solution of the Cauchy problem using iterated Tikhonov regularization[END_REF][START_REF] Constantinescu | Identification du chargement thermique transitoire par contrôle optimal[END_REF][START_REF] Delattre | Application du contrôle optimal à l'identification d'un chargement thermique[END_REF][START_REF] Koslov | An iterative method for solving the Cauchy problem for elliptic equations[END_REF][START_REF] Maniatty | Investigation of regularization parameters and error estimating in inverse elasticity problems[END_REF][START_REF] Marin | Boundary element solution for the Cauchy problem in linear elasticity using singular value decomposition[END_REF][START_REF] Marin | BEM first-order regularisation method in linear elasticity for boundary identification[END_REF]. Finally, the identification of homogeneous constitutive properties is increasingly often made on the basis of measurements taken on structures, for which simplifying assumptions such as constant states of strain or stress are invalid, and inverse techniques are then developed for that purpose [START_REF] Constantinescu | On the identification of elastoviscoplastic constitutive laws from indentation tests[END_REF][START_REF] Forestier | 3D inverse analysis model using semianalytical differentiation for mechanical parameter identification[END_REF][START_REF] Gavrus | An inverse analysis using a finite element model for identification of rheological parameters[END_REF][START_REF] Mahnken | Parameter identification for viscoplastic models based on analytical derivatives of a least-square functionals and stability investigations[END_REF][START_REF] Tardieu | On the determination of elastic coefficients from indentation experiments[END_REF].

This article is devoted to some of the above-mentioned inverse problems, namely the identification of (i) distributions of elastic moduli, (ii) model parameters, (iii) buried cracks or other geometrical objects. These inverse problems will be considered mainly for threedimensional elastic media under equilibrium or dynamical conditions, and also for thin elastic plates. The main goal is to overview some recent results, in an effort to bridge the gap between studies of a mathematical nature and problems defined from engineering practice. Accordingly, emphasis will be given to formulations and solution techniques which are well suited to general-purpose numerical methods for solving elasticity problems on complex configurations, in particular the finite element method [START_REF] Brenner | The mathematical theory of finite element methods[END_REF][START_REF] Reddy | Introduction to the finite element method[END_REF] and the boundary element method [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Chen | Boundary element methods[END_REF]. The important role of the variational principles of elasticity, which in particular provide the foundations of the above-mentioned numerical solution methods, will be highlighted throughout this article. Additionnally, investigations of a more mathematical nature will also be reviewed.

The article is organised as follows. An overview of basic theory and equations of linear elasticity under the small strain hypothesis, including the virtual work principle and variational formulations of equilibrium problems, is presented in Section 2. Then, Section 3 describes strategies for the identification of distributions of elastic moduli or cracks exploiting the virtual work principle as an observation equation, with emphasis on the reciprocity gap concept. The virtual work principle is also an effective tool for setting up parameter sensitivity analyses and computing gradients of cost functions associated to identification problems, as shown in Section 4. Then, Section 5 is devoted to cost functions and parameter identification techniques based on the error in constitutive equation (ECE). Formulations for both threedimensional bodies and plates are presented, and the ability of the energy density function associated to the ECE to outline the geometrical support of defects is discussed. Finally, Section 6 is devoted to geometrical sensitivity tools, based on adjoint solutions, for defect identification. A concise formulation for cavity shape sensitivity is followed by more recent results concerning crack shape sensitivity and a presentation of the topological derivative associated with wave scattering in the limit of vanishingly small objects.

A brief review of the typical orders of magnitude involved in elastic solid bodies will close this introductory section. Values of elastic constitutive parameters for common isotropic materials are given in table [START_REF] Allix | A new identification strategy for dynamic problems with uncertain measurements[END_REF], where the Young modulus E and the Poisson ratio ν are defined in terms of the basic experiment performed by applying a traction force F to both extremities of a cylindrical bar. Under this experiment, the axial length stretches from -ℓ 0 to ℓ while the cross-section shrinks from S 0 to S, and one sets E = (F/S 0 )/(ℓ/ℓ 0 -1) and 1. Constitutive parameters of some common isotropic elastic materials: E Young modulus, ν Poisson ratio, ρ mass density; σ Y and ε Y define the elastic limit, i.e. are the stress and strain levels beyond which the material is no longer elastic.

ν = -( S/S 0 -1)/(ℓ/ℓ 0 -1). Linear elasticity is usually valid when strains are small (typical magnitudes are 10 -3 or less), and is also restricted to stress levels below a certain threshold σ Y beyond which irreversible constitutive properties, e.g. plasticity, set in (see Table 1).

Strains and displacements can be measured directly, e.g. using strain gages, whereas stresses can only be measured indirectly. Classical strain gages are reliable up to 10 -6 but offer only a "pointwise" measurement, in practice over an area of a few square millimeters. Modern technology based on laser interferometry or image correlation techniques [START_REF] Berthaud | Experiment/computation interactions by using digital image correlation[END_REF] are reliable up to 10 -5 but allow measurements of practically continuous fields over extended areas. Such experimental techniques yield rich experimental data and are therefore well-suited to identification problems. The importance of the latter techniques is increasing as inversion techniques specifically exploiting availability of field quantities (either on the boundary or over part of the domain itself) become accessible.

Review of governing equations

Fundamental field equations for three-dimensional elasticity

The deformation of an elastic body, occupying in its undeformed state the region Ω ⊂ R 3 bounded by the surface S, is usually described in terms of a vector displacement field u(x, t) (x ∈ Ω) which is such that the deformation process moves a small material element lying at x to its new position x + u(x, t). The linearized elasticity theory [START_REF] Gurtin | The Linear Theory of Elasticity[END_REF] is established on the assumption of small strains, namely |∇u(x, t)| ≪ 1. In that case, the changes in metric induced by the deformation are described by the linearized strain tensor ε(x, t), defined as a differential operator on u by: ε[u](x, t) = (∇u(x, t) + ∇u T (x, t))/2

(1)

This equation is often referred to as the compatibility equation for small deformations. The strain ε(x, t) is a symmetric second-order tensor. The material is characterized by two constitutive parameters: its mass density distribution ρ(x), associated with the kinetic energy

T (u) = 1 2 Ω ρ| u| 2 dV
(where the dot denotes time differentiation) and the fourth-order tensor of elastic moduli C(x), hereafter referred to as the elasticity tensor, associated with the elastic strain energy

E(u) = 1 2 Ω ε[u] : C : ε[u] dV (2)
The elasticity tensor C defines a positive definite quadratic form over the 6-dimensional space of symmetric second-order tensors. Therefore, C has the following symmetries:

C ijkℓ = C kℓij = C jikℓ (1 ≤ i, j, k, ℓ ≤ 3) (3)
and hence has at most 21 independent coefficients. In the simplest situation of isotropic elasticity, C depends on only two independent moduli. For instance, C can be expressed in terms of the Lamé coefficients (λ, µ):

C ijkℓ = λδ ij δ kℓ + µ(δ ik δ jℓ + δ jk δ iℓ ) (4)
Other commonly used elastic parameters are the Young modulus E and the Poisson ratio ν, which are related to the Lamé constants by

E = µ(3λ + 2µ) λ + µ ν = λ 2(λ + µ) (5)
The stress tensor σ describes internal forces: the traction vector

p[n](x, t) = σ(x, t).n(x) (6)
is such that p[n] dS is the elementary force applied on a infinitesimal surface patch dS of unit normal n located at x ∈ Ω. The fundamental balance equation of the dynamics of deformable bodies (an extension of Newton's second law to a small material element) is then:

div σ(x, t) + f (x, t) -ρ(x) ü(x, t) = 0 (7)
where f (x, t) is a given distribution of body forces. The constitutive assumption of linearized elasticity, adopted in this article, postulates that the stress tensor σ(x, t) depends linearly on the linearized strain tensor, i.e.:

σ(x, t) = C(x) : ε[u](x, t) (8)
On combining the three field equations ( 1), ( 7) and ( 8) and eliminating ε and σ, the displacement field is found to be governed by the partial differential equation

[A C u](x, t) + f (x, t) -ρ(x) ü(x, t) = 0 (9)
with the elasticity operator A C defined by:

A C u := div (C : ε[u]) = div (C : ∇u) (10)
(where the last equality stems from the constitutive symmetries (3)). Equation ( 9) is the analog for linear elasticity of the hyperbolic linear wave equation. Besides, a well-posed elastodynamic problem features initial conditions

u(x, 0) = u 0 (x) u(x, 0) = v 0 (x) (x in Ω) (11)
and boundary conditions on S, for instance displacements ξ and tractions φ prescribed on complementary portions S u and S p = S \ S u of S:

u(x, t) = ξ(x, t) (x on S u , t ∈ [0, T ]), (12)
p C [u](x, t) = φ(x, t) (x on S p , t ∈ [0, T ]) (13)
where the boundary traction operator u → p C [u], the counterpart for elasticity of the normal derivative operator, is defined by

p C [u] = (C : ∇u).n = (C : ε[u]).n (14)
The traction operator will be simply denoted p[u] when there is no ambiguity about C.

Direct problems

The typical direct problem of elastodynamics is the initial-boundary value problem for the unknown displacement field u(x, t) defined by the field equation ( 9), the initial conditions [START_REF] Barbone | Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions[END_REF] and the boundary conditions ( 12) and ( 13). The geometry (Ω), the physical characteristics (C, ρ) of the elastic body, the structure of boundary conditions, the prescribed values (ξ, φ) on the boundary, and the initial data u 0 , v 0 are assumed to be known. Of course, well-posed initial and boundary conditions other than ( 12), ( 13) and ( 11) could be considered as well.

The direct elastic equilibrium problem has a similar structure, except that all field variables are time-independent, so that the field equation ( 9) reduces to the analog for linear elasticity of the Laplace equation, i.e.

[A C u](x) + f (x) = 0 (15)
and initial conditions such as [START_REF] Barbone | Elastic modulus imaging: on the uniqueness and nonuniqueness of the elastography inverse problem in two dimensions[END_REF] are not needed. The direct equilibrium problem for thin plates is defined in Section 5.3.1. Finally, when time-harmonic motions are considered, i.e.:

u(x, t) = Re[u(x)e iωt ], (16)
(where for simplicity the same notation is used for space-time fields and their frequencydomain counterparts) the complex-valued unknown field u(x) solves the analog for linear elasticity of the Helmholtz equation:

[A C u](x) + ρω 2 u(x) + f (x) = 0 (17)
together with boundary conditions such as [START_REF] Barcilon | Inverse problem for vibrating beam[END_REF] and [START_REF] Barcilon | On the multiplicity of solutions of the inverse problem for a vibrating beam[END_REF], where the prescribed data ξ, φ is also complex-valued and obeys the time-harmonic convention [START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF].

The principle of virtual work

The principle of virtual work [START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF] is a fundamental principle of the dynamics of continua, from which the balance equation [START_REF] Andrieux | Sur l'identification de fissures planes via le concept d'écart à la réciprocité en élasticité[END_REF] can then be deduced. Conversely, it can be derived as a weak form of ( 7) treated as a first principle of continuum mechanics: taking the inner product of ( 7) by an arbitrary trial field w(x) (sometimes termed virtual velocity and assumed continuous in this article) leads to

- Ω σ(x, t) : ε[w](x) dV + Ω f (x, t).w(x) dV + S p[n](x, t).w(x) dS = Ω ρ(x) ü(x, t).w(x) dV (∀w) (18)
which constitutes the virtual work principle for deformable continuous media. Equation [START_REF] Ben Ameur | Level set methods for geometric inverse problems in elasticity[END_REF] is general in that it does not refer to specific constitutive properties (e.g. elasticity).

On substituting the constitutive equation ( 8) into [START_REF] Ben Ameur | Level set methods for geometric inverse problems in elasticity[END_REF], a weak formulation of the governing equation ( 9) for the displacement field u(x, t) is obtained. For instance, the following identity holds for elastic equilibrium problems in the absence of body forces:

Ω ε[u](x) : C(x) : ε[w](x) dV = S p C [u](x).w(x) dS (∀w) (19)

Variational formulations

Solutions to equilibrium problems in elasticity have well-known characterizations as minimizers of energy functionals [START_REF] Oden | Variational methods in theoretical mechanics[END_REF]. For a body endowed with a given elasticity tensor C and a well-posed set of boundary conditions of the form (12-13), the displacement field u solving the elastic equilibrium problem defined by ( 12), [START_REF] Barcilon | On the multiplicity of solutions of the inverse problem for a vibrating beam[END_REF], and ( 15) minimizes over the space of kinematically admissible displacements the potential energy

u = arg min v∈C(ξ,Su) W C (v) with W C (v) = 1 2 Ω ε[v] : C : ε[v] dV - Ω f .v dV - Sp φ.v dS (20)
while the stress field σ solution to the same equilibrium problem minimises over the space of statically admissible stress fields the complementary energy

σ = arg min s∈S(φ,Sp) W ⋆ C (s) with W ⋆ C (s) = 1 2 Ω s : C -1 : s dV - Su u.(s.n) dS (21)
In ( 20) and ( 21), the spaces of admissible fields are defined by

C(ξ, Σ) = {u : u = ξ on Σ} (22) S(φ, Σ) = {σ : div σ + f = 0 in Ω, σ.n = φ on Σ} ( 23
)
where Σ ⊆ S is a generic surface. The potential energy and the complementary energy are dual in the sense of the Fenchel duality theorem [START_REF] Chavent | Primal-dual formulations for parameter estimation problems[END_REF][START_REF] Oden | Variational methods in theoretical mechanics[END_REF]. As a direct consequence of the Fenchel duality theorem one can express the solution pair (u, σ) as a minimizer over C(ξ, S u )×S(φ, S p ) of the sum of the potential and complementary potential energy:

(u, σ) = arg min (v,s)∈C(ξ,Su)×S(φ,Sp) [W C (v) + W ⋆ C (s)] . (24)
Moreover, this minimum is equal to zero, i.e.:

W C (u) + W ⋆ C (σ) = 0. (25)

The virtual work principle as an observation equation

The distribution of displacements and in-plane strains can be measured over the surface of a body by optical means. Therefore, displacement and in-plane strain fields on the surface are experimentally available. Inversion strategies based on the assumed availability of continuous data have therefore a practical relevance. In such situations, the virtual work principle allows to formulate observation equations, i.e. mathematical relationships between the observations and the unknown quantities, in a direct and effective way. Most of this section is devoted to cases of overspecified data on the boundary, in which case observation equations are formulated from the virtual work principle in the form of a reciprocity gap. This approach is considered in connection with the identification of distributed elastic moduli (section 3.1) and cracks (section 3.3). The virtual fields method will also be briefly described for completeness (section 3.4).

The reciprocity gap for elastic moduli identification

In this section, the reciprocity gap concept is presented in connection with the identification of an unknown distribution of elastic moduli from complete non-intrusive experiments conducted in quasi-static conditions. Let us consider a body Ω endowed with unknown elastic moduli C ⋆ (x) to be identified from measurements. Known forces φ applied on S induce a displacement u ⋆ in Ω, and the trace ξ of u ⋆ on S is measured. The displacement u ⋆ and elasticity tensor C ⋆ must therefore satisfy the field equation of elastic equilibrium

div (C ⋆ : ε[u ⋆ ]) = 0 in Ω (26)
and be linked to the boundary data pair (ξ, φ) through the overdetermined boundary conditions

p C ⋆ [u ⋆ ] = φ and u ⋆ = ξ on S (27)
Naturally, several experimental pairs (ξ (ℓ) , φ (ℓ) ), each associated with a displacement field u ⋆(ℓ) , may be considered instead of just one.

Λ C (ξ) = p C [u] where u solves A C u = 0 (in Ω) , u = ξ (on S) (28)
The issue of whether the knowledge of Λ C ⋆ uniquely determines C ⋆ is not fully settled. For general anisotropic elastic media, the following argument suggests that the answer is negative.

Let y = Ψ(x) denote a diffeomorphism such that Ψ(x) = x (x ∈ S) and Ψ(Ω) = Ω.

Introducing the mapping defined by x = Ψ -1 (y) into the domain integral in [START_REF] Ïm | Identification of certain polynomial nonlinear structures by adaptive selectively-sensitive excitation[END_REF], one obtains

S [Λ C u](x).w(x) dS x = Ω ∇ x u(x) : C(x) : ∇ x w(x) dV x = Ω [∇ y (u • Ψ -1 )(y)] : L(y) : [∇ y (w • Ψ -1 )(y)] dV y = S Λ L (u • Ψ -1 )(y).w(y) dS y ( 29
)
where the symmetry properties (3) of C have been used and L(y) is defined by

L ijkℓ (Ψ(x)) = |det∇Ψ| -1 (x)C imkn (x)Ψ j,m (x)Ψ ℓ,n (x) (30)
In particular, the last equality in [START_REF] Bui | Sur quelques problèmes inverses élastiques en mécanique de l'endommagement[END_REF] holds true because the identity div L(y) : [∇ y (u • Ψ -1 )(y)] = 0 can be established by first considering the trial fields w that vanish on S. If L(y) defined by ( 30) is a legitimate elasticity tensor, then it follows from (29) that two distinct elasticity tensors C(y) and L(y) produce the same DtN map. Definition [START_REF] Bui | Spatial localization of the error of constitutive law for the identification of defects in elastic bodies[END_REF] implies that L(y) possesses the major symmetry L ijkℓ = L kℓij and that the bilinear form ε → ε : L(y) : ε is positive. However, the necessary minor symmetry L ijkℓ = L jikℓ is not assured by definition [START_REF] Bui | Spatial localization of the error of constitutive law for the identification of defects in elastic bodies[END_REF], and requires the satisfaction of

C imkn (x)Ψ j,m (x)Ψ ℓ,n (x) = C jmkn (x)Ψ i,m (x)Ψ ℓ,n (x) (1 ≤ i, j, k, ℓ ≤ 3) (31)
A careful examination of these symmetry-enforcing constraints, taking again into account the symmetry properties (3) of C, reveals that they amount to 18 independent homogeneous linear equations for 21 independent elastic moduli. Therefore, one can find fourth-order tensors C having the symmetry properties (3) such that the fourth-order tensor L linked to C by ( 30) is an elasticity tensor. Such pairs (C, L) of tensors therefore constitute potential examples of different elasticity tensors that lead to the same DtN map. For this argument to be complete, it is however necessary to prove that the null space of the set of equations ( 31) contains at least one tensor C that defines a positive definite quadratic form over the symmetric second-order tensors, which to our best knowledge is still an open question. The foregoing argument, given in [START_REF] Constantinescu | Sur l'identification des modules élastiques[END_REF], is an extension to elasticity of the line of reasoning used by Kohn and Vogelius [START_REF] Kohn | Relaxation of a variational method for impedance computed tomography[END_REF] to prove that two distinct distributions of anisotropic electric conductivities can have the same DtN map. The latter situation is simpler in that the conductivity tensors are of order two instead of four, and hence the symmetry requirement [START_REF] Bui | Diffraction acoustique inverse de fissure plane : solution explicite pour un solide borné[END_REF] is not involved. In the case of isotropic elasticity, where C depends on just two independent moduli (e.g. the Lamé coefficients (λ, µ) in the representation (4)), Nakamura and Uhlmann first published [START_REF] Nakamura | Global uniqueness for an inverse boundary problem arising in elasticity[END_REF] a theorem stating the identifiability of (λ, µ) from the knowledge of Λ λ,µ . However, a technical error was subsequently found by these authors [START_REF] Nakamura | Global uniqueness for an inverse boundary value problem arising in elasticity (erratum)[END_REF] and Eskin and Ralston [START_REF] Eskin | On the inverse boundary value for linear isotropic elasticity[END_REF]. As a result, the above identifiability result holds only in weaker forms [START_REF] Eskin | On the inverse boundary value for linear isotropic elasticity[END_REF][START_REF] Nakamura | Global uniqueness for an inverse boundary value problem arising in elasticity (erratum)[END_REF] requiring additional prior information or constraints on (λ, µ). For instance, one such result states that there exists ε > 0 and m ∈ N such that

|∇µ i | C m ( Ω) < ε and Λ λ 1 ,µ 1 = Λ λ 2 ,µ 2 implies λ 1 = λ 2 and µ 1 = µ 2 on Ω.

The reciprocity gap functional.

To define a reciprocity gap, consider two elastic bodies occupying the same region Ω ∈ R 3 and characterized by two distinct distributions of elastic moduli C ⋆ (x) and C(x), and two displacement fields u ⋆ (x) and w(x) such that [A C ⋆ u ⋆ ](x) = 0 and [A C w](x) = 0, respectively, in Ω. Then, equation [START_REF] Ïm | Identification of certain polynomial nonlinear structures by adaptive selectively-sensitive excitation[END_REF] is written (i) as the weak form of [A C ⋆ u ⋆ ](x) = 0 with w as the virtual field, and (ii) as the weak form of [A C w](x) = 0 with u ⋆ as the virtual field. Taking into account the symmetry properties (3) of C and C ⋆ and the overdetermined boundary data [START_REF] Bonnet | On the estimation of the geometrical support of modelling defects using the distributed error in constitutive equation[END_REF], one arrives at the identity

S w.φ -ξ.p C [w] dS = Ω ∇u ⋆ : (C ⋆ -C) : ∇w dV := R(C ⋆ -C; u ⋆ , w) (32)
which defines the reciprocity gap R(C ⋆ -C; u ⋆ , w). The terminology comes from the fact that the identity (32) can also be established from the Betti reciprocity theorem [START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF] applied to the two states w and u ⋆ . The Betti theorem being based on the assumption that both states refer to the same (symmetric) elasticity tensors, a "reciprocity gap" occurs, here in the form of the domain integral in [START_REF] Bui | Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap[END_REF] when that assumption is not correct. For any selection of the elastostatic virtual field w (sometimes also termed "adjoint field"), the left-hand side of ( 32) is known, while the right-hand side depends on the unknown moduli C ⋆ . In fact, since u ⋆ itself depends on C ⋆ , equation [START_REF] Bui | Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap[END_REF] is nonlinear in C ⋆ . It is therefore natural to consider its linearized counterpart, the linearized reciprocity gap.

Identification of elastic moduli perturbations using the linearized reciprocity gap.

If C ⋆ differs from C by a small perturbation, i.e. C ⋆ = C + δC, the approximate observation equation obtained from [START_REF] Bui | Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap[END_REF] by retaining the leading O(|δC|) contribution is linear in δC:

S [ξ.p C [w] -w.φ] dS = Ω ∇u : δC : ∇w dV + o(|δC|) (33)
where u solves

[A C u](x) = 0 (in Ω) , σ[u].n = φ (on S) (34)
and is therefore known beforehand. The linearized reciprocity gap (LRG) identity ( 33) is useful for either theoretical investigations of the linearized inverse problem or numerical inversion of data. In particular, uniqueness results for the linearized inverse problem can be established, by generalizing a technique initially proposed for the identification of isotropic perturbations of the isotropic electrostatic conductivity coefficient from boundary measurements in a famous paper by Calderon [START_REF] Calderon | On an inverse boundary value problem[END_REF][START_REF] Isaacson | Comment on Calderon's paper: "on an inverse boundary value problem[END_REF]. The main idea in Calderon's approach was to select probing fluxes and trial potential fields so that the right-hand side in the electrostatic equivalent of the LRG identity (33) yields the spatial Fourier transform of the sought conductivity perturbation. This objective can be achieved as well for the present elastic linearized inverse problem by means of the following procedure.

The reference elasticity tensor C is assumed isotropic, i.e. of the form (4). Both fields u and w in [START_REF] Burczy Ński | Sensitivity analysis and identification of material defects in dynamical systems[END_REF] solve the field equation (34a) of elastic equilibrium, and for that reason can be expressed by means of the Galerkin representation [START_REF] Gurtin | The Linear Theory of Elasticity[END_REF], i.e. in the form

u G (x) = 2(1 -ν)∆g(x) -∇div g(x)
where ∆∆g = 0 (35

)
where g is a biharmonic vector potential and ν = λ/(2λ + 2µ) is the Poisson ratio. Here, particular choices of u and w are defined through vector potentials g set up by means of an extension of the Calderon procedure. Let

z = 1 2 (y + iy ⊥ ) z = 1 2 (y -iy ⊥ ) (36)
where y ∈ R 3 , y.y ⊥ = 0, |y| = |y ⊥ |; the overbar denotes the complex conjugation. As a consequence of these definitions, one has z + z = y, z.z = z. z = 0, z. z = |y| 2 /2; also, for a given choice of y, the vector y ⊥ must lie in the plane orthogonal to y and is therefore characterized by one angle ϑ in that plane. Hence, the complex vector parameter z is defined by the couple (y, ϑ). Using these definitions, potentials g[z](x) of the following form are introduced:

g[z](x) = ϕ[z](x)e -iz.x (37)
for which it is worth noting that ∆ x e -iz.x = 0 as a consequence of the property z.z = 0 (note that the vector function ϕ may depend on the parameter z). The requirement that g[z] be biharmonic, equation (37b), is found to translate into the following equation for ϕ[z](x):

∆ x ∆ x ϕ[z] -4i∇ x (∆ x ϕ[z]).z -4(∇ x ∇ x ϕ[z]) : (z⊗z) = 0 (38)
Now, let ϕ[z] be any vector function satisfying [START_REF] Charalambopoulos | The linear sampling method for the transmission problem in three-dimensional linear elasticity[END_REF], g[z] the corresponding Galerkin potential given by ( 37) and u G [z] the displacement associated to g[z] by [START_REF] Calderon | On an inverse boundary value problem[END_REF]. The corresponding linearized strain tensor ε G [z], defined by (1), has therefore the form

ε G [z](x) = (1 -ν)(∇ x + ∇ T x )∆ x g[z](x) -∇ x ∇ x div g[z](x) = Φ[z](x)e -iz.x (39)
where the second-order tensor field Φ[z](x) is a known combination of z and of ϕ[z](x) and its gradients up to order 3. Now, let the probing and trial fields be chosen as

u 0 (x) = u G [z](x) w(x) = u G [ z](x) (40)
The left-hand side of the LRG identity [START_REF] Burczy Ński | Sensitivity analysis and identification of material defects in dynamical systems[END_REF] is then a known quantity H(y, ϑ), and (33) takes the form of a spatial Fourier transform:

H(y, ϑ) = Ω Φ[z](x) : δC(x) : Φ[ z](x)e -iy.x dV x = F Φ[z](•)⊗Φ[ z](•) (y) ⋆ F[δC(•)](y) (41)
where F[f ](y) denotes the spatial three-dimensional Fourier transform of f (x) and the notation ⋆ combines convolution and inner product. Therefore, one observes that performing experiments for all y ∈ R 3 and one value of ϑ provides one scalar relationship on δC(x) of the form [START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF]. It is worth noting that [START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF] is not in general a convolution equation because the parameter z depends on y but is unaffected by the Fourier transform being performed. The procedure proposed by Ikehata [START_REF] Ikehata | Inversion formulas for the linearized problem for an inverse boundary value problem in elastic prospection[END_REF] for the identification of isotropic perturbations δλ(x), δµ(x) of an isotropic reference medium characterized by λ(x), µ(x) is a particular instance of the above approach based on two specific choices of ϕ[z], namely

ϕ (1) [z] = 2 z/|y| 2 ϕ (2) [z] = 2( z.x) z/|y| 2 (42)
Equation ( 41) then provides the two identities

H (1) (y, ϑ) = - |y| 4 2 F[δµ](y) (43)
H (2) (y, ϑ) = -(1 -2ν) 2 |y| 4 F[δλ + δµ](y) - |y| 4 4 F (z.x)( z.x)δµ (y) (44)
where

H (i) (y, ϑ) (i = 1, 2) denote the values taken by H(y, ϑ) for ϕ[z] = ϕ (i) [z],
respectively. Thus, equations ( 43) and ( 44) yield the Fourier transforms of δµ(x), and of δλ(x) with δµ(x) known, respectively. Note that equations ( 43) and ( 44) do not depend on ϑ, which was to be expected as a consequence of the assumed constitutive isotropy. Investigation of the identifiability of small anisotropic perturbations δC of isotropic elastic moduli is amenable to the same treatment. In particular, Galerkin representations [START_REF] Calderon | On an inverse boundary value problem[END_REF] can still be used as the reference medium is assumed isotropic. To the best knowledge of the authors, only partial results of this type are currently known. One such result [START_REF] Bui | Sur quelques problèmes inverses élastiques en mécanique de l'endommagement[END_REF] stems from the observation that the choice ϕ[z] = ϕ (1) [z] leads to

Φ[z](x) : δC(x) : Φ[ z](x) = -(z⊗z) : δC(x) : ( z⊗ z) = -(y⊗y -y ⊥ ⊗y ⊥ ) : δC(x) : (y⊗y -y ⊥ ⊗y ⊥ )/16 -(y⊗y ⊥ + y ⊥ ⊗y) : δC(x) : (y⊗y ⊥ + y ⊥ ⊗y)/16 (45)
Then, for a fixed choice of y, let (a, b, c) denote an orthonormal frame such that y = yc, and put y ⊥ = y(a cos ϑ+b sin ϑ). Equation [START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF] then is a fourth-degree polynomial in cos ϑ, sin ϑ with all the monomials being of even total degree, and hence can be recast in the form

Φ[z](x) : δC(x) : Φ[ z](x) = A 0 + 2 p=1
A p cos 2pϑ + B p sin 2pϑ :: δC(x) [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF] where A 0 , A p and B p are fourth-order tensors that do not depend on ϑ. Therefore, using ϕ[z] = ϕ (1) [z] in [START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF] for all y and all ϑ allows to find at most five independent linear combinations of the elastic coefficients δC ijkℓ (x). Therefore, this set of relations obviously does not allow to identify the 21 independent coefficient of the most general anisotropic perturbations δC(x). In particular, one has Tr(Φ[z](x)) = 0, so that the perturbation of the compressibility modulus cannot be identified using ϕ[z] = ϕ (1) [z].

An identity of the form ( 46) is also found from the choice ϕ[z] = ϕ (2) [z], but this time with a summation range of 1 to 3, so that up to seven independent combinations of the elastic coefficients δC ijkℓ (x) can be resolved. The generalization of this procedure to the identifiability of general anisotropic perturbations δC is currently under investigation by the present authors.

Identification of cracks using the reciprocity gap

Crack identification problems consist in identifying a crack (or a set of cracks) from a set of overdetermined force-displacement boundary measurements. Such problems can be formulated within different physical contexts such as electrostatics, elastomagnetism, acoustics or elastodynamics. They have important applications in e.g. nondestructive testing or the identification of seismic faults. Many strategies have been proposed for solving crack identification problems, some of which are discussed in the section 6 devoted to geometrical sensitivity techniques. In particular, reciprocity gap (RG) functionals can be formulated and applied to crack identification problems. This section aims at summarizing a few recent results concerning the RG-based approach to crack identification.

Let us consider a body Ω, with elastic constitutive properties characterized by a known moduli tensor C, containing a traction-free embedded crack defined by the open surface Γ This corresponds to the idealized crack model, where both crack faces Γ ± are the same geometrical surface Γ with opposite orientations, the displacement is allowed to jump across Γ, and the possibility of contact between crack faces induced by crack closure is disregarded (see however [START_REF] Stavroulakis | Nondestructive elastostatic identification of unilateral cracks through BEM and neural networks[END_REF], where crack identification techniques incorporating unilateral constraints associated with crack closure are proposed). The cracked body is probed by means of known tractions φ applied on the external boundary S, inducing an elastodynamic state in Ω \ Γ characterized by the displacement u(x, t), which solves the boundary-initial value (direct) problem

[A C u](x, t) -ρ(x) ü(x, t) = 0 in Ω (47) u(x, 0) = 0 , u(x, 0) = 0 in Ω p C [u](x, t) = φ(x, t) on S , p C [u](x, t) = 0 on Γ ±
The supplementary data needed for the identification of the unknown crack is assumed to consist of a measurement ξ of the displacement over the entire external boundary, i.e.

u(x, t) = ξ(x, t) (on S) (48)
An appropriate definition of the RG (32) can be established by considering another solution w(x, t) of the homogeneous Navier equation [START_REF] Constantinescu | Sur l'identification des modules élastiques[END_REF], this time defined for the uncracked body (i.e. such that [[w]] = 0 across Γ). Then, formulating the virtual work principle [START_REF] Ben Ameur | Level set methods for geometric inverse problems in elasticity[END_REF] for the unknown u(x, t) with w(x, t) as the virtual field, then for the unknown w(x, t) with u(x, t) as the virtual field, integrating the resulting equations over the time interval 0 ≤ t ≤ +∞ and subtracting the identities thus obtained leads to

R(Γ; u, w) = +∞ 0 Γ [[u]].p[w] dV dt = +∞ 0 S [ξ.p[w] -w.φ] dS dt + Ω\Γ [u. ẇ -w. u] dV (49)
where [[u]] = u +u -denotes the crack opening displacement, i.e. the displacement jump across Γ, and the unit normal on Γ is defined as n = n -. Similar identities are available in the literature in connection with other types of field equations. In particular, the elastostatic RG identity is obtained from ( 49) by removing the integrations in time and suppressing the last term involving time derivatives.

In several situations, a sequential application of the RG identity to three sets of wellchosen test fields w has been shown to permit a complete identification of a planar crack. The first two sets are constructed so that R(Γ; u, w) yield the normal vector to Γ and the distance of the crack plane P to the origin, respectively. Once P is known, the third set is defined so as to express the crack opening displacement [[u]] in the form of either an expansion on a L 2 -orthonormal basis or a Fourier transform on P. Explicit formulae for such sets of test fields have been established so far for the identification of planar cracks in the context of electrostatics [START_REF] Andrieux | Reciprocity principle and crack identification[END_REF][START_REF] Bannour | A semi-explicit algorithm for the reconstruction of 3D planar cracks[END_REF], time-domain acoustics [START_REF] Bui | Diffraction acoustique inverse de fissure plane : solution explicite pour un solide borné[END_REF], transient heat conduction [START_REF] Ben Abda | Reciprocity principle and crack identification in transient thermal problems[END_REF] and elastostatics [START_REF] Andrieux | Sur l'identification de fissures planes via le concept d'écart à la réciprocité en élasticité[END_REF]. This strategy is also applicable to the case of multiple planar cracks lying in the same plane P.

Another recently developed approach [START_REF] Bui | Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap[END_REF] is based on an instantaneous version of the reciprocity gap

R(Γ, t; u, w) = Γ [[u]].p[w] dS = S [ξ.p[w] -w.φ] dS (50)
where w(x, t) are test fields chosen such that the stress field σ[w](x, t) is a travelling Dirac impulse. The application of the RG identity (50) measures the instantaneous mechanical work done by the tractions of test-field on the the crack opening displacement [[u]], referred to as the instantaneous RG. As a consequence, the RG ( 50) is zero until the wavefront of the test field reaches the crack. This remark permits the definition of several numerical techniques to identify the position of a crack. If the test fields are plane waves, the approach outlined above leads to the identification of the convex hull of an unknown object, or a set thereof. In figure 1 we represent the identified convex hull of two cracks by the wavefronts of the test-field. The colours of the incoming wavefronts denote the value of the instantaneous RG and the plot is stopped when the instantaneous RG becomes nonzero.

Direct applications of the virtual work principle

In some situations, measurements of field variables may be available not only on the boundary but also over the whole sample body. This is in particular the case for measurements made on thin samples (for which plane stress conditions may be assumed). In addition, experimental techniques based on X-ray (micro)tomography allow the non-invasive measurement of displacement fields inside three-dimensional deformable solids [START_REF] Flannery | Threedimensional microtomography[END_REF]. In such situations, the virtual work principle provides a direct link between the experimental data and the unknown quantities, which are usually related to constitutive properties, defects or damage. For instance, consider the situation where the elasticity tensor C ⋆ (x) is unknown and the kinematical response (u, ε) to a loading φ applied over the entire boundary S under quasi-static conditions is known over Ω. The measured strain field ε m can be linked to C ⋆ (x) by means of the virtual work identity [START_REF] Ïm | Identification of certain polynomial nonlinear structures by adaptive selectively-sensitive excitation[END_REF]:

Ω ε m (x) : C(x) : ε[w](x) dV = S φ(x).w(x) dS (∀w) (51)
where the (continuous) virtual fields w(x) can be arbitrarily chosen. Each choice of w(x) yields a linear scalar equation for C ⋆ (x). This is the basis of the so-called virtual fields method [START_REF] Gr Édiac | Numerical issues in the Virtual Fields Method[END_REF][START_REF] Gr Édiac | Special virtual fields for the direct determination of material parameters with the virtual fields method. 1-Principle and definition[END_REF][START_REF] Gr Édiac | Special virtual fields for the direct determination of material parameters with the virtual fields method. 2-Application to in-plane properties[END_REF]. To enhance the effectiveness and efficiency of this technique, a variety of techniques have been proposed for the construction of families of virtual fields w(x) tailored for specific classes of constitutive parameter identification problems. A related approach based on the availability of displacement fields rather than strain fields, the equilibrium gap method, can be found in [START_REF] Claire | A finite element formulation to identify damage fields[END_REF].

The virtual work principle applied to parameter sensitivity analysis

Least-squares functionals

In a variety of practical cases one can assume that the quantity to be identified (e.g. elastic moduli, boundary tractions of displacements, cracks. . . ) can be adequately represented by means of a finite number B of scalar parameters, collectively denoted as a B-vector b ∈ R B . Moreover, the available experimental data is also often discrete in nature due to a variety of practical limitations. In such cases, it is not possible, even within a reasonable degree of approximation, to use strategies such as the ones presented in the previous section, which rely on the theoretical assumption of complete data. Let d ∈ R D denote the vector of experimental data values. In the present context, each measurement d I (1 ≤ I ≤ D) can usually be considered as a scalar quantity of the form D I u produced by the application of a linear operator D I to the real displacement u. For instance:

D I u = u(x s ).n(x s )
is the normal displacement recorded at a sensor located at point x s on the boundary, while

D I u = S m σ kℓ [u]n ℓ dS = S m p k [u] dS = F k
defines the k-component of the resultant force over some measurement surface S m . These measurements result from the application of a known excitation φ (e.g. a distribution of forces over the boundary).

⋆ = arg min b∈B J (b) J (b) = J(u[φ, b], b; d) (52)
which embodies a definition of the "best fit" between measured and simulated data, and possibly some additional prior information. For example, regularized least-squares cost functions of the form

J(u, b; d) = 1 2 M i=1 |D I u -d I | 2 + αP(b) (53)
are widely used. Other functionals arise from more sophisticated principles such as Bayesian approaches to inversion [START_REF] Tarantola | Inverse problem theory[END_REF].

The minimization of J (b) can be performed by means of any of a vast array of minimization algorithms. In the present context of elasticity (and more generally of the mechanics of deformable solid bodies), as well as in many other situations where the direct problem is formulated in terms of boundary-or initial-value problems, the simulated measurements D I u[φ, b] depend on b in an implicit way. This remark has the following consequences:

(a) Each evaluation of J (b) requires one direct solution. Such direct solutions may be computationally expensive, e.g. for complex three-dimensional configurations or dynamical problems, or when the direct problem is nonlinear due to e.g. inelastic constitutive behaviour or contact. (b) This in turn often makes search methods such as evolutionary algorithms, which usually require large number of cost function evaluations, computationally impractical.

Therefore, more traditional gradient-based methods are often used for performing the minimization [START_REF] Cowin | Eigentensors of linear elastic materials[END_REF]. To evaluate the gradient ∇ b J (b), one may resort to (i) numerical differentiation, (ii) direct analytical differentiation, or (iii) an adjoint solution. Approach (i), which entails setting up and solving at least B new direct problems on new configurations obtained by adding small but finite perturbations to each parameter b i in turn, is straightforward to implement but usually time-consuming and also prone to inaccuracies because of the ill-posed character of numerical differentiation. Approach (ii) requires solving B "derivative problems" resulting from the differentiation of the governing equations of the direct problem with respect to each b i (1 ≤ i ≤ B); it substantially improves on (i) in that all B derivative problems are governed by the same operator as the original direct problem. It is in particular well suited to direct problems involving nonlinearities and evolution in time [START_REF] Forestier | 3D inverse analysis model using semianalytical differentiation for mechanical parameter identification[END_REF]. Approach (iii), which replaces the solution of B derivative problems with that of one adjoint problem, is the most efficient when the direct problem is linear, while being applicable to other situations as well. In section 4.2, the adjoint solution approach is presented for quasi-static inversion problems in elasticity. The virtual work principle of continuum mechanics will appear to be very helpful in setting up computationally efficient sensitivity formulae, based on adjoint states, for computing the gradient ∇ b J (b). An extension of this approach to identification problems involving contact, i.e. an unilateral constraint, is then described in section 4.3

Adjoint solution for the identification of elastic moduli

As an example, the adjoint state approach is formulated in this section for the problem of identifying distributions of elastic moduli from elastostatic boundary measurements. Let u = u[C] denote the displacement created in the body endowed with a trial elastic tensor C(x) by the application of known forces φ over S. By virtue of the virtual work principle, u ∈ V solves the weak formulation

Ω ε[u] : C : ε[w] dV - S φ.w dS = 0 ∀w ∈ V ( 54
)
where V is the space of displacement fields having a bounded strain energy over Ω (in practice, additional constraints should be imposed on u so as to prevent spurious rigid-body motions, a secondary issue which is here left aside for the sake of simplicity).

Assuming that displacement measurements ξ are available over a measurement subset S m ⊆ S, a generic cost functional J (C) is considered:

J (C) = J(u[C], C) = S m ϕ(u[C]; ξ) dS + Ω ψ(C) dV (55)
where the first integral measures a distance between u[C] and the measurement ξ (e.g. ϕ(u; ξ) = 1 2 |u -ξ| 2 ) and the second integral may be used to express some prior information on the sought moduli distribution. Note that the format (55) allows for incomplete data, and even pointwise measurements, whereas RG functionals such as [START_REF] Bui | Numerical identification of linear cracks in 2D elastodynamics using the instantaneous reciprocity gap[END_REF] require complete data.

The first variation of J (C) is given by

δJ (C) = S m ∂ϕ ∂u .δu dS + Ω ∂ψ ∂C : δC dV ( 56
)
where δu is the variation induced by δC. The governing problem for the latter is found by differentiating the weak formulation [START_REF] Dassios | Low frequency scattering[END_REF]:

Ω ε[δu] : C : ε[w] dV + Ω ε[u] : δC : ε[w] dV = 0 ∀w ∈ V (57)
Now, let (57) be restricted to trial function w such that A C w = 0. In that case, the virtual work principle with δu as the virtual field yields

Ω ε[δu] : C : ε[w] dV - S q.δu dS = 0 (58) having set q = σ[w].n = p[w]
and where the constitutive symmetry (3) has been used. Subtracting the latter identity from (57), one obtains

S q.δu dS = - Ω ε[u] : δC : ε[w] dV
Comparing the above identity with equation ( 56), the variation of J (C) is therefore expressed as

δJ (C) = Ω ε[u] : δC : ε[w] dV + Ω ∂ψ ∂C : δC dV ( 59
)
where the adjoint solution w solves the problem

A C w = 0 (in Ω), p[w] = - ∂ϕ ∂u (u[C]) (on S m ), p[w] = 0 (on S \ S m ) (60)
or, equivalently, the weak formulation

Ω ε[w] : C : ε[v] dV + S m ∂ϕ u ∂u (u[C], ξ).v dS = 0 ∀v ∈ V (61)
The parameter sensitivity formula (59) treats C(x) as the main unknown. In practice, the sought distribution of moduli is often described as a known function Ĉ of a finite set of unknown parameters b, i.e. C(x) = Ĉ(x; b), in which case equation ( 59) can be used to formulate partial derivatives with respect to individual parameters b i of Ĵ (b) = J ( Ĉ(•; b)):

∂ Ĵ ∂b i = Ω ε[u] : ∂ Ĉ ∂b i : ε[w] dV + Ω ∂ψ ∂C : ∂ Ĉ ∂b i dV ( 62
)
The optimal control viewpoint. The parameter sensitivity formula (59) could alternatively have been established by introducing the Lagrangian

L(u, w, C) = J(u, C) + Ω ε[u] : C : ε[w] dV - S φ.w dS (63)
where the direct problem (in weak form) is treated as a constraint, the trial field w acting as a distributed Lagrange multiplier. Then, the constrained minimization of J (C) is recast as a saddle point search for L(u, w, C), characterized by the stationarity equations

∂L ∂u .δu = 0 (∀δu) ∂L ∂w .δw = 0 (∀δw) ∂L ∂C .δC = 0 (∀δC) (64)
For an assumed value of C, the first two stationarity equations are readily seen to define the direct problem [START_REF] Dassios | Low frequency scattering[END_REF] and the adjoint problem [START_REF] Feij Óo | An application of shape optimization in the solution of inverse acoustic scattering problems[END_REF]. Formula (59) expressing the variation of J (C) is then obtained by substituting the direct and adjoint solutions into the third stationarity equation in [START_REF] Forestier | 3D inverse analysis model using semianalytical differentiation for mechanical parameter identification[END_REF]. This Lagrangian-based approach to the formulation of adjoint problems and parameter sensitivity formulae is applicable to many other situations and widely used.

An alternative approach consists in solving at once the complete set of equations [START_REF] Forestier | 3D inverse analysis model using semianalytical differentiation for mechanical parameter identification[END_REF]. It is especially convenient when all equations (64) are linear, which may happen in other types of inverse problems amenable to the same type of formulation, as the inverse problem can then be solved by means of one calculation. However, this approach leads to a global matrix formulation which is not standard, so that an additional programming effort is needed for the implementation of this approach in standard FEM or BEM packages. An illustration of the application of this technique to determine an unknown boundary temperature from interior temperture measurements is presented in the thermostatic case in [START_REF] Constantinescu | Identification du chargement thermique transitoire par contrôle optimal[END_REF][START_REF] Delattre | Application du contrôle optimal à l'identification d'un chargement thermique[END_REF]. If the complete system (64) is nonlinear, its solution normally requires an iterative method.

Indentation test and adjoint state for unilateral contact

The indentation test is used for the determination of various types of mechanical parameters, in particular those associated with constitutive properties of materials. The test consists in pressing a punch (the indentor) on the surface of a deformable sample. During the test both the indentation depth U (t) and the resultant force F (t) applied by the indentor are recorded (the loci (U (t), F (t)) defining the indentation curve). Even for linear elastic materials, the indentation curve is nonlinear due to the unilateral contact conditions and (depending on the properties of the sample surface) the frictional forces. The practical usefulness of the test lies in its simplicity and the fact that no specimens have to be prepared. However, no closed-form solutions are available for the interpretation of indentation curves (see [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: study of variational inequalities and finite element methods[END_REF]). Hence, to infer the sought material properties from an indentation curve requires inversion algorithms where the corresponding direct problem consists in (numerically) simulating the indentation curve associated with assumed (i.e. trial) material parameters.

For definiteness, consider the identification of elastic properties that depend on a finite set b of constitutive parameters, i.e. C = Ĉ(x; b), by slowly indenting a sample with a rigid punch. Such experimental procedures are actually carried out on e.g. thin films or layered materials. The inverse problem consists therefore in recovering b from an indentation curve (U (t), F (t)). A more general discussion of this type of identification problem is given in [START_REF] Tardieu | On the determination of elastic coefficients from indentation experiments[END_REF].

Let the indentation depth U (t) be prescribed, under quasi-static conditions, and assume that the contact is frictionless (i.e. that only normal forces can develop on the contact surface Γ C ). The elastic deformation of the sample at time instant t, under quasi-static conditions, can be obtained by seeking the displacement field u that minimizes the potential energy subject to the non-penetration inequality constraint, i.e.:

u = arg min v∈V 1 2 Ω ε[v] : Ĉ : ε[v] dV subject to u n -g -U (t) ≤ 0 on Γ (65)
where u n = u.n, Γ ⊂ S is the potential contact surface (i.e. it is assumed that Γ C (t) ⊆ Γ, where Γ C (t) is the actual contact surface at time t), V is again the space of displacements having a bounded strain energy over Ω, and g is the initial gap, measured along the direction normal to Γ, between Γ and the punch. Due to the non-penetration inequality constraint, the actual contact surface Γ C (t) is not known a priori. Necessary (first-order) stationarity conditions for the minimization problem are:

Ω ε[u] : Ĉ : ε[w] dV - Γ pw n dS = 0 (∀w ∈ V) (66)
p ≤ 0 , u n -g -U (t) ≤ 0 , (u n -g -U (t))p = 0 (on Γ) (67)
where the contact pressure p = n.p[u] = n.σ [u].n = n.( Ĉ : ∇u).n appears as the Lagrange multiplier associated with the non-penetration inequality constraint. Once problem ( 66)-( 67) is solved, the actual contact surface Γ C (t) and the contact pressure distribution p for given indentation depth U (t) and elastic properties Ĉ are known.

Let the identification be performed on the basis of the cost function

J (b) = J(F [U, b]) = 1 2 T 0 [F [U, b](t) -F (t)] 2 dt ( 68
)
where the computed resultant contact force F [U, b](t) on the actual contact surface Γ C , given in terms of the contact pressure p associated with the solution u[U, b] of ( 66)-( 67) by

F [U, b](t) = Γ C (t) p(x, t) dS (69)
is compared to the measured resultant contact force F (t). The variation of J (b) induced by a variation δb is:

δJ (b) = T 0 F [U, b](t) -F (t) δF [U, b](t) dt δF [U, b](t) = Γ C (t) δp(x, t) dS (70)
where δp(x, t) is the variation of contact pressure induced by δb. On differentiating the first-order optimality conditions ( 66) and (67c), the variation of the solution to the direct problem ( 66)-( 67) is found to be governed by the weak formulation

Ω ε[δu] : C : ε[v] dV + Ω ε[u] : δC : ε[v] dV - Γ C δpv n dS = 0 (∀v ∈ V) (71)
(u n -g -U (t))δp + pδu n = 0 (on Γ) (72)
where, as in Section 4.2, δC = (∂ Ĉ/∂b).δb, and provided strict complementarity holds (i.e. that the inequalities are strict almost everywhere) in [START_REF] Gao | On the inversion of residual stresses from surface measurements[END_REF]. Now, like in section 4.2, only trial functions v such that Av = 0 are considered. Identity ( 58) holds for such v, and can be combined with ( 71) to obtain

Ω ε[u] : ( ∂ Ĉ ∂b δb) : ε[v] dV + S q.δu dS - Γ δpv n dS = 0
Since (72) together with the complementarity conditions [START_REF] Gao | On the inversion of residual stresses from surface measurements[END_REF] imply that δu n = 0 on Γ C (t) and δp = 0 on Γ \ Γ C (t), this identity becomes

Ω ε[u] : δC : ε[v] dV + S\Γ C (t) q.δu dS + Γ C (t) q ⊥ .δu ⊥ dS - Γ C (t)
δpv n dS = 0 [START_REF] Gr Édiac | Special virtual fields for the direct determination of material parameters with the virtual fields method. 2-Application to in-plane properties[END_REF] where z ⊥ = z -(z.n)n for any z ∈ R 3 . Let the adjoint state be defined as the elastostatic state (v, q) satisfying the same homogeneous boundary conditions than u on S \ Γ C together with the boundary conditions

v n = [F (u, b) -F (t)] q ⊥ = 0 (on Γ C (t)) (74)
which ensure that the second and third integrals in (73) vanish. On adding the resulting identity to [START_REF] Geymonat | Identification of mechanical properties by displacement field measurement: a variational approach[END_REF], the sensitivity of the cost function ( 68) is given in terms of the direct and adjoint solutions at all times t as

δJ (b) = T 0 Ω ε[u] : δC : ε[v] dV dt (75)
This technique has been extended to the indentation of materials with nonlinear constitutive properties and applied to experimental data [START_REF] Constantinescu | On the identification of elastoviscoplastic constitutive laws from indentation tests[END_REF]. It is worth noting that the adjoint problem is defined in terms of a bilateral constraint (prescribed normal displacement). It is therefore linear, and in particular simpler than the direct contact problem. A mathematical discussion of optimality conditions in optimal control problems where the state equations feature unilateral constraints of the no-penetration type can be found in [START_REF] Bergounioux | On the structure of Lagrange multipliers for stateconstrained systems[END_REF].

Formulations based on the error in constitutive equation

As mentioned earlier, many inverse problems in solid mechanics, and in other areas of physics as well, revolve around the identification of constitutive parameters, i.e. parameters appearing in the constitutive relations modelling the physical behaviour of materials. Among the key issues in that area is the identification of distributed parameters. The sets of all admissible displacements and stresses are defined e.g. by ( 22) and ( 23) with reference to equilibrium equations and boundary conditions, and admissible strains are then obtained from the compatibility equation [START_REF] Allix | A new identification strategy for dynamic problems with uncertain measurements[END_REF]. Unknown constitutive parameters appear in the constitutive equations, which provide the link between strains and stresses. It is therefore natural to consider cost functions based on a notion of error in constitutive equation (ECE). This concept is, either implicitly or explicitly, at the core of several studies aimed at the identification of conductivity-type distributed parameters. ECE-based cost functions are involved in e.g. [START_REF] Wexler | An impedance computed tomography algorithm and system for ground water and hazardous waste imaging[END_REF] (identification of permeability distributions) or [START_REF] Kohn | Numerical implementation of a variational method for electric impedance tomography[END_REF] (identification of electrostatic conductivity distributions). Their convexity properties have been studied extensively for the electrostatic case in [START_REF] Kohn | Determining conductivity by boundary measurements[END_REF][START_REF] Kohn | Determining conductivity by boundary measurements II. Interior results[END_REF][START_REF] Kohn | A Variational Method for Parameter Identification[END_REF], where the ECE for anisotropic conductivities is shown to be the quasiconvexification of that defined for isotropic conductivities. Independently, the concept of ECE has first been introduced in elasticity by Ladevèze and Leguillon [START_REF] Ladev Èze | Error estimates procedures in the finite element method and applications[END_REF] in connection with error estimation in finite element computations. A more general definition, applicable to non-linear and history-dependent constitutive properties, has been since formulated in terms of the so-called Drucker stability inequality [START_REF] Ladev Èze | Nonlinear computational structural mechanics. New approaches and non-incremental methods of calculation[END_REF]. The concept of ECE in solid mechanics has proved to be very fruitful not only in connection with its initial motivation (namely deriving error estimation techniques in computational solid mechanics) but also in parameter identification problems.

This section aims at presenting ECE-based formulations for parameter identification in elasticity in some detail. The basic concept is formulated, together with a parameter sensitivity formula, in section 5.1. Then, section 5.2 is devoted to an identification algorithm based on alternating directions. The extension of these concepts is then presented in section 5.3, and several variations around the ECE in elastostatics are summarized in section 5.4. Finally, an ECE-based approach for model updating from vibrational measurements is described in section 5.5, and its capability for prior localization of defects is discussed.

Error in constitutive equation in elastostatics

The variational principles of elasticity are very useful in defining cost functions based on the error in constitutive equation (ECE), which are well suited to the identification of distributed parameters. The main properties of such cost functions are now discussed in more detail and illustrated, in connection with elastostatics (section 5.1.1) and elastic vibrations (section 5.5).

ECE functionals.

Let us again consider the problem of reconstructing the unknown distribution of elastic moduli C ⋆ from a series of overspecified boundary data pairs of displacements and traction (ξ, φ).

One can introduce and motivate ECE functionals by considering the sum E(v, s, C) of the potential energy [START_REF] Bergounioux | On the structure of Lagrange multipliers for stateconstrained systems[END_REF] and the complementary energy [START_REF] Berthaud | Experiment/computation interactions by using digital image correlation[END_REF] for an assumed elasticity tensor C and a well-posed (i.e. not overspecified) set of boundary conditions of the form (12-13):

E(v, s, C) = W C (v) + W ⋆ C (s) = 1 2 Ω ε[v] : C : ε[v] dV + 1 2 Ω s : C -1 : s dV - S v.(s.n) dS (76)
If the admisible fields are such that v ∈ C(φ, S p ) and s ∈ S(ξ, S u ), with the spaces C(•, •) and S(•, •) defined by ( 22) and ( 23), an integration by parts and some algebraic manipulations permit to recast E(v, s, C) as a functional measuring the gap in the constitutive law, in either of the equivalent forms

E(v, s, C) = 1 2 Ω |C -1/2 : s -C 1/2 : ε[v]| 2 dV (77) = 1 2 Ω s -C : ε[v] : C -1 : s -C : ε[v] dV (78)
In particular, in view of ( 78), the statement [START_REF] Bonnet | Sensitivity analysis for shape perturbation of cavity or internal crack using BIE and adjoint variable approach[END_REF] means that an admissible pair (u, σ) solves an elastic equilibrium problem if and only if they are related through the elastic constitutive equation ( 8), as expected. Equivalently, the absolute minimum E(u, σ, C) = 0 is reached only when the assumed elasticity tensor C is the correct one. These remarks suggest to introduce the following definition for E(C), the error in constitutive equation (ECE) functional for the identification of a distribution of elastic moduli:

E(C) = arg min (v,s)∈C(ξ,S)×S(φ,S) E(v, s, C) (79)
in which an important difference with [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] is that overspecified boundary data of the form ( 27) is now involved. In general, i.

[u D,N ]) = 0 (in Ω), u D = ξ ; C : ε[u N ] .n = φ (on S) (80)
One can then easily show through integration by parts and algebraic manipulations that E(C) can be expressed in any of the following equivalent forms:

E(C) = 1 2 Ω ε[u D ] -ε[u N ] : C : ε[u D ] -ε[u N ] dV (81) = 1 2 Ω σ[u D ] -σ[u N ] : C -1 : σ[u D ] -σ[u N ] dV (82) = 1 2 S σ[u D ].n -φ .(ξ -u N ) dS (83)
Note that E(C) expressed by either of the above formulae depends on C explicitly, but also implicitly through u D and u N .

Gradient of the ECE functional.

The gradient of E(C) is given by

δE(C) = ∂E ∂v • δu D + ∂E ∂s • δσ[u N ] + ∂E ∂C • δC (v,s,C)=(u D ,σ[u N ],C)
in order to take into account the implicit dependence on C of E(C) through u D and u N . Besides, u D and u N solve the weak formulations

∂E ∂v • δv = 0 , ∀δv ∈ C(S, 0) ; ∂E ∂s • δs = 0 , ∀δs ∈ S(S, 0)
Hence, the gradient of E(C) is expressed in terms of the solutions u D and u N as

δE(C) = ∂E ∂C (u D , σ[u N ], C) = 1 2 Ω ε[u D ] : δC : ε[u D ] -ε[u N ] : δC : ε[u N ] dV (84)

Identification algorithm based on alternating directions

The lack of mathematical results for the properties of E in elasticity notwithstanding, this functional has been used to define an alternating directions algorithm for the numerical reconstruction of distributed elastic moduli [START_REF] Constantinescu | On the identification of elastic moduli from displacementforce boundary measurements[END_REF], along the lines proposed in [START_REF] Kohn | Numerical implementation of a variational method for electric impedance tomography[END_REF] for isotropic conductivities. Starting from an initial guess C (0) (x), a sequence C (i) (x) is computed, the iterate C (i+1) (x) being obtained from the previous iterate C (i) (x) through the following steps:

(i) compute u D(i) and u N(i) by solving the problems (80) with C = C (i) , and then σ N(i) = C (i) : u N(i) ;

(ii) compute the new iterate C i+1 from

C i+1 = arg min C E(C, u D(i) , σ N(i) )
To simplify step (ii), it is useful [START_REF] Constantinescu | On the identification of elastic moduli from displacementforce boundary measurements[END_REF] to express the elasticity tensor in the form

C(x) = 6 k=1 c k (x) ζ k (x) ⊗ ζ k (x) (85)
in terms of its eigenelastic moduli c k > 0 and its eigentensors ζ k [START_REF] Cowin | Eigentensors of linear elastic materials[END_REF], chosen such that ζ k : ζ ℓ = δ kℓ in order to define an orthonormal family of symmetric second-order tensors. Then, one obtains from a simple algebraic calculation that

E(C, u D(i) , σ N(i) ) = 1 2 6 k=1 c k ε (i) k : ε (i) k + c -1 k σ (i) k : σ (i) k -2ε (i) k : σ (i) k with ε (i) k = ε[u D(i) ] : (ζ k ⊗ ζ k ) and σ (i) k = (C (i) : ε[u N(i) ]) : (ζ k ⊗ ζ k ).
The minimization of step (ii) with respect to the eigenmoduli c k is straightforward and yields

c (i+1) k = σ (i) k : σ (i) k ε (i) k : ε (i) k 1/2 (86)
It is important to note that this procedure assumes a given set of eigentensors ζ k , and therefore does not attempt to find the ζ k that are optimal with respect to E(C, u D(i) , σ N(i) ). It is therefore applicable only within certain a priori assumptions regarding material symmetries, so that the ζ k are known beforehand.

If material isotropy is assumed, there are only two distinct eigenmoduli c 1 = λ + 2µ/3 (known as the bulk modulus) and c 2 = 2µ, the shear modulus. The updating formulae [START_REF] Johnson | Contact mechanics[END_REF] for step (ii) become c (i+1) 1

=

Tr(σ

(i) k ) Tr(ε (i) k ) c (i+1) 2 = Dev(σ (i) k ) : Dev(σ (i) k ) Dev(ε (i) k ) : Dev(ε (i) k ) 1/2 (87)
where Dev(s) := s -1 3 Tr(s)I is the deviatoric part of the symmetric second-order tensor s. The previously described algorithm is based on the solution of standard Neumann and Dirichlet elasticity problems [START_REF] Hugues | The finite element method -linear static and dynamic finite element analysis[END_REF]. It can therefore be easily coded within a conventional finite element analysis environment. The examples presented here, programmed within the CAST3M [START_REF]A general purpose code for solving partial differential equations by the finite element method[END_REF] object-oriented finite element package, concern the identification of a copper inclusion of a rectangular or coin shape in an aluminum matrix. Two inclusion shapes (rectangle and wedge) are considered, as indicated respectively on figures 2 and 3. Plane strain conditions are assumed. Both materials are assumed to have anisotropic elasticity properties with cubic symmetry, so that three elastic moduli distributions have to be identified from a series of simulated boundary measurements. The latter result from applying a parabolic distribution of forces over 5 consecutive nodes on the boundary. This loading pattern is moved along the boundary, each possible configuration defining a separate (simulated) experiment. The boundary displacements induced by these loadings have been artificially polluted by white noise of amplitude up to 10%. The finite element mesh for both examples has 24 × 24 regularly spaced nodes.

Extension to elastic plate problems

In this section, an extension of the previously discussed algorithm to elastic moduli identification in thin elastic plates from overdetermined boundary measurements is presented.

Fundamental equations for thin elastic plates. A thin plate is a planar solid occupying a domain of the form

Ω = ω × [-h/2, h/2],
where ω ⊂ R 2 is the planar mean surface and h ≪ Diam(ω) is the thickness. Only flexural motions of plates are considered here, leaving aside stretching motions (i.e. in-plane deformations) which are governed by equations similar to those of section 2.1. Under the Love-Kirchhoff kinematic assumption, the displacement and strain for such motions are expressed in terms of u( x), the out-of-plane displacement of the mean surface, also referred to as the deflection:

u( x, x 3 ) = u( x)e 3 -x 3 ∇u , ε( x, x 3 ) = -x 3 ∇ ∇u ( x ∈ ω, - h 2 ≤ x 3 ≤ h 2 )
where x is the position vector in the planar region ω and ∇ denotes the two-dimensional gradient with respect to x. The equilibrium field equations are

div div M ( x) + f ( x) = 0 ( x ∈ ω) ( 88
)
where f is a distribution of forces applied normally to the plate surface, div is the twodimensional divergence with respect to x, and M denotes the bending moment tensor, i.e. the counterpart of the stress tensor σ, defined in terms of Cartesian components by

M ij ( x) = h/2 -h/2 x 3 σ ij ( x, x 3 ) dx 3 (i, j = 1, 2)
The elastic constitutive equation takes the form

M = D : ( ∇ ∇u) with D = h 3 12 C PS
where D is the fourth-order flexural rigidity tensor and C PS denotes the elasticity tensor for plane stress, i.e. that associated with the linear relationship σ ij = C PS ijkℓ ε kℓ (i, j, k, ℓ = 1, 2) obtained after elimination of ε 13 , ε 23 , ε 33 through enforcing σ 13 = σ 23 = σ 33 = 0 in the threedimensional constitutive equation σ = C : ε (where C is assumed to be independent on x 3 ). In the most general case, the bending of anisotropic elastic plates is therefore described by six independent moduli. The symmetry and positive definiteness of D follow from those of C.

The equilibrium of elastic plates is therefore governed by the fourth-order self-adjoint partial differential equation

A D u( x) + f ( x) = 0 A D u := div div [D : ( ∇ ∇u)] ( x ∈ ω) (89)
Boundary conditions consist in prescribing two quantities at any point x ∈ ∂ω: (i) either the deflection u or the shear force Q, and (ii) either the normal derivative r or the normal moment M N , and hence have the form

u = u D (on S d ), Q = Q D (on S f ); r = r D (on S r ), M N = M D (on S m ) (90)
where S d ∩ S f = ∅, S d ∪ S f = ∂ω, S r ∩ S m = ∅ and S r ∪ S m = ∂ω and with the definitions

r = u ,i n i M N = M ij n i n j Q = n i (2M ij,j -M ij,k n j n k ) ( x ∈ ∂ω)

ECE functional.

Like in three-dimensional elasticity, a potential energy W and a complementary energy W ⋆ can be defined, this time by

W C (u) = 1 2 ω ( ∇ ∇u) : D : ( ∇ ∇u) dV - ω f u dV - S f M D u ds - Sm Q D r ds(91) W ⋆ C (M ) = 1 2 ω M : D -1 : M dV - S d u D Q ds - Sr r D M N ds (92)
and the sum of W and W ⋆ yields the error in constitutive equation:

E(u, M , D) = W C (u) + W ⋆ C (M ) = 1 2 ω ( ∇ ∇u) : D : ( ∇ ∇u) + M : D -1 : M dV - ∂ω [uQ + rM N ] ds = 1 2 Ω M -D : ( ∇ ∇u) : D -1 : M -D : ( ∇ ∇u) dV (93)
Assuming that the overdetermined data is such that (u D , r D , M D , Q D ) are known on the complete boundary ∂ω, the error in constitutive equation (ECE) functional is defined by

E(D) = arg min u∈C,M ∈S E(u, M , D) (94)
C(u D , r D , ∂ω) = u | u = u D and ∂u ∂n = r D on ∂ω S(M D , Q D , ∂ω) = M | div div M = 0 in ω, M N = M D and Q = Q D on ∂ω
where the above definition of C and S, the sets of kinematically and statically admissible fields, is consistent with the overdetermined nature of the boundary data and are the counterparts for plates of equations ( 22) and ( 23). Identifiability and uniqueness issues have been investigated by Ikehata [START_REF] Ikehata | An inverse problem for the plate in the Love-Kirchhoff Theory[END_REF][START_REF] Ikehata | The linearization of the Dirichlet to Neumann map in anisotropic plate theory[END_REF]. Using the Calderon approach, he has established [START_REF] Ikehata | An inverse problem for the plate in the Love-Kirchhoff Theory[END_REF] that the DtN map uniquely determines the flexural rigidity tensor D( x) if constitutive isotropy is assumed. For anisotropic plates, he has shown [START_REF] Ikehata | The linearization of the Dirichlet to Neumann map in anisotropic plate theory[END_REF] that flexural rigidity tensors belong to either of two classes, uniqueness being valid for one but not for the other.

A numerical reconstruction technique for distributed flexural rigidities based on the minimization of the ECE functional [START_REF] Kress | Integral equation methods in inverse obstacle scattering[END_REF], similar to that presented in section 5.2 for threedimensional elasticity, has been proposed and implemented for isotropic thin elastic plates. Again, alternating directions were used. Each iteration involves the solution of two plate bending problems for the current value of D (defined in terms of kinematic data (u D , r D ) and static data (M D , Q D ) on ∂ω, respectively) followed by an explicit updating of D from the curvature tensor ∇ ∇u and the bending moments M just computed. The latter step is again formulated in terms of the eigenelastic stiffnesses c k of D in the representation

D = 3 k=1 c k ζ (k) ⊗ ζ (k)
where ζ (k) are the corresponding eigentensors. The general-purpose FEM code CAST3M has, as with the plane-strain example of Section 5.2, been used as a basis for the implementation. In the following example, a square plate is divided into 2 × 20 × 20 triangular finite elements of DKT (Discrete Kirchhoff Triangle) type [START_REF] Hugues | The finite element method -linear static and dynamic finite element analysis[END_REF]. Synthetic data was generated by solving the plate bending problem for the true distribution of D for nine cases of applied forces, in this case concentrated loads applied at nine different locations. The synthetic measured displacement data consists in the deflection u at the location of the concentrated load and the deflection distribution on ∂ω. In some cases the normal derivative r and the twisting moment have has also been used as data, but this additional information did not produce substantial changes in the results. The results corresponding to the identification of a square inclusion in a plate are displayed in figure 4.

Other ECE-based functionals in elastostatics 5.4.1. Modified ECE functional.

There are many other ways of incorporating the error in constitutive equation into cost functions for the purposes of parameter identification. One natural extension of (79) consists of relaxing the satisfaction of overdetermined boundary data (ξ, φ) in the definition of spaces of admissible fields. One can then consider functionals of the form where the admissible fields now in data ( φ) H(v, ξ, C) defined

H(C) = min (v,s, ξ, φ) H(v, s, ξ, φ, C) (v, s) ∈ C( ξ, S) × S( φ, S) (95)
H(v, φ, C) 1 Ω -: : -1 -: + S 2 -ξ| ξ B | -2 φ ( 96
)
with B and • ξ | | indicating norms. the minimization respect (v, of s, ξ, φ φ, C) corresponds the definition The ECE results finding pair s) with respect ECE the accuracy of boundary data, cases (A, B) = (0, 0) and (A, B) = (∞, ∞) indicating complete disregard and exact enforcement of this data, respectively. Adjusting them by means of the L-curve method [START_REF] Hansen | Rank-deficient and discrete ill-posed problems[END_REF] has not been done so far but is certainly worthy of investigation. Error functionals similar to [START_REF] Kress | Integral equation methods in inverse obstacle scattering[END_REF] have been proposed in [START_REF] Allix | A new identification strategy for dynamic problems with uncertain measurements[END_REF] in connection with the identification of constitutive parameters from dynamical experiments and in [START_REF] Geymonat | Identification of mechanical properties by displacement field measurement: a variational approach[END_REF] for incorporating displacement fields measured in the whole domain instead of just on the boundary.

To illustrate this idea, consider an elastic solid (in the plane-strain framework) occupying the unit square, the reference shear modulus and Poisson ratio being µ = 1 Pa and ν = 0.3, respectively. A FEM discretization is introduced on the basis of a regular mesh made of 35×35 square four-noded isoparametric elements. Simulated data for a square defect occupying the four finite elements located in rows (10 ,11) and columns (11 ,12), and whose elastic moduli are µ def = 2µ and ν = 0.4, are computed for six parabolic distributions of normal loads applied at six different locations. The force distribution is treated as exact data, and B = 50 is used throughout, whereas the displacement recorded on the boundary may be polluted. Contour maps of the modified ECE functional H(C) and of the least-squares cost function

J (C) = 1 2 S |u N -ξ| 2 dS ( 97
)
with u N defined by [START_REF] Hugues | The finite element method -linear static and dynamic finite element analysis[END_REF], both plotted as a function of the location of the center of a trial fourelement square defect (with the trial distribution of moduli C defined accordingly), are shown in figure 5. The modified ECE functional appears to have better convexity in the vicinity of the "true" defect than J (C). In addition, in the presence of (simulated) noise on the displacement data, H(C) with A = 1 is better behaved in the vicinity of the true defect than both the J (C) and H(C) with A = 50.

ECE as a penalty term.

The major drawback of the least squares functionals are generally bad stability properties, in the sense that small data errors induce large errors in the solution. The classical technique for restoring stability is regularization, a subject covered by a vast literature (see [START_REF] Engl | Regularization of inverse problems[END_REF][START_REF] Tikhonov | Solutions to ill-posed problems[END_REF]). Here, we wish confine ourselves interesting of ECE the of identification problems [START_REF] Chavent | Primal-dual formulations for parameter estimation problems[END_REF]. authors that state in problems form be in the of energy-based penalty to up functions of the form (using notation as in equation ( 52) and where η is a penalty parameter) exploiting the fact that u = u[φ, b] is the minimizer of the potential energy W C(b) (v). However, the value of the minimum W C(b) (u) is problem-dependent, making this approach somewhat impractical. The approach proposed in [START_REF] Chavent | Primal-dual formulations for parameter estimation problems[END_REF] defines the penalty term as the sum of the potential and complementary energies rather than just the potential energy, thus taking advantage of properties ( 24), [START_REF] Bonnet | Sensitivity analysis for shape perturbation of cavity or internal crack using BIE and adjoint variable approach[END_REF]. Such penalty term is in fact the error in constitutive equation. Accordingly, the primal-dual formulation [START_REF] Chavent | Primal-dual formulations for parameter estimation problems[END_REF] of the identification problem is then defined in terms of the cost function

η (b) min J(v, b) W (v) 0.
J PD η (b) = min v∈C,s∈S J(v, d, b) + 1 η [W C(b) (v) + W ⋆ C(b) (s)] ( 99
)
In it is known from ( 25) that the minimum attainable by the penalty term is always zero, and therefore problem-independent. Formulations ( 99) and ( 95) are clearly linked, with coefficients (A, B) in ( 95) playing the same role as the penalty parameter η in (99).

Error in constitutive equation in dynamics

The reconstruction of distributed parameters (such as the flexural stiffness or the mass density) of mechanical structures from vibrational data, i.e. measured values of eigenfrequencies and eigenmodal displacements, is a class of inverse problem of engineering interest, especially in connection with updating finite element (FE) models of mechanical structures, i.e. correcting FE models so that they agree best with measurements on the real structure. Theoretical studies [START_REF] Barcilon | Inverse problem for vibrating beam[END_REF][START_REF] Barcilon | On the multiplicity of solutions of the inverse problem for a vibrating beam[END_REF][START_REF] Barcilon | Inverse problem for the vibrating beam in the free-clamped configuration[END_REF] show in particular that the knowledge of all eigenfrequencies for one set of boundary conditions is usually insufficient for a reconstruction of such characteristics, even assuming that the latter are in the form of a small unknown perturbation of a known reference value. On the other hand, mechanical structures of engineering interest are too complex for the reconstruction inverse problem to be amenable to the same kind of mathematical analysis.

Structural FE model updating using measured vibrational data has been the subject of many investigations during the recent years. Estimators of modeling errors and cost functions based on the ECE play a prominent role [START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF][START_REF] Deraemaeker | Reduced bases for model updating in structural dynamics based on constitutive relation error[END_REF][START_REF] Feuardent | Amélioration des modèles par recalage : application aux structures spatiales[END_REF][START_REF] Ladev Èze | Error estimates procedures in the finite element method and applications[END_REF][START_REF] Ladeveze | Parametric correction of finite element models using modal tests[END_REF][START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF]. ECE has also been used in conjunction with Gaussian inversion [START_REF] Ben Abdallah | Inversion gaussienne appliquée à la correction paramétrique de modèles structuraux[END_REF][START_REF] Bonnet | Nonlinear Gaussian inversion applied to structural parameter identification[END_REF]. Other approaches for FE model updating include references [START_REF] Andriambololona | Optimisation des essais et recalage de modèles structuraux[END_REF][START_REF] Ïm | Identification of certain polynomial nonlinear structures by adaptive selectively-sensitive excitation[END_REF][START_REF] Link | Requirements for the structure of analytical models used for parameter identification[END_REF][START_REF] Zimmermann | Structural damage detection using a minimum rank update theory[END_REF] and the survey paper [START_REF] Mottershead | Model updating in structural dynamics : a review[END_REF].

ECE-based functional.

Consider for definiteness the free vibrations of a linearly elastic solid occupying the domain Ω and endowed with the elasticity tensor C ⋆ (x) and mass density ρ ⋆ (x) distributions. The governing problem for such motions has the form

A C ⋆ u ⋆ + ρ ⋆ ω 2 u ⋆ = 0 (in Ω) u = 0 (on S u ) p[u] = 0 (on S p ) (100)
and is known to have a countable set of real eigenvalues 0 < ω 1 ≤ ω 2 ≤ . . ., each of finite multiplicity (0 < ω 1 holds whenever the boundary conditions prevent rigid-body motions).

Assume that imperfect values ρ(x), C(x) of the true distributions of mass density ρ ⋆ (x) and elasticity tensor C ⋆ (x) are known. The model updating problem then consists of finding corrections ∆ρ = ρ ⋆ρ, ∆C = C ⋆ -C so as to match vibrational measurements. Accordingly, let ω and ξ denote measured values of an angular eigenfrequency ω and its corresponding eigendisplacement u ⋆ . In ECE-based updating techniques, the distributions ρ ⋆ , C ⋆ are sought so as to minimize the error criterion

J (C, ρ) = measured ω min (v,Γ,s)∈X H ω(v, γ, s, C, ρ) (101)
where the summation is taken over all eigenmodal measurements (ω, ξ) and the product set X of admissible fields is defined by

X = (v, γ, s) | v ∈ C(0, S u ), γ ∈ C(0, S u ), s ∈ D(0, γ, S p ) (102)
with C(0, S u ) is again defined according to [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] while D(0, γ, S p ) is the set of stresses dynamically admissible with respect to γ, i.e.

D(φ, γ, Σ) = {σ | div σ + γ = 0 in Ω, σ.n = φ on Σ} (103)
Moreover, the modified error in constitutive equation H ω(v, γ, s, C, ρ) is defined by

H ω(v, γ, s, C, ρ) = α 2 E(v, s, C) + 1 -α 2 I ω(v, γ, ρ) + γ 2 D a(v-ξ, v-ξ) dV ( 104
)
where a(z, z) is an arbitrarily chosen energy-type bilinear form defined on the displacement measurement area, E(v, s, C) is the error in constitutive equation defined by ( 78), I ω is the error in dynamical constitutive equation defined by

I ω(v, γ, ρ) = Ω 1 ρω 2 |γ + ρω 2 v| 2 dV, (105)
γ > 0 and 0 ≤ α ≤ 1 are weighting coefficients. It can be shown [START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF][START_REF] Ladeveze | Parametric correction of finite element models using modal tests[END_REF] that J (C, ρ) is zero if and only if C = C ⋆ and ρ = ρ ⋆ . Thus, a primal-dual formulation similar to (99) could conceivably be set up, with γ then becoming a penalty parameter allowed to assume arbitrarily large values. The minimum in ( 101) is restricted to the set of admissible fields X , eq. ( 102). Accordingly, introduce the Lagrangian L:

L = H ω(v, γ, s, C, ρ) + Ω (s : ε[w] + γ.w) dV (106)
which incorporates the dynamical constraint s ∈ D(0, γ, S p ), equation ( 103), in weak form and where the test function w ∈ C, i.e. the Lagrange multiplier field, belongs to C(S u , 0). The first-order variation of L then reads:

δL = L ,v • δv + L ,s • δs + L ,Ga • δγ + L ,C • δC + L ,ρ δρ (107)
where, for each measured frequency:

L ,s • δs = α Ω (C -1 : s -ε[v] + 1 α ε[w]) : δs dV L ,Ga • δγ = (1 -α) Ω 1 ρω 2 γ + v + 1 1 -α w .δγ dV L ,v • δv = α Ω (C : ε[v] -s) : ε[δv] dV + (1 -α) Ω [γ + ρω 2 v].δv dV + γ D a(v-ξ, δv) dV and L ,C • δC = α Ω ε[v] : δC : ε[v] -(C -1 : s) : δC : (C -1 : s) dV (108)
L ,ρ δρ = (1 -α) Ω ω2 |v| 2 - 1 ρ 2 ω2 |γ| 2 δρ dV (109)
To solve the partial minimization problem involved in the definition (101) of J (C, ρ), one has to enforce the first-order optimality conditions L ,s = L ,Γ = L ,v = 0 with the help of the above equations. From the first two conditions, the minimizers s and γ are given by

s = C : ε v - 1 α w γ = -ρω 2 v 1 (1 -α)
.w [START_REF] Marin | Boundary element solution for the Cauchy problem in linear elasticity using singular value decomposition[END_REF] in terms of v and the Lagrange multiplier w. Then, on enforcing L ,v = 0 together with the constraint ( 103), (v, w) are found to be governed by the coupled equations:

Ω α ε[w] : C : ε[v ′ ] -ρω 2 w.v ′ dV + γ D a(v-ξ, v ′ ) dV = 0 ∀v ′ ∈ C(0, S u ) (111) Ω ε v- 1 α w : C : ε[w ′ ] -ρω 2 v+ 1 (1 -α) w .w ′ dV = ∀w ′ ∈ C(0, S u ) (112)
As a result, the cost function ( 101) is found to have the explicit expression

J (C, ρ) = γ 2 D a(v-ξ, v-ξ) dV + 1 2 Ω 1 α ε[w] : C : ε[w] + ρω 2 1-α |w| 2 dV ( 113
)
in terms of the solution (v, w) to ( 111)- [START_REF] Marklein | Linear and nonlinear inversion algorithms applied in nondestructive evaluation[END_REF]. In addition, the sensitivity of J (C, ρ) to perturbations δC and δρ is obtained directly from (107), upon substituting v, w, γ, s into ( 108) and ( 109), as

δJ = Ω ε 2v - 1 α w : δC : ε[w] -ω2 2v + 1 1-α
w .w δρ dV (114)

Prior localization of defects

The distributions ε[w] : C : ε[w] and ρ|w| 2 associated with estimators of modeling errors based on the ECE have the ability to provide an a priori estimation of the geometrical support of localized modelling errors prior to inversion. This useful feature is supported by many numerical experiments [START_REF] Bonnet | On the estimation of the geometrical support of modelling defects using the distributed error in constitutive equation[END_REF][START_REF] Chouaki | An updating of structural dynamic model with damping[END_REF][START_REF] Deraemaeker | Reduced bases for model updating in structural dynamics based on constitutive relation error[END_REF][START_REF] Ladeveze | Parametric correction of finite element models using modal tests[END_REF][START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF], exploiting either synthetic or real experimental data, and similar results are available for static problems [START_REF] Bui | Spatial localization of the error of constitutive law for the identification of defects in elastic bodies[END_REF]. This observation can, again, be given a partial formal explanation, as follows. Choose α = 1 and γ = -ρω 2 v (i.e. the error in elastic constitutive relation), so that the cost function H ω involves the error in elastic constitutive relation E but not the error in dynamical equation I ω, and assume an ideal situation where the measured mode ξ is known exactly and on the whole structure Ω. Then, in the limiting case γ → ∞, the system (111)-( 112) reduces to

v = ξ and Ω ε[w] : C : ε[w ′ ] dV = Ω ε[v] : C : ε[w ′ ] -ρω 2 v.w ′ dV (∀w ′ ∈ V)
Since v = ξ is the measured eigenmode corresponding to the eigenvalue ω and the unknown stiffness and mass distributions, C ⋆ and ρ ⋆ , the second equation above becomes

Ω ε[w] : C : ε[w ′ ] dV = - Ω ε[ξ] : ∆C : ε[w ′ ] -∆ρω 2 ξ.w ′ dV (∀w ′ ∈ V) (115)
and w therefore appears to be governed by a static equilibrium problem associated with the body force distribution ∆ρω 2 ξ and the initial strain distribution C -1 : ∆C : ε[ξ]. Denoting by D ⊂ Ω the geometrical support of the corrections ∆C and ∆ρ, the solution w of (115) has the representation

w( x) = D ε[ξ](x) : ∆C : ∇ x G( x, x) -∆ρω 2 ξ(x).G( x, x) dV x (116)
where G( x, x) is the elastostatic Green's tensor for elastic moduli C and suitable homogeneous boundary conditions on ∂Ω. From this result, one can show that the ECE density ε[w] : C : ε[w] evaluated at x away from the defect support D behaves like r -6 for stiffness defects (i.e. ∆C = 0) and r -4 for mass defects (i.e. ∆ρ = 0), where r denotes characteristic distance from x to D (the main argument used is the fact that, for threedimensional bodies, G( x, x) ∼ |x -x| -1 ). This in particular corroborates the fact, generally observed in numerical experiments, that the geometrical localization using ECE density is better for stiffness defects than for mass defects.

The analysis for the two-term ECE (0 < α < 1) is similar, with the difference that the Green's tensor now solves

div [C : ∇G( x, •)] -ρω 2 G( x, •) + δ xI = 0 (117)
One has again G( x, x) ∼ |x -x| -1 , however, so that the previous qualitative conclusions hold again. The a priori localization properties are in practice good even with few and/or inexact eigendisplacement measurements. A general quantitative analysis of this fact for imperfect or partial data, and of the influence of the number of measured eigenfrequencies, which would provide a better understanding of the still largely empirical localization properties, is not currently available. Two illustrative example of the a priori localization capabilities of the ECE density are now presented.

Example 1 (synthetic data).

A finite element model of a cooling tower (a component of a power plant) in which a crack is simulated is considered. The eigendisplacement field ξ associated with a vibrational mode has been computed, together with the distributed ECE, under the assumption that a complete measurement of ξ is available. Figure 6 presents the FE mesh and (through the distribution of Young moduli) the simulated defect, while the distribution of ECE over the mesh is depicted on Figure 7. Thus, it is apparent that in the most favorable situation of a complete and exact measurement of the field ξ (disregarding the FE discretization error), the distribution of ECE provides a very accurate estimation of the geometrical support of the defect.

Example 2 (experimental data

). An experimental setup, made for a study conducted at EDF, the French electric power company, consists of an elastic thin shell structure as described in figure 8, with Young modulus E = 1, 97 10 11 N/m 2 , Poisson ratio ν = 0, 3, specific mass ρ = 8000 kg/m 3 . Eigendisplacements in the direction normal to the shell are measured at 48 sensor locations (4 rings of 12 equally spaced sensors located on the cylindrical part of the structure, respectively 22,5cm, 45cm, 67,5cm and 90cm above the disk-shaped bottom). The test structure rests on three elastic supports located at (r = 0.3 m, θ = 0, 2π/3, 4π/3), whose stiffness is estimated at k = 10 7 N/m along all coordinate directions. Experimental eigenfrequency and eigendisplacement data has been recorded for two artificial "defects" created by matter removal (rectangular holes) in the cylindrical part. A finite element model of the structure has been set up using the general-purpose code CAST3M, in which the cylindrical part is made of 4 identical rings 24 eight-noded shell elements. Each node supports six unknowns (3 displacements and 3 rotations), for a total of 3822 unknowns. The sensor locations coincide with mesh nodes.

VAL -ISO >-8.80E+03 < 6.67E+04 -8.21E+03 -4.67E+03 -1.13E+03 2.41E+03 5.95E+03 9.49E+03 1.30E+04 1.66E+04 2.01E+04 2.36E+04
The figure 9 depicts the values of the ECE integrated over each finite element (normalized to max=1) obtained for the "small" and "large" defects (respectively located on elements 61-62 and 61-62-63-64), using experimental data for 13 eigenmodes (i.e. 13 frequencies and 13 × 48 displacements). The highest values of the distributed ECE are seen to correspond to the defects and their neighbourhood (in a FEM mesh sense). One sees that, in spite of the discrete nature of the measurements and of the inevitable experimental errors, the distribution of ECE still points correctly to a region of the structure that contains the actual defect. . ECE distributions (normalized to max=1) for the "small" (elements 61-62, left) and "large" (elements 61-62-63-64, right) defects, using real data for 13 measured eigenmodes.

Geometrical sensitivity techniques for defect identification

Defect identificaton problems are often formulated in terms of the minimization of a cost function featuring the experimental data and (possibly) prior information. Such cost functions are nonconvex and exhibit local minima. Despite that fact, traditional minimization methods are usually preferred to global search techniques such as evolutionary algorithms, the latter being in most cases infeasible due to the prohibitive computational cost of large numbers of direct elastic scattering solutions. To perform optimally, gradient-based optimization techniques are used in conjunction with shape sensitivity formulations. Such formulations are now well-established, especially for dealing with perturbations of smooth shapes, see e.g. [START_REF] Dems | On a class of conservation rules associated with sensitivity analysis in linear elasticity[END_REF][START_REF] Kirsch | The domain derivative and two applications in inverse scattering theory[END_REF][START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF][START_REF] Sokolowski | Introduction to shape optimization[END_REF]. A typical derivation, using the adjoint solution method, is provided in section 6.1 for establishing cavity shape sensitivity formulae expressed in either domain integral form (suitable for domain discretization techniques such as the FEM) or boundary integral form (suitable for the BEM). When integral equation techniques are applied to crack identification problems, the previous result is not applicable and a specific sensitivity formula, established in section 6.1, is needed. Still, the stand-alone use of gradient-type minimization for such purposes is not satisfactory for its success is strongly dependent on a reliable prior information about the geometry of the hidden object. This has prompted the development of alternative techniques, which may be used in isolation or as a preliminary step for choosing adequate initial guesses in standard optimization schemes. One such technique, under current development and presented in Section 6.4, is based on the concept of topological derivative and aims primarily at providing a preliminary elastic-wave imaging of solid bodies, allowing a selection of suitable initial guesses in terms of location and size of defects. Another is the linear sampling, which is not addressed in this article (see [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF] for a recent investigation in elastodynamics, and the references therein). For completeness, it is also important to mention techniques based on level set methods, not addressed in this article, whose implementation also involve shape sensitivity concepts [START_REF] Ben Ameur | Level set methods for geometric inverse problems in elasticity[END_REF][START_REF] Burger | Levenberg-Marquardt level set methods for inverse obstacle problems[END_REF][START_REF] Ramananjaona | Shape reconstruction of buried obstacles by controlled evolution of a level set: from a min-max formulation to numerical experimentation[END_REF].

In this section, geometrical sensitivity techniques are presented in connection with the problem of determining the shape and position of an unknown object (cavity or crack) embedded in the elastic body from elastodynamic experimental data. Time-harmonic conditions are assumed for simplicity, although a treatment in the time domain could have been presented as well.

A defect bounded by the surface Γ is embedded in an elastic body Ω ∈ R 3 externally bounded by the closed surface S (i.e. Ω denotes the defect-free reference body). Let B and Ω Γ = Ω \ (B ∪ Γ) denote the interior of the defect(s) and domain occupied by the flawed body, respectively. The defect is taken to be either a cavity (i.e. Γ is a closed surface and |B| = 0) or a crack (i.e. B = ∅ and Γ is an open surface across which displacement jumps are expected); multiple defects can of course be represented with a non-connected surface Γ. Under these assumptions, the direct problem for the displacement u is of the form

(A C + ρω 2 )u = 0 in Ω Γ ; u = ξ on S u , p[u] = φ on S p , p[u] = 0 on Γ (118)
Equivalently, it is governed by the weak formulation

Ω Γ [σ[u] : ∇w -ρω 2 u.w] dV - Sp φ.w dS - Su [(u -ξ).q + p.w] dS = 0 (119)
where (w, q) denote unconstrained trial displacements and tractions. The problem of determining the shape and position of an unknown defect using experimental data, e.g. ultrasonic measurements, is considered. The lack of information about Γ is compensated by overspecified boundary data on S. Here, the available experimental data is assumed of the form

u(x) = ξ(x) (x on S m u ⊆ S u ) p(x) = φ(x) (x on S m p ⊆ S p ) (120)
where the measurement surfaces S m u and S m p do not necessarily cover the whole of S; in particular, pointwise measurements can be represented in the above fashion with the aid of appropriate Dirac distributions.

Here, the formulation of the defect identification in terms of the minimization of a cost function is considered. The simplest choice for the latter, as usual, is a least-squared misfit function, e.g.

J (Γ) = J(u Γ , p Γ , Γ) = 1 2 S m p (ξ -u Γ ) 2 dS + 1 2 S m u (φ -p Γ ) 2 dS (121)
where (u Γ , p Γ ) refer to the solution of the direct problem [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF] for a given Γ. For the sake of generality, generic cost functions defined by

J(u, p, Γ) = Sp ϕ u (u, x) dS + Su ϕ p (p, x) dS + Γ ψ(x) dS (122)
which obviously generalize [START_REF] Oden | Variational methods in theoretical mechanics[END_REF], are considered in this section (the last term, explicit in terms of Γ, may be used to express prior information).

In the context of elastic-wave imaging, the computational cost of each direct solution is very high for realistic three-dimensional situations. This precludes the use of global search techniques like genetic algorithms which entail a large number of forward simulations. The need to keep the number of direct computations as low as possible suggests instead to stick with classical gradient-based optimization algorithm. The minimization of J with respect to Γ using such methods needs in turn, for efficiency, the ability to evaluate the gradient of the functional J with respect to perturbations of the shape of Γ, in addition to J (Γ) itself.

Adjoint formulations for cavity shape sensitivity

Any small perturbation of Γ can be described by means of a domain transformation which does not affect the external boundary S, i.e. of the form x b = x + bθ(x) where b > 0 is a time-like parameter and the transformation velocity field θ(x) is such that θ = 0 on S.

Denoting by

⋆ f = f ,b + ∇f.
θ the Lagrangian derivative of some field variable f , the derivative at b = 0 of generic integrals over a domain V or a surface Σ have the well-known form

d db V f dV = V ( ⋆ f + f div θ) dV d db Σ f dS = Σ ( ⋆ f + f div S θ) dS (123)
where div S θ = div θn.∇θ.n is the surface divergence of θ (n: outward unit normal vector). Also, recall that (∇u) ⋆ = ∇ ⋆ u -∇u.∇θ. To establish shape sensitivity formulae for cost functions of the form [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], one may either use an optimal control approach whereby a Lagrangian is introduced so as to incorporate the direct elastodynamic problem as an equality constraint, or exploit a reciprocity identity together with a direct differentiation of [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] with respect to the defect shape. The latter approach is chosen here. The sensitivity of the cost functions to perturbations of the cavity shape is obtained in two forms . One involves domain integrals and is therefore wellsuited to domain discretization method such as the finite element method (section 6.1.1), while the other features boundary integrals and is better suited to solution techniques based on boundary integral equation methods, which are frequently used for this kind of inverse problem (section 6.1.2). In ( 127), (w, q) are components of an arbitrary elastodynamic states without body forces. The adjoint solution is defined as the particular elastodynamic state that fulfills the well-posed set of boundary conditions

w = ∂ϕ p ∂p (x on S u ) , q = - ∂ϕ u ∂u (x on S p ) , q = 0 (x on Γ) (128)
On substituting these boundary conditions into [START_REF] Robertson | Boundary identifiability of residual stress via the Dirichlet to Neumann map[END_REF] and adding the resulting equation to equation [START_REF] Ramananjaona | Shape reconstruction of buried obstacles by controlled evolution of a level set: from a min-max formulation to numerical experimentation[END_REF], one obtains the following sensitivity formula:

dJ db = Ω Γ σ[w Γ ] : ∇u Γ -ρω 2 w Γ .u Γ div θ -σ[w Γ ].∇u Γ + σ[u Γ ].∇w Γ : ∇θ dV + Ω Γ ∇u Γ : (∇C.θ) : ∇w Γ dV + Γ (∇ψ.θ + ψdiv S θ) dS (129)
The sensitivity formula (129) provides a very efficient means for computing derivatives of the cost function [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], because (i) as usual in adjoint solution methods, only one adjoint solution is needed, irrespective of the number of geometrical parameters which are used in practice for the representation of the unknown defect shape,and (ii) the adjoint solution is governed by the same field equations as the direct solution, allowing re-use of the (built and factored) discrete matrix operator involved in the direct solution. Related sensitivity formulae are given in e.g. [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF][START_REF] Feij Óo | An application of shape optimization in the solution of inverse acoustic scattering problems[END_REF].

Sensitivity formula in boundary integral form.

A sensitivity formula in boundary integral form can in fact be established from [START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF] by observing that the following identity holds for the direct and adjoint solutions (or, more generally, for any pair of elastodynamic states without body forces):

σ[w] : ∇u -ρω 2 w.u div θ -σ[w].∇u + σ[u].∇w : ∇θ -∇u Γ : (∇C.θ) : ∇w = div σ[w] : ∇u -ρω 2 w.u θ -σ[w].∇u + σ[u].∇w .θ (130)
Application of identity [START_REF] Santos | On the solution of an inverse scattering problem in seismic whiledrilling technology[END_REF], together with the divergence formula and the traction-free conditions p Γ = q Γ = 0 on Γ, to (129) leads to the following sensitivity formula, equivalent to [START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF] but expressed in terms of boundary integrals:

dJ db = Γ [σ[w Γ ] : ∇u Γ -ρω 2 w Γ .u Γ ]θ n dS + Γ (∇ψ.θ + ψdiv S θ) dS (131)
To ensure that formula [START_REF] Schmerr | Ultrasonic flaw sizing inverse problems[END_REF] can actually be evaluated using only the boundary traces of the direct and adjoint solutions, the bilinear form σ[w] : ∇u must be expressed in terms of the surface gradients ∇ S u, ∇ S w, taking p = q = 0 into account in the process:

σ[u Γ ] : ∇w Γ = µ 2ν 1 -ν div S u Γ div S w Γ + 1 2 (∇ S u Γ + ∇ T S u Γ ) : (∇ S w Γ + ∇ T S w Γ ) -(n.∇ S u Γ ).(n.∇ S w Γ ) (132)
The evaluation of shape sensitivities using ( 131)-( 132) is very useful in cases where the direct and adjoint elastodynamic solutions are computed by means of the boundary element method. Such implementation is presented for three-dimensional elastodynamics in [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF][START_REF] Nintcheu Fata | A Computational Basis for Elastodynamic Cavity Identification in a Semi-Infinite Solid[END_REF]

Adjoint formulations for crack shape sensitivity

Defect identification is often concerned with finding cracks. Crack identification by means of reciprocity gap functionals was presented in section 49. If one instead considers formulating the problem in terms of minimizing an objective function or solving the set of observation equations by a Newton-type method, then formulae for crack shape sensitivity are useful. Ideealized cracks are geometrically singular defects. As a result, the elastodynamic strains and stresses are singular as well, with an asymptotic behaviour of O(d -1/2 ) in the limit d → 0 (where d is the distance to the crack front ∂Γ). One is therefore led to enquire whether the sensitivity formulae ( 129) and ( 131) remain valid for crack identification problems and, if not, to establish valid substitute formulae.

The density of the domain integral in (129) has a O(d -1 ) singularity, which is integrable. Besides, an examination of the steps leading to [START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF] reveals that all domain integrals involved in the derivation are integrable. In particular, ∇u and ∇ ⋆ u have the same O(d -1/2 ) crack front singularity, whereas the partial derivative ∂(∇u)/∂b has a O(d -3/2 ). It was therefore essential to express the derivative weak formulation [START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF] in terms of the Lagrangian derivative of the direct solution. From these considerations, it can be concluded that the sensitivity formula in domain integral form, equation [START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF], remains valid when the cavity is collapsed onto a crack.

In contrast, formula (131) certainly does not hold when Γ is a crack. For example, one has θ n (x) = 0, and hence dJ /d = 0 if ψ = 0, for any in-plane stretching transformation of the crack (whereas for a cavity θ n (x) = 0 implies that its shape does not change). This contradicts e.g. well-known results of fracture mechanics [START_REF] Freund | Dynamical fracture mechanics[END_REF] showing that the mechanical potential energy of a cracked solid is affected by in-plane crack extensions. Besides, identity [START_REF] Santos | On the solution of an inverse scattering problem in seismic whiledrilling technology[END_REF] is not applicable for cracks because the divergence on the right-hand side has a O(d -2 ) singularity, which is not integrable.

To find a sensitivity formula in boundary integral form applicable to cracks, one has to resort to a limiting process, outlined now. The body Ω is partitioned into Ω = Ω ε ∪ (D ε \ Γ), where D ε = {x, dist(x, ∂Γ) ≤ ε} for some sufficiently small ε > 0 is a tubular neighbourhood of the crack front ∂Γ bounded by the tubular surface Σ ε and Ω ε = Ω \ Dε . On introducing this splitting into Eq. ( 129), applying the divergence formula for the contribution of Ω ε and invoking the fraction-free boundary condition assumed to hold on Γ for both the

0 (K[u Γ ] I K[w Γ ] I + K[u Γ ] II K[w Γ ] II ]) dt ds (140)
where, in particular, the SIFs K[u Γ ] I etc. are now functions of both s and t. Note that the last integral in [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], which is unaffected by whether time-harmonic or transient conditions are assumed, has been omitted in [START_REF] Tarantola | Inverse problem theory[END_REF] for the sake of brevity. To take it into account in a time-domain cost function of the form [START_REF] Tarantola | Inverse problem theory[END_REF], one has simply to append the last term of formulae ( 129), ( 131) and ( 135) to ( 138), ( 139) and [START_REF] Weglein | Inverse scattering series and seismic exploration[END_REF], respectively. The sensitivity formula [START_REF] Weglein | Inverse scattering series and seismic exploration[END_REF] is applied in [START_REF] Bonnet | Sensitivity analysis for shape perturbation of cavity or internal crack using BIE and adjoint variable approach[END_REF][START_REF] Burczy Ński | Sensitivity analysis and identification of material defects in dynamical systems[END_REF] to 2D crack identification problems formulated by means of a time-domain BEM approach, with the 2D frequency-domain situation is considered in [START_REF] Rus | Crack shape sensitivity by the adjoint variable method using a boundary-only integral[END_REF].

Topological derivative

Cost functions such as [START_REF] Oden | Variational methods in theoretical mechanics[END_REF] or [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] are highly non-convex. Thus, even using efficient sensitivity computation techniques such as those outlined thus far in this section, gradientbased optimization techniques for defect identification thus have obvious drawbacks: one run provides just one local minimum, which may depend on the (blind) choice of an initial guess. These pitfalls can be circumvented to some extent by prior information, which, if in good agreement with the physical reality, may lead to improvements for both the choice of initial guess and the cost function properties.

Another approach consists in defining computational procedures that allow a preliminary probing of the medium, thus performing a global search. At least two directions are currently investigated for that purpose. One is the linear sampling method. Initially developed in the context of acoustic and electromagnetic wave imaging of unbounded media [START_REF] Colton | The linear sampling method in inverse electromagnetic scattering theory[END_REF][START_REF] Colton | A simple method for solving inverse scattering problems in the resonance region[END_REF], it has also been extended to elastic-wave imaging [START_REF] Charalambopoulos | The linear sampling method for the transmission problem in three-dimensional linear elasticity[END_REF][START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF]. In addition, a version of the linear sampling method applicable to elastodynamic near-field data has very recently been proposed [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF]. Another direction, to which this section is devoted, revolves around the preliminary computation of the topological derivative.

The topological derivative of J , T (x o ), synthesizes the sensitivity of J with respect to the creation of an infinitesimal cavity at a chosen location x o inside the reference, i.e. cavityfree counterpart of the probed body. The concept of topological derivative first appeared in [START_REF] Eschenauer | Bubble method for topology and shape optimization of structures[END_REF] and [START_REF] Schumacher | Topologieoptimierung von Bauteilstrukturen unter Verwendung von Lochpositionierungskriterien[END_REF] in the context of shape optimization of mechanical structures, wherein the spatial distribution of T (x o ) was used as a criterion for the removal of "excess" material through regions where T < 0. Recently, its rigorous mathematical formulation has been established within the framework of elastostatic problems and Laplace equation [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF].

Although initially defined in connection with optimization problerms, the concept of topological derivative is also applicable to defect identification problems. Here, its usefulness as a means of facilitating a subsequent minimization-based solution by providing a rational basis for selecting a reliable initial guess in terms of its topology, approximate size and location is considered. In [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], this idea was considered in the context of inverse elastic scattering pertaining to semi-infinite and infinite domains, where the availability of suitable fundamental solutions made it possible to establish explicit expressions for T (x o ), while developments along similar lines for 2D elastostatics are presented in [START_REF] Gallego | Identification of cracks and cavities using the topological sensitivity boundsary integral equation[END_REF][START_REF] Jackowska-Strumillo | Topological optimization and inverse problems[END_REF].

It is instructive to compare the topological derivative concept, as applied to the direct and inverse scattering of small objects by waves, to other asymptotic approaches. In cases featuring only one characteristic length (one scatterer embedded in an unbounded medium an illuminated by a plane wave), the topological derivative approach essentially provides (up to a scaling factor) the lowest-order moment of the normalized scattering amplitude in the theory of low-frequency direct and inverse scattering [START_REF] Charalambopoulos | Inverse scattering via low-frequency moments[END_REF][START_REF] Dassios | Energy functionals in scattering theory and inversion of low frequency moments[END_REF][START_REF] Dassios | Low frequency scattering[END_REF]. However, in situations where other characteristic lengths (size of a finite body, radius of curvature of wave fronts) are involved in addition to the vanishing size of the scatterer, the topological derivative approach differs from the low-frequency approach. Also relevant is the systematic method for deriving asymptotic expansions of the relevant field variable, and its application to the identification of small inhomogeneities, proposed in [START_REF] Alves | Boundary integral formulae for the reconstruction of imperfections of small diameter in an elastic medium[END_REF][START_REF] Ammari | Reconstruction of small inhomogeneities from boundary measurements[END_REF][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]. Asymptotic expansions of cost functions, and in particular the value of its topological derivative, can then be obtained as a post-treatment, in a manner reminiscent of the direct approach in sensitivity analysis.

In what follows, a procedure for the computation of the topological derivative applicable to arbitrary elastic bodies is presented. It is based on the use of an adjoint solution, in order to avoid the computation of the field sensitivities and thereby to yield the topological derivative in a computationally optimal fashion. 6.4.1. Topological derivative for elastodynamic scattering. To this end, let B ε (x o ) = x o + εB, where B ⊂ R 3 is a fixed bounded open set of boundary S and volume |B| which contains the origin, define the region of space occupied by a cavity of (small) size ε > 0 containing a fixed sampling point x o . Following [START_REF] Garreau | The topological asymptotic for PDE systems: the elasticity case[END_REF][START_REF] Sokolowski | On the topological derivative in shape optimization[END_REF], one is interested in the asymptotic behavior of J (Ω ε ) for infinitesimal ε > 0, where Ω ε = Ω \ B ε (x o ). With reference to this limiting behavior, the topological derivative T (x o ) of the cost functional J (Ω) at x o for a cavity-free body is defined through the expansion:

J (Ω ε ) = J (Ω) + δ(ε)|B|T (x o ) + o(δ(ε)) (ε ≪ Diam(Ω), B ε (x o ) ⊂ Ω) (141)
where the function δ(ε), such that lim ε→0 δ(ε) = 0, characterizes the asymptotic behaviour of J (Ω ε ) and will be specified as a result of the analysis to follow. One may note that this definition is not restricted to spherical infinitesimal cavities (for which B is the unit ball, S the unit sphere and |B| = 4π/3). In general, the value T (x o ) is expected to depend on the shape of B.

With reference to [START_REF] Wexler | An impedance computed tomography algorithm and system for ground water and hazardous waste imaging[END_REF], the evaluation of J (Ω ε ) requires the knowledge of the elastodynamic solution u ε to the direct problem [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF] with B replaced by B ε ≡ B ε (x o ). To this end, it is convenient to decompose the total displacement field u ε as

u ε = u + ũε (142)
where ũε denotes the scattered field and u is the free field defined as the response of the defect-free (reference) solid Ω due to the given excitation (ξ, φ).

To investigate the asymptotic behaviour of J (Ω ε ) defined by [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF], one may start with the expansion

J (Ω ε ) = J (Ω) + Re Sp ∂ϕ u ∂u . ũε dS + Su ∂ϕ p ∂p .p ε dS + o(| ũε | Sp , |p ε | Su ) (143)
Again, the reciprocity identity will provide a very simple approach to choose a suitable adjoint state and establish a formula for the topological derivative T (x o ). Following the same pattern as in section 6.1, let (w, q) denote the trace on ∂Ω ε of an arbitrary elastodynamic state without body forces defined on Ω. Then, the elastodynamic reciprocity theorem applied to the boundary traces ( ũε , p ε ) and (w, q) leads to the identity Su p ε .w dS -Sp ũε .q dS -Γε [ ũε .q + p.w] dS = 0 (144) having used the fact that u ε = 0 on S u , p ε = 0 on S p and p ε = -p on Γ ε by virtue of the decomposition [START_REF] Zimmermann | Structural damage detection using a minimum rank update theory[END_REF]. Now, the adjoint solution is defined as the elastodynamic state in Ω that fulfills the well-posed set of boundary conditions

w = - ∂ϕ p ∂p (x on S u ) , q = ∂ϕ u ∂u (x on S p ) (145)
Substituting these boundary conditions into (144) and adding the resulting equation to (143) leads to

J (Ω ε ) = J (Ω) -Re Γε [p.w + ũε .q] dS + o(| ũε | Sp , |p ε | Su ) (146)
To evaluate the asymptotic behavior of the integral in the above formula, the virtual work principle applied for u in the domain B ε with w as trial function yields

- Γε p.w dS = Bε [σ : ∇w -ω 2 u.w] dV = ε 3 |B|[σ : ∇w -ω 2 u.w](x o ) + o(ε 3 ) (147)
(where the minus sign in front of the boundary integral appears because the surface integral over Γ ε involves the inward unit normal to Γ ε ). Besides, the scattered displacement on Γ ε is known [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] to have the asymptotic behaviour

ũε (εxb) = εσ kℓ (x o )V kℓ (xb) + o(ε) (xb ∈ S, x = εxb ∈ Γ ε ) (148)
where the V kℓ (z) solve the six canonical elastostatic exterior problems

A C V kℓ = 0 (in R 3 \ B) , σ[V kℓ ].n = - 1 2 (n k e ℓ + n ℓ e k ) (on S)
Finally, from the asymptotic behaviour (148), it is easy to prove by means of an integral representation formula [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Kupradze | Dynamical problems in elasticity[END_REF] that

| ũε | Sp = O(ε 3 ) , |p ε | Su = O(ε 3 ) (149)
On gathering (147), ( 148) and ( 149) into (146), the precise expressions for the topological derivative T (x o ) and the asymptotic behaviour δ(ε) introduced in expansion (141) are obtained as:

T (x o ) = Re σ[w] : D : ∇u -ρω 2 w.u (x o ) δ(ε) = ε 3 (150)
where the constant fourth-order tensor D is defined by

D ijkℓ = 1 2 (δ ik δ jℓ + δ iℓ δ jk ) - 1 B S [λV aa i δ kℓ + 2µV kℓ i ](xb)n j (xb) dS
The expression (150) of the topological derivative is valid for any shape of the small cavity B ε . For spherical cavities (i.e. when B ε is a ball of radius ε and S is the unit sphere), the canonical solutions can be evaluated analytically, and the tensor D is found as a result to have the closed-form expression

D ikjℓ = 3(1 -ν) 2(7 -5ν) 5(δ ik δ jℓ + δ iℓ δ jk ) + 5ν -1 1 -2ν δ ij δ kℓ
Numerical experiments based on the topological derivative in elastodynamics recently appeared in connnection with identification of cavities [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF][START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF] and of penetrable elastic inclusions [START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF].

Topological derivative for acoustic scattering.

A formula similar to (150) can be established for the linear acoustic case [START_REF] Guzina | Small-obstacle asymptotics and scalar wave imaging[END_REF]. Consider the scattering of acoustic waves by a small penetrable inclusion B ε (x o ) embedded in a reference medium characterized by the constant wavenumber k and the mass density ρ, the inclusion material being characterized by the wavenumber k/γ and the mass density ρ/β (with γ, β > 0). The total field u (e.g. the acoustic pressure) is then governed by

(∆ + k 2 )u = 0 (in Ω ε ), (∆ + γ 2 k 2 )u = 0 (in B ε ) (151)
together with the standard transmission conditions on the interface Γ ε (x o ) and boundary conditions such that values of u and of p = ∂u/∂n are prescribed on S u and S p = ∂Ω \ S u , respectively. Equivalently, u is also governed by a generalization of the Lippman-Schwinger integral equation [START_REF] Martin | Acoustic scattering by inhomogeneous obstacles[END_REF]. Consider cost functions of the form

J (ε) = J(u ε , p ε ) = Sp ϕ u (u ε , x) dS + Su ϕ p (p ε , x) dS (152)
where (u ε , p ε ) solve the transmission problem for an assumed inclusion B ε with constitutive properties defined by (γ, β). One then finds

J (ε) = J (0) + ε 3 |B|T (x o ) + o(ε 3 ) (153)
where the topological derivative T (x o ) is given this time by

In (154), u is the acoustic free-field while the adjoint field w solves (∆ + k 2 )w = 0 in Ω, w = -∂ϕ p ∂p on S u , q = ∂ϕ u ∂u on S p and the second-order tensor D(B, β, γ) is known for any inclusion shape B and constitutive parameters (β, γ) [START_REF] Guzina | Small-obstacle asymptotics and scalar wave imaging[END_REF]. For the simplest case where B is the unit sphere, one has

D ij (B, β, γ) = 3(1 -β) 2 + β δ ij
The limiting situation β = 0 in (154) yields the expression of the topological derivative for the case of a hard (i.e. rigid) obstacle of vanishing size ε.

Numerical examples in linear acoustics.

To illustrate the foregoing notions, sample results from numerical experiments, for frequency-domain acoustics, are now presented. Similar experiments are currently being conducted for the preliminary imaging of cavities and elastic inclusions using elastodynamic data. The topological derivative T (x o ) for the cost function J (Γ, β, γ) and wavenumbers k = 1 and k = 4 in the reference medium is displayed in figures 11 and 12, respectively. Again, the values obtained for T (x o ) in a series of horizontal planes are consistent with the actual location of the "true" inclusion, despite the fact that the latter is of finite size while the asymptotic formula (153) only holds in the limit ε → 0.

Conclusion and further developments

In this article, several types of inverse problems arising in the linear theory of elasticity have been considered, revolving around the identification of material parameters, distributions of elastic moduli, and geometrical objects such as cracks and inclusions. The underlying thread of a large part of the discussion was the fact that useful tools for the formulation, analysis and solution of inverse problems arising in linear elasticity, namely the reciprocity gap and the error in constitutive equation, stem from variational and virtual work principles, i.e. fundamental principles governing the mechanics of deformable solid continua. In addition, the virtual work principle appears to be instrumental in establishing computationally efficient formulae for parameter or geometrical sensitivity, in connection to the adjoint solution method. Sensitivity formulae have been presented for various situations, especially in connection with contact mechanics, cavity and crack shape perturbations, thus enriching an already extensive repertoire of such results. Finally, the concept of topological derivative and its implementation for the identification of cavities or inclusions have been expounded. Many of the techniques described in this article are especially well suited to continuous data, i.e. in connection with experimental techniques allowing to measure fields (in practice, discrete data sampled on a very fine grid and over an extensive portion of the body or its boundary). In particular, recent optical techniques now allow to measure displacements fields or in-plane strain fields on the boundary, and microtomographic techniques allow to measure displacements inside a body. Temperature field data provided by infrared thermography can also be useful due to the thermomechanical coupling present in most solids.

Real-world mechanical systems feature complex geometries and limited data, whereas rigorous analyses of fundamental issues such as existence or uniqueness of solutions, or convergence properties of inversion algorithms are usually available only for idealized situations. Advances in the complexity of mechanical systems for which such fundamental results hold are obviously desirable, but it is expected that many real-world complex engineering systems will remain out of the reach of such analyses. Nevertheless, most of the topics presented in this article were selected on the basis of their ability to be implemented within the framework of the general-purpose numerical solution techniques used in solid mechanics, namely the finite element method and the boundary element method.

Some directions for further work directly related to topics presented in this article include (i) a more systematic development and testing of the reciprocity gap approach as a method for the identification of cracks, (ii) the study of convexity properties of functionals based on the error in constitutive equations, and develop similar functional in connection with the identification of nonlinear constitutive properties and (iii) the development of topological sensitivity techniques for time-domain formulations and in refined forms (in particular based on higher-order expansions with respect to defect size) and their integration in defect identification strategies.

Figure 1 .

 1 Figure 1. Identification of the convex hull of a set of two straight cracks using the instantaneous reciprocity gap.

  For a given choice of the parameter vector b within some admissible set B ⊆ R B , the computed displacement field u[φ, b] uniquely solves the direct problem defined by the excitation φ and the parameter vector b. The inverse problem basically consists in finding values of b such that the simulated measurements D I u[φ, b] agree with the actual measurements d I , with consistency with additional a priori information possibly required. In practice, the inverse problem is often formulated as the minimisation of an cost function b

  e. when C = C ⋆ , the overspecified boundary data (ξ, φ) are incompatible and E(C) > 0. If C = C ⋆ , (ξ, φ) are compatible and E(C) = 0. Let u D and u N denote the solutions of the Dirichlet and Neumann boundary value problems defined in terms of a given elasticity tensor C and overspecified boundary measurements (27), i.e. by div (C : ε

Figure 2 .

 2 Figure 2. Identificatication of a distribution of Poisson ratio: "true" distribution used to create the synthetic data (left); identified distribution after 5 iterations (center) and after 32 iterations (right).

Figure 3 .

 3 Figure3. Identificatication of a distribution of Poisson ratio: "true" distribution used to create the synthetic data (left); identified distribution after 2 iterations for noisy data (center) and after 14 iterations for noise-free data (right).

20 Figure 4 .

 204 Figure 4. Indentification of a distribution of the bending stiffness in an isotropic plate: "true" distribution used to create the synthetic data (left) and identified distribution after 20 iterations using the second eigenmodule (right).

3 Figure 5 .

 35 Figure 5. Contour maps of J (C) (top), H(C) with (A, B) = (50, 50) (middle) and H(C) with (A, B) = (1, 50) (bottom), plotted as functions of the location of the center of a trial four-element square defect. The graphs in the left and right columns correspond to no noise and a simulated relative noise of amplitude 1% on the synthetic data, respectively. The white square indicates the location of the "true" four-element defect.

Figure 6 .

 6 Figure 6. Error in constitutive equation and free vibrations (synthetic example): FE mesh and distribution of Young modulus showing the simulated defect.

Figure 7 .

 7 Figure 7. Error in constitutive equation (ECE) and free vibrations (synthetic example): FE mesh and distributed ECE computed from one perfectly known modal displacement.

Figure 8 .

 8 Figure 8. Description of the test structure and FEM model.

Figure 9

 9 Figure 9. ECE distributions (normalized to max=1) for the "small" (elements 61-62, left) and "large" (elements 61-62-63-64, right) defects, using real data for 13 measured eigenmodes.

6. 1 . 1 .

 11 Sensitivity formula in domain integral form. First, the directional derivative of the cost function[START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] in a defect shape transformation defined by the transformation velocity θ(x) (x ∈ Γ) is given, by virtue of identity[START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF], by dJ db = Sp ∂ϕ u ∂u .p Γ ) denote the Lagrangian derivative of the solution to the direct problem. The partial derivative of a real-valued function with respect to a complex vector is conventionally defined as∂ϕ ∂w ≡ ∂ϕ ∂w R i ∂ϕ ∂w I w R = Re(w) , w I = Im(w)(125)Upon application of the differentiation formulae[START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF] to the weak formulation[START_REF] Nintcheu Fata | A Computational Basis for Elastodynamic Cavity Identification in a Semi-Infinite Solid[END_REF], the Lagrangian derivative ( [σ[u Γ ] : ∇wρω 2 u Γ .w]div θ dV + Ω Γ σ[u Γ ].∇w + σ[w].∇u Γ .∇θ dV + Ω Γ ∇u Γ : (∇C.θ) : ∇w dV = 0 (127)

T

  (x o ) = Re ∇w.D(B, β, γ).∇u + (βγ 2 -1)k 2 wu (x o )

Example 1 : 1 . 1 .

 111 three hard spherical objects. A set of three hard spherical objects (with respective centers (-3, 2, -3), (2, 3, -3) and (0, -4, -4) and radii 1, 1 and 2) are embedded in the half-space x 3 ≤ 0. The surrounding medium is characterized by the wavenumber k = The free-field u is created by the successive application of four unit point sources x S = (±5, ±5, 0) on the external surface S = {x 3 | x 3 = 0} (a homogeneous Neumann

3 Figure 10 .

 310 Figure 10. Identification of a set of three hard spherical objects (with respective centers (-3, 2, -3), (2, 3, -3) and (0, -4, -4) and radii 1, 1 and 2) in an acoustic half-space: distribution of T in the x 3 = -3 (horizontal) plane. The horizontal contours of the three objects are shown (white circles).

- 5 Figure 11 .

 511 Figure 11. Identification of a penetrable obstacle of ellipsoidal shape (center (-1, 2, -3), semiaxes 1, 0.6, 0.4 aligned with the global Cartesian frame, wave velocity contrast γ = 0.7, mass density contrast β = 0.45) embedded in an acoustic half-space: distribution of T in nine regularly spaced horizontal planes x 3 = a with a = -1, . . . , -5. The horizontal contour of the "true" penetrable obstacle is shown (white ellipse). The same color scale is used in all nine graphics. Case k = 1.

- 5 Figure 12 .

 512 Figure 12. Identification of a penetrable obstacle of ellipsoidal shape (center (-1, 2, -3), semiaxes 1, 0.6, 0.4 aligned with the global Cartesian frame, wave velocity contrast γ = 0.7, mass density contrast β = 0.45) embedded in an acoustic half-space: distribution of T in nine regularly spaced horizontal planes x 3 = a with a = -1, . . . , -5. The horizontal contour of the "true" penetrable obstacle is shown (white ellipse). The same color scale is used in all nine graphics. Case k = 4.

Table

		E (GPa)	ν	ρ (kg/m 3 ) σ Y (GPa)	ε Y
	aluminium	71	0.34	2.6	150 -400 .002 -.006
	steel	210	0.29	7.8	200 -1600 .002 -.007
	titanium	105	0.34	4.5	700 -900 .006 -.009
	marble	26	0.3	2.8	10	.0004
	glass	60	0.2 -0.3 2.5 -2.9	1200	.02

⋆ u Γ dS + Su ∂ϕ p ∂p . ⋆ p Γ dS + Γ (∇ψ.θ + ψdiv S θ) dS(124)where the derivative d/d is implicitly evaluated at b = 0, the domain transformation is assumed to leave S unchanged, and ( ⋆ u Γ , ⋆

⋆u, ⋆ p) of the solution to the direct problem[START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF] solves the weak formulationΩ Γ σ[ ⋆ u Γ ] : ∇wρω 2 ⋆ u Γ .w dV -Su ⋆ p Γ .w dS + Ω Γ ∇u Γ : (∇C.θ) : ∇w dV + Ω Γ σ[u Γ ] : ∇wρω 2 u Γ .w div θ + σ[u Γ ].∇w + σ[w].∇u Γ .∇θ dV = 0 (126)where (w, q) again denote unconstrained trial fields and the terms containing ⋆ w, ⋆ q have been ignored (they merely reproduce the direct problem in weak form and thus collectively vanish).The derivatives( ⋆ u Γ , ⋆ p Γ) can be computed by solving the derivative weak formulation[START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF] and the result substituted into[START_REF] Ramananjaona | Shape reconstruction of buried obstacles by controlled evolution of a level set: from a min-max formulation to numerical experimentation[END_REF]. This is known as the direct approach to sensitivity. When applied to the computation of the gradient of J (Γ) it entails the solution of one derivative problem[START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF] for each geometrical parameter involved.A more efficient approach when (as here) the direct problem is linear consists in eliminating the derivatives ( ⋆ u, ⋆ p) by means of an adjoint solution, i.e. a suitable choice of the trial fields (w, q). To this purpose, and following the general approach of sections 4.2 and 4.3, consider the restriction of the weak formulation[START_REF] Reynier | Sur le contrôle de modélisations éléments finis: recalage à partir d'essais dynamiques[END_REF] to trial fields (w, q) such that

(A C + ω 2 )w = 0 and q = σ[w]. On using the constitutive symmetry σ[ ⋆ u] : ∇w = σ[w] : ∇ ⋆ u, applying the divergence formula to its first integral, equation (126) yields the identitySp ⋆ u Γ .q dS + Γ ⋆ u Γ .q dS -Su ⋆ p Γ .w dS + Ω Γ

K[w] I,II,III (s) are known as the stress intensity factors (SIFs). Expansions similar to[START_REF] Sokolowski | Introduction to shape optimization[END_REF] hold for the adjoint solution w, with SIFs K[w] I,II,III . The universal angular functions f I (θ), f II (θ), g I (θ), g II (θ) can be found in textbooks on fracture mechanics such as[START_REF] Unger | Analytical fracture mechanics[END_REF]. Since by virtue of these expansions σ[u Γ ] : ε[w Γ ] = O(1/r), the integral over D ε in (133) vanishes in the limit ( dV = r(1 + O(r)) drdθ ds in D ε ). Besides, it can be verified that [[σ[u SΓ ] : ∇w S Γ ]] = 0, i.e. that the most singular contributions cancel out in the integral overΓ ε . Hence, [[σ[u Γ ] : ∇w Γ ]] = O(d 1/2), and the integral over Γ ε becomes in the limit ε → 0 the corresponding weakly singular (convergent) integral over Γ. Finally, under mild smoothness assumptions on the closed curve ∂Γ and the velocity field θ, one has:Σε [σ Γ : ∇w Γρ uΓ . ẇΓ ]θ n -[p Γ .∇w Γ + q Γ .∇u Γ ].θ dS = ∂Γ θ n (s) π -π σ[u S Γ ] : ∇w S Γ (ε, θ, s)εdθ ds -∂Γ θ(s) • π -π p[u S Γ ].∇w S + p[w S Γ ].∇u S (ε, θ, s)εdθ ds + O(ε 1/2 )The integrals in the right-hand side yield a finite contribution in the limit ε → 0, which can be evaluated in a straightforward way using expansions[START_REF] Sokolowski | Introduction to shape optimization[END_REF]. This results in the following expression of ⋆ J(Γ):⋆ J (Γ) = Γ θ n σ[u Γ ] : ∇w Γρω 2 u Γ .w Γ dS -1 µ ∂Γ θ ν K[u Γ ] III K[w Γ ] III ds

Then, the counterparts of the sensitivity formulae[START_REF] Salenc ¸on | Handbook of continuum mechanics (general concepts, thermoelasticity)[END_REF],[START_REF] Schmerr | Ultrasonic flaw sizing inverse problems[END_REF] and[START_REF] Stavroulakis | Nondestructive elastostatic identification of unilateral cracks through BEM and neural networks[END_REF] are respectively⋆ J (Γ) = T 0 Ω Γ σ[w Γ ] : ∇u Γ + ρ üΓ .w Γ div θ σ[w Γ ].∇u Γ + σ[u Γ ].∇w Γ : ∇θ dV dt(138)⋆ J (Γ) = T 0 Γ σ[w Γ ] : ∇u Γ + ρ üΓ .w Γ θ n dS dt(139)⋆ J (Γ) = Γ

direct and the adjoint solutions, one obtains

where

denotes the jump of f across Γ. Now, the limiting form when ε → 0 of Eq. ( 133) is sought. Recall the well-known expansions of the mechanical fields near the fixed crack front, isotropic elasticity being assumed:

where (r, θ) denote local polar coordinates, attached to a point x(s) of ∂Γ characterized by its arc length s, in the plane orthogonal to ∂Γ and emanating from x(s). The functions

where ν is the unit normal to ∂Γ lying in the tangent plane of Γ and pointing outward of Γ, and θ ν = θ.ν. Again, the identity ( 132) is useful.

Shape sensitivity formulae in the time domain

Without going again through a detailed derivation, sensitivity formulae similar to ( 129), ( 131) and ( 135) can be established in connection with elastodynamics in the time domain. Assuming measurements have been obtained over a time interval 0 ≤ t ≤ T , the generic cost function J (Γ) = J(u Γ , p Γ , Γ) has in that case the format

instead of [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF]. The adjoint solution (w Γ , q Γ ) fulfills the boundary conditions [START_REF] Rus | Crack shape sensitivity by the adjoint variable method using a boundary-only integral[END_REF] and the final conditions w(x, T ) = ẇ(x, T ) = 0