Marc Bonnet
email: bonnet@lms.polytechnique.fr

Bojan B Guzina
email: guzina@wave.ce.umn.edu

Elastic-wave identification of penetrable obstacles using shape-material sensitivity framework ⋆

Keywords: Shape-material sensitivity, Elastodynamics, Identification, Inclusion, Boundary element method, Constrained optimization

Elastic-wave identification of penetrable obstacles using shape-material sensitivity framework

Marc Bonnet, B. B.

Introduction

Elastic-wave sensing of penetrable (i.e. deformable) heterogeneities in a solid matrix is a long-standing problem in mechanics with applications to nondestructive material testing, seismic prospecting, medical diagnosis, and un-derground object identificaton. In the context of seismic exploration, comprehensive three-dimensional (3D) mapping of subterranean structures commonly entails the interpretation of a large number, often thousands, of motion measurements using elastodynamic or acoustic models based on domain discretization, see e.g. [START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF]. In contrast, this investigation focuses on the mapping of objects buried in a known reference solid, from only a limited number of remote measurements. In such instances, boundary integral equation (BIE) formulations [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF] provide a direct mathematical link between the observed waveforms and the geometry and material characteristics of a hidden object, and therefore allow to exploit effectively the limited data.

Ahtough inverse scattering in general has been the subject of intensive mathematical and computational research [START_REF] Cakoni | Qualitative methods in inverse scattering theory[END_REF][START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF][START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF][START_REF]Scattering and Inverse Scattering in Pure and Applied Science[END_REF], only limited efforts have so far been devoted to the wave-based reconstruction of homogeneous elastic inclusions. Two-dimensional BIE formulations of the inclusion identification problem were proposed in [START_REF] Lee | Identification of geometric shapes and material properties of inclusions in two-dimensional finite bodies by boundary parameterization[END_REF][START_REF] Marin | Identification of material properties and cavities in two-dimensional linear elasticity[END_REF] (elastostatics) and [START_REF] Rus | Solution of identification inverse problems in elastodynamics using semi-analytical sensitivity computation[END_REF] (elastodynamics). Volume integral equations of the Lippmann-Schwinger type, which entail domain discretization over a region in which flaws are a priori expected, are used in e.g. [START_REF] Wiskin | Inverse scattering from arbitrary two-dimensional objects in stratified environments via a Green's operator[END_REF] with seismic waves idealized as acoustic waves and [START_REF] Abubakar | Iterative forward and inverse algorithms based on domain integral equations for three-dimensional electric and magnetic objects[END_REF][START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] in conjunction with the contrast source inversion concept for elastic and electromagnetic waves, while fast solution methods are proposed in [START_REF] Hohage | Fast numerical solution of the electromagnetic medium scattering problem and applications to the inverse problem[END_REF]. Mathematical results on identifiability are given in e.g. [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF]. More recently, approximate identification methods based on the small-inclusion asymptotics [START_REF] Ammari | Reconstruction of elastic inclusions of small volume via dynamic measurements[END_REF][START_REF] Ammari | A boundary integral method for computing elastic moment tensors for ellipses and ellipsoids[END_REF][START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF][START_REF] Kang | Identification of elastic inclusions and elastic moment tensors by boundary measurements[END_REF], specialized analytical solutions [START_REF] Baganas | Detection of spherical inclusions in a bounded, elastic cylindrical domain[END_REF], or energy considerations [START_REF] Alessandrini | Evaluating the volume of a hidden inclusion in an elastic body[END_REF] were proposed for preliminary "scanning" of solid bodies.

The focus of this investigation is the development of a computational platform for the 3D identification of penetrable elastic inclusions via an elastodynamic BIE framework. This approach rests upon the full three-dimensional elastodynamic model, with no resort to approximations such as the Born linearization or small-inclusion asymptotics. The inclusions are assumed homogeneous and bonded to the surrounding reference medium; as such, they are characterized by their boundary (i.e. the surface that separates them from the reference medium), elastic tensor, and mass density. For identificatiom purposes, the inverse problem is reduced to the minimization of a cost functional representing the misfit between experimental observations (values of displacements at sensor locations) and their simulations for an assumed inclusion configuration. The latter are based on a coupled system of regularized boundary integral equations [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF]. For computational efficiency of the gradient search technique employed by the inverse solution, the shape sensitivity of the featured cost functional is evaluated via an adjoint field approach which, besides the matter of elegance, is computationally much more efficient than finitedifference evaluations. This is accomplished by generalizing upon the shape sensitivity approach proposed in [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF] for elastic-wave void identification and in [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] for the inverse scattering of acoustic waves. To completely characterize penetrable elastic defects, the material sensitivity of the cost functional is derived using two alternative methodologies: (i) a direct differentiation approach based on the material-parameter derivative of the governing BIE pair, or (ii) adjoint field approach. A similar BIE-based treatment of shape-material sensitivity has been recently proposed [START_REF] Zacharopoulos | Three dimensional reconstruction of shape and piecewise constant region values for Optical Tomography using spherical harmonic parameterisation and a Boundary Element Method[END_REF] for optical tomography featuring the scalar Helmholtz equation with a complex wavenumber, while e.g. [START_REF] Akçelik | Shape determination for deformed electromagnetic cavities[END_REF][START_REF] Park | Shape identification for natural convection problems using the adjoint variable method[END_REF] present other applications of adjoint-based shape sensitivity analyses.

Making use of the message-passing interface, shape and material sensitivities derived in this study are implemented in a data-parallel code and integrated into a nonlinear optimization algorithm towards the solution of the 3D inverse problem. Preliminary identification studies on sample inclusion configurations demonstrated that the unconstrained (quasi-Newton) minimization algorithm, that was successfully used in previous studies [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] for 3D void reconstruction, performs unsatisfactorily for the geometric-material identification problem at hand. For this reason, a more robust algorithm based on the augmented-Lagrangian cost functional [START_REF] Nocedal | Numerical Optimization[END_REF] has been developed as a means to deal with physical inequality constraints, without increasing the dimension of the parametric space, and ensure that all iterates (in terms of the inclusions' geometric and material parameters) be physically admissible. The elastodynamic transmission problem, whose repeated solution is required by the minimization, is thus always well-posed. The numerical results, which employ a direct boundary element method (BEM) for the primary, adjoint and material-sensitivity solutions, demonstrate the feasibility of identifying the geometry and material characteristics of penetrable (elastic) defects hidden in a semi-infinite solid from only a limited number of waveform measurements taken on the (traction-free) surface.

Direct and inverse scattering by elastic inclusions

Consider an inverse scattering problem where the reference homogeneous solid Ω, containing a bonded inclusion Ωtrue with boundary Γ true , is probed by elastic waves. The reference medium, whose external boundary (available for testing) is denoted by S, is characterized by its elastic tensor C and mass density ρ; the respective material characteristics of the inclusion are denoted as Ĉtrue and ρtrue . The ensuing shape and material sensitivity analyses (sections 3 to 5) are carried out for the reference solid of an arbitrary shape, whereas the computational treatment and numerical results presented thereafter (sections 6, 7) assume a semi-infinite configuration whereby S denotes the traction-free surface of the half-space, and isotropic elastic properties in Ω -and Ωtrue .

Inverse problem. To identify the geometry Ωtrue and material characteristics Ĉtrue and ρtrue of the hidden defect, time-harmonic excitations are applied in the form of volume (f ) and surface (g) force densities having respective supports V ⊂ Ω and S t ⊂ S, and displacements u D over the complementary surface S u = S\S t . The implicit time-harmonic factor exp(iωt) where ω denotes the angular frequency of excitation is omitted hereon. Letting Ω denote a trial inclusion bounded by Γ and Ω -= Ω \ (Γ ∪ Ω) be the region surrounding the obstacle, the prescribed excitation (f , g, u D ) gives rise to elastodynamic displacement fields u = u[ Ω, Ĉ, ρ] in Ω -and û = û[ Ω, Ĉ, ρ] in Ω.

For identification purposes, the displacement u obs induced in the flawed solid by (f , g, u D ) is monitored over the measurement surface S obs ⊂ S t (other possibilities, e.g. finite sets of measurement points, being also allowed by the ensuing treatment). Ideally, a defect configuration ( Ωtrue , Ĉtrue , ρtrue ) such that

u[ Ωtrue , Ĉtrue , ρtrue ] = u obs (on S obs ) (1)
is sought, where u on S obs is understood in the sense of the trace. In practice, due to many factors (e.g. incomplete and/or inexact measurements, modelling uncertainties), the inclusion is sought so as to minimize a misfit cost functional

J ( Ω, Ĉ, ρ) = S obs ϕ u[ Ω, Ĉ, ρ](ξ), u obs (ξ), ξ dS ξ , (2)
where function ϕ, which quantifies the misfit between the predicted and observed displacements, is assumed to be differentiable with respect to its arguments. For example, the usual least-squares measure of misfit is defined through 2ϕ(u,

u obs , ξ) = | u -u obs | 2 .
Forward problem. Let L and L denote the Navier partial differential operator respectively associated with the reference solid and inclusion, i.e.

Lu ≡ div (C : ∇u) + ρω 2 u, Lû ≡ div ( Ĉ : ∇û) + ρω 2 û. (3)
In ( 3) and thereafter, all quantities defined with reference to the inclusion or its constitutive parameters are indicated with a hat symbol. Moreover, the nabla symbol ∇ denotes the gradient operator, with the convention ∇(•) = (•) ,ℓ ⊗e ℓ (the comma denoting a partial derivative), while the column symbol ':' indicates a two-fold inner product between tensors, e.g. (C : ∇u

) ij = C ijkℓ u k,ℓ .
The predicted displacement featured in cost functional (2) solves the forward problem

(a) Lu + f = 0 (in Ω -), (b) Lû = 0 (in Ω), (4)
u = û, t + t = 0 (on Γ), (5)
(a) t = g (on S t ), (b) u = u D (on S u ), (6)
where equations ( 4), ( 5) and ( 6) respectively state the elastodynamic field equations, the perfect-bonding interfacial transmission conditions, and the external boundary conditions. In ( 5) and ( 6), t ≡ (C : ∇u)•n and t ≡ ( Ĉ : ∇û)•n are the boundary tractions relative to the reference medium and the inclusion, respectively, with n = -n and n denoting the unit normal exterior to Ω -.

If the reference domain Ω extends to infinity in any direction, as is the case for the semi-infinite solid examined later, u must in addition satisfy a suitable radiation condition at infinity (this can be relaxed so as to allow e.g. scattering of incident plane waves, see Appendix A.2). Moreover, (4) implicitly carries the assumption Ω ∩ V = ∅, i.e. the forward problem ( 4)-( 6) is considered only for a inclusion that is separated from the body force support.

Two alternative formulations of problem ( 4)-( 6) are now summarized: (i) the weak formulation, upon which the general results in terms of sensitivity formulas are established (Secs. 3 to 5), and (ii) the BIE formulation used here as a basis for the computational treatment and numerical results (Secs. 6, 7).

Weak formulation. The forward transmission problem ( 4)-( 6) can be recast in weak form whereby (u, û) ∈ V(u D ) must satisfy

A (u, û), (w, ŵ) -F (w) = 0 ∀(w, ŵ) ∈ V(0), (7)
with the function spaces V, the symmetric bilinear form A, and the linear form F defined by

V(u D ) = (w, ŵ) w ∈ H 1 loc (Ω -) 3 , ŵ ∈ H 1 ( Ω) 3 , w = u D on S u , w = ŵ on Γ , (8a) A (u, û), (w, ŵ) = Ω - a(u, w) dV + Ω â( û, ŵ) dV (8b) F (w) = V f •w dV + St g•w dS (8c)
and the bilinear energy densities a in Ω -and â in Ω given by

a(u, w) = ∇u : C : ∇w -ρω 2 u•w, (9a) â
( û, ŵ) = ∇û : Ĉ : ∇ ŵ -ρω 2 û• ŵ. ( 9b
)
BIE formulation (semi-infinite solid). Wave propagation and scattering in an elastic half-space is a suitable idealization for a number of applications such as non-destructive material testing and seismic exploration. BIE formulations, which deal effectively with unbounded domains, are a natural framework for such configurations. With reference to the Cartesian frame {O, ξ 1 , ξ 2 , ξ 3 }, let the host domain Ω be semi-infinite (ξ 3 ≥ 0) and bounded by the traction-free surface S = {ξ | ξ 3 = 0}. Let U (x, ξ) and T (x, ξ) denote the (half-space) elastodynamic Green's tensors, defined such that U iℓ and T iℓ (i, ℓ = 1, 2, 3) respectively denote the ith component of the displacement and traction at ξ ∈ Ω resulting from a unit time-harmonic point force applied at x ∈ Ω in the ℓth direction, with T iℓ vanishing identically on S. Similarly, let Û (x, ξ) and

Γ u(ξ) -u(x) •[ T (x, ξ)] 1 dS ξ + Γ u(ξ)•[ T (x, ξ)] 2 dS ξ + Γ t(ξ)• Û (x, ξ) dS ξ = 0, x ∈ Γ (11)
written respectively for the "matrix" Ω -and inclusion Ω, in terms of the traces (u, t) on Γ of the exterior field. The free field u F featured in [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] is the solution to (4), [START_REF] Ammari | Reconstruction of elastic inclusions of small volume via dynamic measurements[END_REF] when no inclusion is present, and is explicitly given by

u F (x) = V f (ξ)•U (x, ξ) dV ξ + S g(ξ)•U (x, ξ) dS ξ , x ∈ Ω. ( 12
)
For regularization purposes, the traction Green's functions in [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] and [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF] are decomposed [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF] into the sum of (frequency-independent) singular parts

[T ] 1 , [ T ] 1 and (frequency-dependent) regular parts [T ] 2 , [ T ] 2 according to T (x, ξ) = [T (x, ξ)] 1 + [T (x, ξ)] 2 , T (x, ξ) = [ T (x, ξ)] 1 + [ T (x, ξ)] 2 .
On solving [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] and [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF] for u and t, the displacement in the (reference) solid surrounding the inclusion is given by the integral representation formula

u(x) = u F (x) + Γ t(ξ)•U (x, ξ) -u(ξ)•T (x, ξ) dS ξ , x ∈ Ω -. ( 13
)
Letting S R denote the sphere of radius R centered at the origin, [START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF] and [START_REF] Cakoni | Qualitative methods in inverse scattering theory[END_REF] rest on the assumption that u is a radiating elastodynamic state in the semiinfinite solid Ω -, whereby u, t satisfy the generalized radiation condition [START_REF] Madyarov | A Radiation condition for layered elastic media[END_REF] lim

R→+∞ S R ∩Ω u(ξ) • T (x, ξ) -t(ξ) • U (x, ξ) dS ξ = 0, x ∈ Ω -. ( 14
)

Differentiation with respect to inclusion perturbations

To quantify the effect of the inclusion's boundary and material-parameter perturbations on the cost function (2), the defect configuration ( Ω; Ĉ, ρ) is assumed to depend on a time-like evolution parameter τ [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], with the unperturbed, 'initial' configuration ( Ω; Ĉ, ρ) conventionally associated with τ = 0. In this study, only the first-order infinitesimal perturbations of ( Ω; Ĉ, ρ), i.e. the first-order "time" derivatives at τ = 0, are considered. Perturbations of the inclusion's constitutive properties can thus be expressed as

Ĉτ = Ĉ + Ĉ′ τ, ρτ = ρ + ρ′ τ, (15)
whereas the shape perturbations of Ω can be synthesized as

x ∈ Ω → x τ = x + θ(x)τ ∈ Ω τ , (16)
where θ(x) is a given (i.e. prescribed) transformation velocity field. In the sequel, θ(x) is assumed to vanish outside of a neighbourhood of Ω, which postulates the existence of a bounded region O satisfying

(a) Ω ⊂ O ⊂ Ω, (b) θ = 0 in Ω \ O, (c) O ∩ S = ∅, (d) O ∩ V = ∅, (17)
with (c) and (d) stemming from the natural assumption that the trial inclusion intersect neither the external surface nor the body-force support.

Total and Lagrangian derivatives of field variables. Consider a generic field variable g which depends on the inclusion configuration ( Ωτ , Ĉτ , ρτ ), e.g. the solution of the forward problem (4)- [START_REF] Ammari | A boundary integral method for computing elastic moment tensors for ellipses and ellipsoids[END_REF]. Such quantity can be represented in the form g(x ; Ωτ , Ĉτ , ρτ ). Let the total derivative, ⋄ g, of g be defined by

⋄ g = lim τ ց0 1 τ g(x τ ; Ωτ , Ĉτ , ρτ ) -g(x ; Ω, Ĉ, ρ) , (18)
i.e. by following the evolution of g at a point x τ moving according to geometric transformation [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF], while the inclusion's material parameters are perturbed according to [START_REF] Colton | Integral Equation Method in Scattering Theory[END_REF]. The contributions to ⋄ g of the inclusion's geometric and material perturbations can be separated through additive decomposition

⋄ g = • g + g ′ , (19)
with the shape sensitivity

• g and material sensitivity g ′ respectively defined by "freezing" the material parameters and the shape of the inclusion in [START_REF] Guzina | Topological derivative for the inverse scattering of elastic waves[END_REF], i.e.

• g = lim τ ց0 1 τ g(x τ ; Ωτ , Ĉ, ρ) -g(x ; Ω, Ĉ, ρ) (20)
g ′ = lim τ ց0 1 τ g(x ; Ω, Ĉτ , ρτ ) -g(x ; Ω, Ĉ, ρ) = ∂g ∂ Ĉ • Ĉ′ + ∂g ∂ ρ ρ′ . (21)
The shape sensitivity thus corresponds to the Lagrangian time derivative of continuum kinematics with the physical time variable replaced with a pseudotime. For the ensuing developments, it is useful to note that the shape sensitivity of a gradient is given [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] by

(∇g) • = ∇ • g -∇g • ∇θ. (22)
Differentiation of integrals. Consider a generic domain integral I of the form

I = O g(x ; Ωτ , Ĉτ , ρτ ) dV,
in Eulerian description, where O denotes the support of geometric perturbation [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]. Upon noting that the shape sensitivity of dV is given by

• dV = [div θ] dV , the total derivative of I [39] can be written as (a) ⋄ I = • I + I ′ , (b)
• I = O • g + gdiv θ dV, (c) I ′ = O g ′ dV. ( 23
)
4 Shape sensitivity using an adjoint solution

Taking into account assumption (17c) whereby S obs remains invariant under perturbation [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF], the shape sensitivity of cost functional (2) takes the form

• J = Re S obs ϕ ,u (u, u obs , ξ)• • u dS , (24)
where the complex-valued function ϕ ,u is defined as

ϕ ,u = ∂ϕ ∂u R -i ∂ϕ ∂u I , (u R = Re(u), u I = Im(u)). ( 25
)
On parameterizing the geometric perturbation ( 16) in terms of a selected shape parameter a and setting τ = aa 0 , formula [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF] yields the derivative of J with respect to a evaluated at a 0 in terms of the derivative solution

• u; this is the essence of the so-called direct differentiation approach. One can however circumvent the actual computation of derivatives • u (one for each shape parameter describing the sought inclusion) by resorting to the adjoint field approach, to whose formulation the remainder of this section is devoted.

Displacement shape sensitivity -weak formulation

Since the weak formulation [START_REF] Baganas | Detection of spherical inclusions in a bounded, elastic cylindrical domain[END_REF] of the forward problem ( 4)-( 6) holds for all perturbed inclusion configurations, the governing weak formulation for the displacement shape sensitivity ( • u,

• û) can be obtained by exploiting the identity

• A (u, û), (w, ŵ) - • F(w) = 0 ∀(w, ŵ) ∈ V(0). (26)
As shown in Appendix A.1, application of the Lagrangian differentiation formula (23b) to (8b) results in

• A (u, û), (w, ŵ) = A ( • u, • û), (w, ŵ) + A (u, û), ( • w, • ŵ)
+

Ω - Lw•(∇u•θ) dV + Ω L ŵ•(∇û•θ) dV + Γ n•E(u, w) + n• Ê( û, ŵ) •θ dS, ( 27
)
where the bilinear tensorial function E(u, w), related to the dynamic Eshelby energy-momentum tensor [START_REF] Gurtin | Configurational forces as basic concept of continuum physics[END_REF], is defined by

E(u, w) = a(u, w)I -(C : ∇w)•∇u -(C : ∇u)•∇w. (28)
(I: second-order identity tensor) and Ê( û, ŵ) is defined similarly in terms of the inclusion characteristics. Invoking assumption (17d), one further finds that the Lagrangian derivative of F (w) is given by

• F(w) = F ( • w). (29)
On substituting ( 27) and ( 29) into [START_REF] Kirsch | An introduction to the mathematical theory of inverse problems[END_REF], invoking the equality obtained by setting (w, ŵ) = (

• w,

• ŵ) in [START_REF] Baganas | Detection of spherical inclusions in a bounded, elastic cylindrical domain[END_REF], and noting that • u = 0 on S u (since the prescribed displacement u D is insensitive to the inclusion shape), the displacement shape sensitivity ( 0) is found to be governed by the weak formulation:

• u, • û) ∈ V(
A ( • u, • û), (w, ŵ) = - Ω - Lw•(∇u•θ) dV - Ω L ŵ•(∇û•θ) dV Γ n•E(u, w) + n• Ê( û, ŵ) •θ dS ∀(w, ŵ) ∈ V(0). ( 30
)

Adjoint solution

The main motivation behind the adjoint state approach is to evaluate the shape sensitivity [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF] in an indirect, and computationally faster, manner by circumventing the actual computation of field sensitivities • u. Interpreting the integral in [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF] as a virtual work and treating the sensitivity • u therein as a test function leads to define the adjoint state (v, v) ∈ V(0) by the weak statement

A (v, v), (w, ŵ) = S obs ϕ ,u (u, u obs , ξ) • w dS ∀(w, ŵ) ∈ V(0), (31)
whose solution (v, v) is the adjoint state. The adjoint transmission problem [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF] can equivalently be stated in strong form as:

(a) Lv = 0 (in Ω -), (b) Lv = 0 (in Ω)            (c) p = ϕ ,u (on S obs ) (d) p = 0 (on S t \ S obs ) (e) v = 0 (on S u ) (f) v = v, p + p = 0 (on Γ) ( 32
)
where p = (C : ∇v)•n and p = ( Ĉ : ∇v)• n are the traction vectors respectively associated with v and v. Alternatively, traces (v, p) on Γ of the solution to problem (32) satisfy integral equations ( 10), [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF] with the free field given by

v F (x) = S obs ϕ ,u u(ξ), u obs (ξ), ξ •U (x, ξ) dS ξ , x ∈ Ω. (33)
The nature (32c) of the adjoint excitation and the assumed boundedness of S obs ensure that (v, p) satisfies the generalized radiation condition (14).

Shape sensitivity formula

Setting (w, ŵ) = (

• u,

• û) in [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF], formula [START_REF] Ikehata | Reconstruction of the shape of the inclusion by boundary measurements[END_REF] for the shape sensitivity

• J becomes • J = Re A (v, v), ( • u, • û) .
Choosing (w, ŵ) = (v, v) in [START_REF] More | Line search algorithms with guaranteed sufficient decrease[END_REF] then readily yields, by virtue of field equations (32a,b) and the symmetry of A(•, •), the following expression for • J , where the displacement shape sensitivity no longer appears:

• J = -Re Γ n•E(u, v) + n• Ê( û, v) •θ dS . ( 34
)
Equation ( 34) is, however, not well suited for applications where the free and adjoint solutions are computed by means of the BEM, as it features displacement gradients on Γ. This is addressed by introducing the decomposition

∇u = ∇ S u + u ,n ⊗ n ( 35
)
of a gradient in terms of its tangential component ∇ S u and the normal derivative u ,n , and expressing the latter in terms of ∇ S u and t by inverting the relationship t = (C : ∇u)•n whereby

u ,n = D•∆t. (36)
Here, the second-order tensor D and the combination ∆t are respectively defined by D = [n•C•n] -1 (exploiting the minor symmetry of C), and

∆t = t -(C : ∇ S u)•n. ( 37
)
By virtue of ( 35)-( 37), tensor function E(u, v) can be written in terms of the interfacial tractions and tangential displacement gradients, so that

n•E(u, v)•θ = ∇ S u : C : ∇ S v -∆t•D•∆p -ρω 2 u•v θ n -(t•∇ S v + p•∇ S u)•θ.
This finally allows to express [START_REF] Nocedal | Numerical Optimization[END_REF] in terms of quantities directly available from the boundary element solution. Utilizing transmission conditions ( 5) and (32f), the desired form of the shape sensitivity result is thus established as

• J = Re Γ ∇ S u : ( Ĉ -C) : ∇ S v -∆ t• D•∆p + ∆t•D•∆p -(ρ -ρ)ω 2 u•v θ n dS , (38)
with ∆ t, ∆p defined by [START_REF] Park | Shape identification for natural convection problems using the adjoint variable method[END_REF] with C, n replaced by Ĉ, n, and

D = [ n• Ĉ• n] -1 .
In the case where either material has isotropic elasticity (which is not a prerequisite for the derivation of ( 38)), tensors C, D or Ĉ, D are given in terms of the shear modulus µ and Poisson's ratio ν of the relevant material (with I denoting the symmetric fourth-order identity tensor) by

C = 2µ ν 1 -2ν I ⊗ I + I , D = 1 µ I - 1 2(1 -ν) n ⊗ n .

Material sensitivity

The material sensitivity of the generic cost function ( 2) is a priori given by

J ′ = Re S obs ϕ ,u (u, u obs , ξ)•u ′ dS . ( 39
)

Domain integral formulation, adjoint solution approach

Proceeding along the lines of Section 4, the domain formulation of the displacement material sensitivity (u ′ , û′ ) is found by considering the material sensitivity of weak statement [START_REF] Baganas | Detection of spherical inclusions in a bounded, elastic cylindrical domain[END_REF], i.e.

A (u ′ , û′ ), (w, ŵ) = - Ω â′ ( û, ŵ) dV, ∀(w, ŵ) ∈ V(0) (40)
where the bilinear density â′ is defined by (9b) with the inclusion characteristics Ĉ and ρ replaced respectively by their sensitivities Ĉ′ and ρ′ .

Expression [START_REF] Petryk | Time derivatives of integrals and functionals defined on varying volume and surface domains[END_REF] suggests that J ′ can be expressed in terms of the previouslydefined adjoint problem [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF]. Indeed, on setting (w, ŵ) = (v, v) in ( 40), (w, ŵ) = (u ′ , û′ ) in [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF], and proceeding as before, one finds that

J ′ = -Re Ω â′ ( û, v) dV . (41)
The domain-integral format of formula ( 41) is obviously impractical for BEMbased applications. For that reason, alternative approaches that facilitate the evaluation of J ′ using BEM forward and adjoint solutions are examined next.

Before proceeding to the BIE specialization of ( 41) and the subsequent numerical results, both focused here on semi-infinite host domains, it should be emphasized that the shape sensitivity formula [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] and its material counterpart [START_REF] Plessix | Waveform inversion of reflection seismic data for kinematic parameters by local optimization[END_REF] are general in the sense that they uniformly apply to finite, semiinfinite or infinite (R 3 ) host domains. The present BIE formulation and its underlying assumption of a semi-infinite "host" define an illustration, not a limitation, of the proposed shape-material sensitivity framework.

Surface integral formulation, direct approach

The direct differentiation approach makes use of ( 39) with the material sensitivity u ′ on S obs evaluated by differentiating the governing boundary integral equations. On expressing the governing pair ( 10) and [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF] in operator form as

T u (x) -U t (x) = u F (x) T u (x) -Û t (x) = 0 (x ∈ Γ), (42)
and keeping in mind that the full-space Green's tensors Û , T depend on the inclusion's material parameters whereas U , T and the free field u F do not, the sensitivities u ′ , t ′ on Γ can be shown to solve the pair of integral equations

T u ′ (x) -U t ′ (x) = 0 T u ′ (x) -Û t ′ (x) = Û′ t (x) -T ′ u (x) (x ∈ Γ) (43)
where integral operators Û ′ , T ′ featured on the right-hand side are defined in terms of the respective Green's tensor derivatives Û ′ , T ′ , see Appendix A.3.

Once equations ( 43) are solved for u ′ , t ′ on Γ, the displacement sensitivity u ′ on S obs follows by taking the material sensitivity of representation formula [START_REF] Cakoni | Qualitative methods in inverse scattering theory[END_REF].

On substituting the resulting expression into (39), the material sensitivity J ′ is finally given, using an operator notation similar to that in [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF], by

J ′ = Re S obs ϕ ,u (u, u obs , ξ)• U obs t ′ (x) -T obs u ′ (x) dS . ( 44
)

Surface integral formulation, adjoint solution approach

As an alternative to the above direct differentiation strategy, an adjoint field approach for the evaluation of J ′ as a surface integral may be formulated as follows. Let the new adjoint state (v, v) be defined as the solution of a transposed system of integral equations, written in weak form as

Γ v(x)• T [w](x) -U[t w ](x) + v(x)• T [w](x) -Û[t w ](x) dS x = S obs ϕ ,u • U obs [t w ](x) -T obs [w](x) dS x , ∀(w, t w ) ( 45
)
where w and t w are trial (vector) functions on Γ. Next, multiplying the first and the second equation in [START_REF] Rus | Solution of identification inverse problems in elastodynamics using semi-analytical sensitivity computation[END_REF] respectively by v(x) and v(x), integrating the result over x ∈ Γ, and adding the identities one obtains

Γ v(x)• T [u ′ ](x) -U[t ′ ](x) + v(x)• T [u ′ ](x) -Û[t ′ ](x) dS x = Γ v(x)• Û′ [t](x) -T ′ [u](x) dS x , ∀(w, t w ). ( 46
)
Setting (w, t w ) = (u ′ , t ′ ) in [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF], subtracting the resulting identity from [START_REF] Wiskin | Inverse scattering from arbitrary two-dimensional objects in stratified environments via a Green's operator[END_REF] and recalling the material sensitivity formula [START_REF] Sokolowski | Introduction to shape optimization[END_REF], material sensitivity J ′ follows in a surface-integral, adjoint-based form as

J ′ = Γ v(x)• Û ′ [t](x) -T ′ [u](x) dS x . (47)
The above material sensitivity formula is in particular well suited for use within a Galerkin BEM framework, as it involves two nested surface integrals.

Computational treatment

In what follows, an inclusion identification method based on the proposed shape-material sensitivity approach and an augmented-Lagrangian cost functional is implemented using a BEM framework and the underlying assumption that the background domain Ω is semi-infinite.

Boundary integral approximation

To illustrate the utility of sensitivities [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] and [START_REF] Sokolowski | Introduction to shape optimization[END_REF] for elastic-wave identification of penetrable defects, let the trial obstacle Ω be isotropic and described in terms of a finite set a = (a 1 , a 2 , . . . , a D , μ, ν, ρ) of D+3 geometric and material parameters, where μ, ν and ρ denote respectively the shear modulus, Poisson's ratio, and mass density of the defect. Introducing an auxiliary notation

J a (a) ≡ J Ω(a 1 , . . . , a D ), Ĉ(μ, ν), ρ (48)
to reflect the defect parametrization and making reference to [START_REF] Guzina | A linear sampling approach to inverse elastic scattering in piecewise-homogeneous domains[END_REF] and [START_REF] Guzina | On the stress-wave imaging of cavities in a semi-infinite solid[END_REF], the sensitivities ∂J a /∂a k used for minimizing J a (a) are computable according to Geometric:

∂J a ∂a k = • J | τ =a k , k = 1, 2, . . . D (49a)
Material:

∂J a ∂a k = J ′ | τ =a k , k = D+1, D+2, D+3. (49b)
Boundary element discretization. For the evaluation of surface integrals over ∂ Ω = Γ that are involved in the computation of • J and J ′ , one may assume Γ = K e=1 E e , where {E e } K e=1 are closed and non-overlapping surface elements. Following the usual approach, each boundary element E e ⊂ Γ is parametrized by a mapping E → E e that introduces local coordinates, η = (η 1 , η 2 ) ∈ E, where E is a polygonal domain in R 2 . The approximating boundary surface, Γ h = K e=1 E h e , can then be generated by interpolating a set of parameter-dependent nodes ξ q (a) ∈ E e with pre-defined mesh connectivity. To this end, the Q-noded approximation E h e of a generic surface element E e ⊂ Γ is written as

ξ(η) = Q q=1 N q (η) ξ q (a), ξ ∈ E h e , ξ q ∈ E e , η ∈ E, (50)
where N q are the relevant shape functions. By virtue of (50) the normal transformation velocity θ n in ( 38) is, for a given parameter a k , approximated as

θ n (ξ) ≡ θ k n (ξ) = Q q=1 N q (η) ∂ξ q ∂a k •n, ξ ∈ E h e , ξ q ∈ E e , k = 1, 2, . . . D (51)
where n denotes the unit normal to the surface. Assuming next the isoparametric representation of elastodynamic quantities, the boundary displacement (u) and traction (t) fields over Γ are interpolated at a generic point ξ ∈ E h e in terms of the nodal displacements u q and tractions t q at ξ q ∈ E e as

u = Q q=1 N q (η) u q , t = Q q=1 N q (η) t q . ( 52
)
Evaluation of shape sensitivities. The computation of • J and J ′ entails solving transmission problems associated respectively with the primary field u, the adjoint field v, and the material-sensitivity field u ′ for each material parameter associated with the inclusion (i.e. a total of five problems under the assumption of isotropic inclusions). With reference to ( 12), ( 33), ( 42), ( 43), (50) and (52), the discrete algebraic systems for these fields can be written as

Primary H Ũ -G T = U F , Ĥ Ũ -Ĝ T = 0, (53a)
Adjoint H Ṽ -G P = V F , Ĥ Ṽ -Ĝ P = 0, (53b)
Material-sensitivity

H Ũ′ -G T′ = 0, Ĥ Ũ′ -Ĝ T′ = Ĝ′ T -Ĥ′ Ũ (53c)
where vectors Ũ and T contain the nodal approximations of the primary (displacement and traction) fields u and t; coefficient matrices H, G, Ĥ and Ĝ approximate the respective integral operators T , U, T and Û in [START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF]; Ṽ and P are associated with the respective adjoint fields v and p; U F and V F collect the nodal values of the free fields ( 12) and (33); vectors Ũ′ and T′ refer respectively to u ′ and t ′ , while Ĥ′ and Ĝ′ are the matrix discretizations of T ′ and Û′ in [START_REF] Rus | Solution of identification inverse problems in elastodynamics using semi-analytical sensitivity computation[END_REF]. The coefficient matrices are the same for all three discretized systems, which allows for a computationally-effective solution of (53b) and (53c) once the primary problem (53a) has been solved.

In the context of [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF], however, the shape sensitivity computation requires not only the primary and adjoint fields on Γ, but also their surface gradients. By virtue of (52), the required surface gradient of u is approximated as

∇ S u = Q q=1 ∇ S N q (η) ⊗ u q (54)
over each boundary element, with a similar expression applying in terms of ∇ S v. To evaluate the tangential derivative ∇ S N q in (54), let the local companion basis {r 1 , r 2 , n} at any point ξ ∈ E h e be defined from the (differentiable) boundary element parametrization ξ = ξ(η) of (50) by

r 1 = Q q=1 ∂N q ∂η 1 ξ q , r 2 = Q q=1 ∂N q ∂η 2 ξ q , n = r 1 × r 2 r 1 × r 2 , (55)
so that vectors r 1 and r 2 are tangent to E h e . The surface gradient of shape functions ∇ S N q then takes an explicit form

∇ S N q = n × ∂N q ∂η 2 r 1 - ∂N q ∂η 1 r 2 q = 1, 2, . . . Q. ( 56
)
Since the triplet {r 1 , r 2 , n} defined by (55) forms a positive-oriented basis in R 3 , reversing the connectivity of a given element (i.e. listing its nodes in opposite order) leads to the sign-reversal of n through swapping of r 1 and r 2 . Hence, a desired orientation for the approximate surface Γ h can be achieved by adequately setting the mesh connectivity. Here, quantities pertaining to Ω - (background) and Ω (inclusion) in sensitivity formula [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF], and consequently matrices H and Ĥ in (53a-c), are defined in terms of the inward and outward orientations of Γ, respectively. Consistent element orientation is thus ensured by using two opposite mesh connectivity tables for Γ h (one "direct" for Ω -, the other "reverse" for Ω). This method allows for systematic generalization towards multiple-inclusion or nested-inclusion configurations. The surface gradients (54) are insensitive to the choice of mesh connectivity orientation.

Parallel computation

Owing to the high computational cost commonly associated with 3D inverse scattering, regularized boundary integral treatment [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF] of the primary, ad- 44) for

• J and J ′ , in a data parallel code using the message-passing interface (MPI) [START_REF] Pacheco | Parallel Programming with MPI[END_REF]. Data-type parallelism normally applies when identical operations are performed concurrently on multiple data items. With reference to [START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF] and (53a), such is the case with repeated, timeconsuming, computation of the elastodynamic Green's tensors U and T (underlying the evaluation of matrices H and G), in situations when the reference domain Ω is semi-infinite. In contrast to the fundamental solution for an infinite domain R 3 which is available in closed form (e.g. Û and T underpinning matrices Ĥ and Ĝ, see (A.8a,b)), elastodynamic Green's tensors for a semiinfinite solid are given as improper integrals [START_REF] Guzina | On the analysis of wave motions in a multi-layered solid[END_REF] whose numerical quadrature entails two to three orders-of-magnitude higher computational effort.

On denoting by n p the number of processes and by N the number of BEM degrees of freedom on Γ, the code is accordingly parallelized by block-distributing the computation of H and G matrices (both of dimension N ×N) where every participating process is assigned approximately N/n p columns of each array. As vectors U F and V F in (53a,b) also involve the (displacement) Green's tensor U for the reference domain, see [START_REF] Bonnet | Sounding of finite solid bodies by way of topological derivative[END_REF] and [START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF], their computation is similarly distributed among the participating processes. For consistency, matrices Ĥ and Ĝ are likewise computed in parallel fashion, even though this does not result in a meaningful reduction of the run time owing to the closed-form nature of Û and T . Once computed for the solution of the primary problem, the LU-factorized "global" coefficient matrix (combining H, G, Ĥ and Ĝ) is stored and reused for solving the adjoint and material-sensitivity problems.

To illustrate the performance of the parallel code the speed-up ratio Q s , i.e. the ratio of the elapsed time of a serial program over that of its parallel counterpart with n p processes, is considered. The computation is performed on the IBM BladeCenter H cluster equipped with 307 LS21 nodes, each containing two dual-core 2.6 GHz Opteron processors sharing 8GB of memory. For the purpose of comparison, the inclusion is described as a nine-parameter ellipsoid (D +3 = 9) characterized by its three centroidal coordinates, three semi-axes, and three material constants. The boundary element mesh approximating the surface of the defect has 650 nodes; the testing configuration is comprised of 25 uniaxial sources and 36 triaxial receivers located on the surface of a semi-infinite solid (see Fig. 1). Table 1 shows the CPU times per evaluation of the cost function J a and its sensitivities ∂J a /∂a k (k = 1, • • • , 9), the speedup ratios Q s , and a measure of the per-processor efficiency Q s /n p for the sample problem in a semi-infinite solid. As a point of reference, computation of the analogous problem when the host domain is infinite (Ω = R 3 ), takes 1 minute and 31 seconds on a single processor. In the ensuing examples, the halfspace calculations are performed with n p = 48 which represents a reasonable compromise between the speed-up ratio and per-processor efficiency.

Defect parametrization

The geometry of the trial defect Ω is, for the ensuing numerical experiments, described in terms of an ellipsoid whose principal axes are aligned with the reference Cartesian frame {O; ξ 1 , ξ 2 , ξ 3 }; its evolution within the host domain Ω is restricted to i) translation and ii) stretch along the principal axes. For problems involving identification of a single isotropic defect, such description entails the use of a nine-dimensional parametric space

a = c 1 d , c 2 d , c 3 d , α 1 d , α 2 d , α 3 d , μ µ , ν, ρ ρ (57)
(i.e. D = 6), which incorporates the defect's centroidal motion (c i , i = 1, 2, 3), principal stretches (α i , i = 1, 2, 3), and material characteristics (μ, ν, ρ). Parameters a k in (57) are moreover defined in dimensionless fashion using material characteristics (µ, ρ) of the reference solid and an arbitrary length scale d. With such definitions, analytical dependence of the nodal coordinates, ξ q = ξ q (a), of the surface mesh on the evolving defect boundary Γ is introduced as an affine deformation of the boundary element mesh for a reference unit sphere S (described by Lagrange coordinates (X 1 , X 2 , X 3 )) so that

ξ q i = c i + α i X q i , X q ∈ S, i = 1, 2, 3 (58)
assuming no summation over index i. On the basis of ( 57) and (58), one finds that the normal transformation velocities θ k n defined by (51) are given by

θ k n , k = 1, . . . , D = n 1 , n 2 , n 3 , ξ 1 -c 1 α 1 n 1 , ξ 2 -c 2 α 2 n 2 , ξ 3 -c 3 α 3 n 3 d. ( 59
)
Since formulas [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] and [START_REF] Sokolowski | Introduction to shape optimization[END_REF] are not restricted to simply-connected defects, one could parametrize the subsurface heterogeneity as multiple defects, using e.g. description (57) for each. The assumed topology then cannot be altered during the minimization process. As to the correct choice of "initial" topology (e.g. in terms of the number of defects), such preliminary information could be obtained from the available measurements using for instance the methods of topological sensitivity [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][START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF] or linear sampling [START_REF] Baganas | A linear sampling method for the inverse transmission problem in near-field elastodynamics[END_REF][START_REF] Guzina | A linear sampling approach to inverse elastic scattering in piecewise-homogeneous domains[END_REF][START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF].

Minimization

As examined earlier, the goal of this study is the 3D identification of "penetrable" subsurface defects via the minimization of cost functional J ( Ω, Ĉ, ρ)

given by ( 2). Here it is important to remember, however, that certain arguments of J , most notably the material characteristics of the trial defect, are subject to inequality constraints that must be enforced to maintain the physical relevance of the solution. In the context of isotropic elasticity assumed in the ensuing examples, one finds that besides ρ > 0, one must have μ > 0 and -1 < ν < 0.5 to sustain the positive definiteness of the strain energy density [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF]. While numerous techniques are available for non-linear minimization subject to inequality constraints [START_REF] Nocedal | Numerical Optimization[END_REF], most of such algorithms provide only "soft" bounds that can be violated during the minimization process. To aid the strict enforcement of the featured inequality constraints on μ, ν and ρ, the nine-dimensional defect parametrization in (57) is restated using the

transformed variables b = (b 1 , b 2 , . . . , b 9 ) ∈ R 9 where b k = a k , k = 1, . . . 6, b k = log(a k ) k = 7, 9, b 8 = log(0.5 -a 8 ) (60)
which ensure that µ > 0, ν < 0.5, and ρ > 0. The cost function is then expressed naturally in the transformed variables through J b (b) = J a (a), with J a defined by (48). On the basis of (60), the required sensitivities of J b can be computed in terms of ∂J a /∂a k given by ( 38), ( 44) and (49b) as

∂J b ∂b k = ∂J a ∂a k , k ≤ 6, ∂J b ∂b k = ∂J a ∂a k a k k = 7, 9 ∂J b ∂b 8 = ∂J a ∂a 8 (a 8 -0.5).
The minimization problem is posed in a constrained fashion, using parametrization (60), as min

b J b (b), C i (b) ≥ 0 ∀i ∈ I, (61)
where the "soft" inequality constraints C i (I being a set of integers) reflect any additional restrictions in terms of defect's centroidal coordinates (c 1 , c 2 , c 3 ), semi-axes (α 1 , α 2 , α 3 ), and material properties (μ, ν, ρ). The inequality constraints used in the ensuing numerical examples are listed in Appendix A.4.

Following [START_REF] Nocedal | Numerical Optimization[END_REF], optimization problem (61) is for practical implementation reduced, using slack variables, to the unconstrained minimization

min b L A (b, λ m , γ m ) (62)
of an augmented Lagrangian

L A (b, λ m , γ m ) ≡ J b (b) + i∈I ψ(C i (b), λ m i ; γ m ), ψ(C, λ, γ) =    -λ C + C 2 /2γ C ≤ γ λ, -γ λ 2 /2 C > γ λ ( 63
)
whose gradient is computable as

∇ b L A (b, λ m ; γ m ) = ∇L A (b) - i∈I|C i (b)≤γmλ m i (λ m i -C i (b)/γ m ) ∇C i (b). ( 64
)
Given the initial penalty parameter γ 0 > 0, tolerance τ 0 > 0, starting point b 0 , and reference vector of Lagrange multipliers λ 0 , the augmented Lagrangian method introduces a sequence (m = 1, 2, . . .) of unconstrained minimization problems with explicit Lagrange multiplier estimates (λ m ) and decreasing penalties (γ m ) that produce a good estimate of the local KKT (Karush-Kuhn-Tucker) solution, b ⋆ , of (61) even when γ is not particularly close to zero [START_REF] Nocedal | Numerical Optimization[END_REF]. This latter feature is highly desirable as it reduces the possibility of ill-conditioning that commonly occurs for vanishing values of the penalty parameter γ. The algorithm terminates when ∇ b L A < τ ⋆ , where τ ⋆ is the user-chosen ultimate tolerance. For any given iterate m, the nonlinar minimization (62) is effected using the BFGS quasi-Newton method [START_REF] Nocedal | Numerical Optimization[END_REF] (with stopping criterion defined by ∇ b L A < τ m ) and an inexact line search based on the strong Wolfe conditions [START_REF] More | Line search algorithms with guaranteed sufficient decrease[END_REF]. Due to the unconstrained character of minimization subproblems (62), constraints (A.14) are "soft" in that they do not by themselves prevent b from reaching physically inadmissible, or merely undesirable, values (e.g. penny-shaped ellipsoids which may lead to an illconditioned BEM solution). To deal wih the problem, the line search algorithm embedded in (62) has been augmented by a step-reduction feature that prevents b from exceeding the soft-bound limits (A.14) by more than 30%.

For the numerical examples presented next, the internal parameters were set to γ 0 = 0.025, τ 0 = 100, λ 0 = 1, γ m+1 /γ m = τ m+1 /τ m = 0.2, and τ ⋆ = 10 -3 .

Results

The effectiveness of the proposed shape-material sensitivity approach as a tool for reconstructing buried penetrable objects is now demonstrated on a set of numerical results. In all examples to follow, the buried obstacle is "illuminated" using N = 4 × 4 = 16 point forces, sequentially applied at locations x n (n = 1 . . . N) in the ξ 3 -direction over the square testing grid (ξ 1 , ξ 2 ) ∈ [-3d, 3d]×[-3d, 3d] in the ξ 3 = 0 plane, where d is a reference length; for each source, the response of the solid is monitored using M = 5×5 = 25 triaxial receivers with locations x m (m = 1 . . . M) arranged (in the same plane) as shown in Fig. 1. The distance function ϕ(u, u obs , ξ) in ( 2) is taken in the least-squares form and given by

ϕ(u, u obs , ξ) = 1 2 N n=1 M m=1 δ(ξ -x m )|u -u obs | 2 .
Problem quantities are normalized using the length scale d together with the shear modulus and mass density (µ, ρ) of the reference solid; in particular, the nondimensional angular frequency ω = ωd(ρ/µ) 1/2 is introduced. In all examples, the reference solid is additionally characterized by ν = 0.35.

Sensitivity evaluation

To verify the numerical implementation, geometric and material sensitivities ∂J a /∂a k (k = 1, 2, . . . , 9) stemming from [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] and ( 44) are compared with their central difference approximations computed using the boundary integral approach of [START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF] and a surface mesh with 650 eight-noded quadrangular elements. The comparison is performed at a = (.1, .1, 3; .5, .5, .5; 2, .35, .9) for a "true" ellipsoid given by a true = (0, 0, 3; .5, .5, .5; 5, .25, 1), assuming an infinite reference domain Ω = R 3 to ensure maximum accuracy for the Green's functions and focus on the performance of the proposed computation scheme.

The frequency of illumination corresponds to ω = 3. From Table 2, one can see that the relative discrepancy between the sensitivity formulas and their central difference approximations (computed using ±4% perturbation on each parameter) does not exceed 0.4%. It is moreover important to note that the speedup (i.e. the reduction in the computational effort) over the central difference approach is approximately 1/(2(D + 3)) where D + 3 is the total number of (geometric and material) design parameters. This estimate stems from the fact that formulas [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] and ( 44) essentially revolve around the solution of one forward problem (since the adjoint and material-sensitivity problems then exploit the existing matrix factorization), whereas central-difference evaluations entail the set-up and solution of two irreducible problems for each a k .

Obstacle reconstruction

Two examples are now presented to illustrate the reconstruction of penetrable defects in an isotropic, semi-infinite solid via shape sensitivities [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF], material sensitivities [START_REF] Sokolowski | Introduction to shape optimization[END_REF], and the constrained minimization approach described in Section 6.4. In both examples, anticipating the non-convexity characterizing most inverse scattering problems, the initial trial defect is placed relatively close to its target. This assumption is made reasonable by the fact that probing techniques [START_REF] Bao | Inverse scattering by a continuation method with initial guesses from a direct imaging algorithm[END_REF][START_REF] Potthast | A survey on sampling and probe methods for inverse problems[END_REF] based on e.g. topological sensitivity [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][START_REF] Guzina | From imaging to material identification: a generalized concept of topological sensitivity[END_REF] or linear sampling [START_REF] Baganas | A linear sampling method for the inverse transmission problem in near-field elastodynamics[END_REF][START_REF] Guzina | A linear sampling approach to inverse elastic scattering in piecewise-homogeneous domains[END_REF][START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF] provide a reliable preliminary information about the defect location and material characteristics. While these techniques can in principle use the same experimental data exploited by the nonlinear minimization, their explicit coupling with the present scheme is beyond the scope of this study.

Sensitivity ∂ c 1 J ∂ c 2 J ∂ c 3 J ∂ α 1 J ∂ α 2 J ∂ α 3 J ∂ μJ ∂ ν J ∂ ρJ Formulae 0.
Hard obstacle. For this example the testing grid, placed on the surface of the half-space (ξ 3 = 0), and the true defect (indicated as "Hard") are shown in Fig. 1. The frequency of illumination is again such that ω = 3, corresponding to a shear wavelength λ S = 2πd/3. The true defect, centered at (a 1 , a 2 , a 3 ) true = (1, .5, 3), is described as an ellipsoid with semi-axes (a 4 , a 5 , a 6 ) true = (.4, .6, .5) and material properties (μ, ν, ρ) true = (4µ, 0.25, 1.1ρ); the initial iterate is placed at a 0 = (-.4, -.2, 2.5; .5, .5, .5; 1.8, .3, 1), see Fig. 1.

To avoid committing the "inverse crime" [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF] whereby the same model is used to synthesize as well as to invert the data in an inverse problem, synthetic observations u obs are generated using a boundary element mesh with 1,460 nodes, whereas the minimization exploits a coarser mesh with 650 nodes. Fig. 2 illustrates the iterative reconstruction process for this example. To aid the physical insight, the bottom-right panel depicts selected iterations in the 3D space, with the surface color of each iterate (i.e. trial defect) corresponding to its shear modulus according to the attached color bar. As can be seen from the display, the solution converges to the global minimum of J after approximately 70 iterations. Not surprisingly, the centroidal coordinates exhibit the fastest convergence, followed by that in terms of the semi-axes and material properties. It should be noted, however, that the synthetic data in this example contain no extraneous perturbations other than those caused by the use of dissimilar BEM meshes. To examine the effect of measurement uncertainties, synthetic observations u obs are next corrupted as ũobs := (1 + ̺χ) u obs (65) over all source-receiver pairs, where ̺ is the noise amplitude and χ ∈ [-1, 1] is a uniform random variable. Table 3 lists the reconstructed defect parameters for the noise amplitude levels of 0, 1, and 2%. The defect reconstruction is seen to be fairly sensitive to measurement noise. For instance, a ̺ = 1% perturbation of experimental data leads to average (absolute) errors of 0.4%, 1.6% and 15% in terms of the defect's centroid, semi-axes, and material properties. While the latter figure may seem excessive, one may recall that perturbation (65) is specified in terms of the total field u. As a point of reference, the induced perturbation in terms of the scattered field u S -u F , which carries all available information about the defect, exceeds 50% for selected source-receiver pairs when ̺ = 1%. In practical situations, this problem may be mitigated using a combination of i) multi-tonal illumination, in combination with accelerated BEM such as the Fast multipole method allowing for higher frequencies in forward simulations [START_REF] Chaillat | A multi-level fast multipole BEM for 3-D elastodynamics in the frequency domain[END_REF][START_REF] Nishimura | Fast multipole accelerated boundary integral equation methods[END_REF], and ii) Bayesian (e.g. maximum likelihood) data analysis, where prior information on the sought inclusion and measurement errors are incorporated into a posterior probability density function [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF].

Soft obstacle. To illustrate the method in cases where the inclusion is more compliant than the background material, this example considers a true defect with centroidal coordinates (a 1 , a 2 , a 3 ) true = (-1.1, -0.3, 2.7), semi-axes (a 4 , a 5 , a 6 ) true = (.6, .8, .6) and material properties (μ, ν, ρ) true = (0.2µ, 0.1, 0.8ρ). The testing arrangement, mimicking that in the previous example, as well as the true ("Soft") defect are again shown in Fig. 1. Each of the 16 sequentiallyapplied point sources acts in the ξ 3 -direction at frequency ω = 2, corresponding to a shear wavelength λ S = πd. The meshes used for the generation of synthetic measurements and the reconstruction respectively feature 650 and 290 nodes. Figure 3 illustrates the reconstruction of the soft defect, assuming a 0 = (-.4, -.2, 2.5; .5, .5, .5; .8, .3, 1) for the initial iterate whose geometry is shown in Fig. 1. Again, the bottom right-hand panel plots the selected iterations, color-coded according to the trial shear modulus μ of each iterate. In this case, the solution converges to the global minimum of J after approximately 50 iterations, with centroidal coordinates leading the way. With reference to the bottom left-hand panel, a notable feature of the minimization process is an activated inequality constraint ν ≤ 0.49, see (A.14), which steers the numerical solution safely away from the incompressible case ν = 0.5.

Conclusions

In this work, 3D inverse scattering of elastic waves involving penetrable solid defects is investigated within the framework of a boundary integral equation method. The inverse problem is reduced to the minimization of a misfit between experimental observations and their simulations for a trial inclusion. To maximize the accuracy and efficiency of gradient-based search algorithms, the shape sensitivity of the cost function is formulated via an adjoint problem approach, extending earlier works on void identification. This is complemented by a novel material sensitivity formulation, developed using two alternative methodologies, namely the direct-differentiation and adjoint-field approaches. The proposed shape and material sensitivity formulas, computable as surface integrals over the trial defect boundary, are implemented and incorporated into a nonlinear optimization algorithm based on an augmented Lagrangian that facilitates the imposition of inequality constraints. From preliminary numerical studies, the latter were found to be critical in avoiding physically inadmissible or computationally inadequate trial inclusion configurations. The effectiveness of the proposed sensitivity formulation is demonstrated on numerical examples dealing with the reconstruction of an ellipsoidal inclusion embedded in a semiinfinite solid. Future work towards a comprehensive computational platform for elastic-wave imaging of penetrable defects will incorporate other necessary components: preliminary defect-indicator function based on e.g. topological sensitivity or linear sampling methods, refined misfit functions based on e.g. Bayesian concepts that allow for stochastic data analysis, and a fast multipole (accelerated) version of the elastodynamic boundary element method. Leading singularity of kernel sensitivities. It is useful to investigate the leading contributions to kernel sensitivities (A.12) for r → 0, e.g. for the purpose of handling singular element integrals. Noting that the corresponding leading contributions for the elastodynamic kernels (A.8a-b) themselves are the corresponding components of the elastostatic Kelvin solution, given by [ Ûiℓ ] 1 = 1 8π μr (1 + γ2 )δ iℓ + (1 -γ2 )r ,i r ,ℓ , (A.13a)

[ Tiℓ ] 1 = 1 4πr 2 3(γ 2 -1)r ,i r ,ℓ r ,j n j + γ2 (r ,ℓ n iδ iℓ r ,j n jr ,i n ℓ ) , (A.13b) the leading contributions to (A.12) for r → 0 are easily obtained from

A.4 Inequality constraints

With reference to parametric descriptions (57) and (60) and the constrained minimization problem (61), a lower and an upper bound is imposed on each parameter b i (i = 1, 2, . . . . 9), resulting in the total of 18 constraints C i (b). On noting that the size of the square testing area S obs in Fig. 1 is 4d × 4d, these inequality constraints are specified so that where µ and ρ are the shear modulus and the mass density of the reference (i.e. background) solid. Physically, the restrictions on c i require that the centroid of the defect is located at least partially "under" the testing area, that the defect is physically separated from the surface of the half-space, and that the maximum search depth be commensurate with the size of the testing grid; the bounds on the ellipsoid's semi-axes are imposed to avoid overly distorted shapes, while the additional restrictions on material parameters are used to i) prevent ill-conditioning of the numerical solution (e.g. in terms of excessively small values of μ), and ii) focus the search on the range of expected values.

Fig. 1 .

 1 Fig. 1. Testing grid and sample defects in a semi-infinite solid ξ 3 > 0 (ν = 0.35)

2 .

 2 Reconstruction of a hard defect (ω = 3, μtrue = 4µ).

  Fig.

3 .

 3 Reconstruction of a soft defect (ω = 2, μtrue = 0.2µ).

(

  where notation ∂ k indicates partial differentiation w.r.t. k) one readily obtains∂ k G ,i (r; k) = -krr ,i G(r; k) ∂ k G ,ij (r; k) = -k rr ,i G ,j (r; k) + δ ij G(r; k) ∂ k G ,ijℓ (r; k) = -k rr ,i G ,jk (r; k) + δ ij G ,k (r; k) + δ ik G ,j (r; k)With the help of the above identities, the partial derivatives of the Green's function (A.8a-b) with respect to the wavenumbers are found to be given byμk S ∂ kS Ûiℓ = (1 + i kS r)r ,i r ,ℓ + (1i kS r)δ iℓ G(r; kS ) -2μ Ûiℓ (A.11a) μk P ∂ kP Ûiℓ = -γ 2 (i kP r + 1)r ,i r ,ℓδ iℓ G(r; kP ) (A.11b) kS ∂ kS Tiℓ = -2 Tiℓ + 2rn i G ,ℓ (r; kP ) -2rr ,ℓ n j G ,ij (r; kS ) -k2 S r δ ℓi r ,n + r ,i n ℓ G(r; kS ) (A.11c) kP ∂ kP Tiℓ = γ2 2rr ,ℓ n j G ,ij (r; kP ) + (2 k2 P -k2 S )rr ,ℓ n i G(r; kP ) + 2r n ℓ G ,i (r; kP )n i G ,ℓ (r; kP ) (A.11d)with γ again given by (A.9), and having set r ,n = r ,j n ,j . Then, noting that Û and T with respect to the material parameters of the inclusion are given in terms of expressions (A.11a-d) byd Ûiℓ dρ = 1 2ρ kP ∂ kP + kS ∂ kS Ûiℓ d Tiℓ dρ = 1 2ρ kP ∂ kP + kS ∂ kS Tiℓ d Ûiℓ dμ = -1 2μ 2 + kP ∂ kP + kS ∂ kS Ûiℓ d Tiℓ dμ = -1 2μ kP ∂ kP + kS ∂ kS Tiℓ d Ûiℓ dν = -γ 2 kP (γ 2 -1) 2 ∂ kP Ûiℓ d Tiℓ dν = -γ 2 kP (γ 2 -1) 2 ∂ kP TiℓThe total material-parameter sensitivities of kernels Û ′ and T ′ used in (43) are finally given in terms of the above expressions by Û ′ = ρ′ d dρ Û

  1ν) 2 r (r ,i r ,ℓδ iℓ ) + O(1), d Tiℓ dν = 1 8π(1ν) 2 r 2 δ ℓi r ,n + r ,i n ℓr ,ℓ n i -3r ,i r ,ℓ r ,n + O(1),which coincide, as expected, with the corresponding material sensitivities of the Kelvin solution (A.13a-b).

d 20 ≤ c 3 -

 203 α 3 , c 3 ≤ 5d, Semi-axes: d 20 ≤ α i ≤ d (i = 1, 2, 3), (A.14)Material properties:µ 100 ≤ μ ≤ 10µ, 0.01 ≤ 0.5ν ≤ 0.5, ρ 100 ≤ ρ ≤ 3ρ,

Table 1

 1 Sample CPU times [min:sec] per evaluation the cost function and its sensitivities.

	n p	2	4	8	12	24	48	96
	CPU time 46:01 23:27 12:41 9:05	5:52	3:24	2:31
	Q s	1.98	3.89	7.20 10.05 17.01 26.85 36.28
	Q s /n p	0.99	0.97	0.90	0.84	0.70	0.56	0.38
	joint, and material-sensitivity problems in (53a-c) is implemented, together
	with formulas (38) and (

Table 2

 2 Sensitivities ∂ a k J ≡ ∂J /∂a k : comparison with central differences

Table 3

 3 Sensitivity of the solution to experimental noise (hard defect)

	Parameter	c 1 /d c 2 /d c 3 /d α 1 /d α 2 /d α 3 /d μ/µ	ν	ρ/µ
	True	1.0	0.50	3.0	0.40	0.60	0.50	5.0	0.25	1.1
	Identified, ̺ = 0%	1.0	0.50	3.0	0.40	0.60	0.50	5.0	0.25	1.1
	Identified, ̺ = 1% 0.996 0.496 3.0 0.408 0.607 0.508 3.76 0.297 1.1
	Identified, ̺ = 2% 0.992 0.492 2.99 0.367 0.583 0.481 5.19 0.016 1.06
				21

T (x, ξ) denote the (full-space) elastodynamic Green's tensors corresponding to the material properties of the inclusion. With these definitions, the forward problem (4)-(6) for a semi-infinite solid can be reformulated in terms of a pair of regularized boundary integral equations[START_REF] Bonnet | Boundary Integral Equations Methods for Solids and Fluids[END_REF][START_REF] Pak | Seismic soil-structure interaction analysis by direct boundary element methods[END_REF]:u(x) + Γ u(ξ)u(x) •[T (x, ξ)] 1 dS ξ + Γ u(ξ)•[T (x, ξ)] 2 dS ξ -Γ t(ξ)•U (x, ξ) dS ξ = u F (x),x ∈ Γ,[START_REF] Bonnet | BIE and material differentiation applied to the formulation of obstacle inverse problems[END_REF]

A Appendices

A.1 Proof of identity [START_REF] Lee | Identification of geometric shapes and material properties of inclusions in two-dimensional finite bodies by boundary parameterization[END_REF] From [START_REF] Guzina | On the analysis of wave motions in a multi-layered solid[END_REF] and the symmetry properties of the elastic tensor C, one obtains with E(•, •) defined by [START_REF] Madyarov | A Radiation condition for layered elastic media[END_REF]. Moreover, one easily checks by that the following identity holds for any pair of (sufficiently smooth) fields u, w:

On using (A.2), applying the divergence theorem together with (17c), and invoking field equation (4a), equation (A.1) becomes

Applying the same treatment to the second integral of (8b) yields a variant of (A.3) wherein Ω -, n, u, w, C, ρ are respectively replaced with Ω, n, û, ŵ, Ĉ, ρ. Expression [START_REF] Lee | Identification of geometric shapes and material properties of inclusions in two-dimensional finite bodies by boundary parameterization[END_REF] then follows from applying (A.3) and its counterpart to (8b) and invoking definition (8b) for interpreting terms featuring a(•, •) or â(•, •).

A.2 Generalization of (30) to arbitrary free fields

The proof of sensitivity result [START_REF] Pelekanos | Contrast source inversion methods in elastodynamics[END_REF] given in Section 4 assumes that either (i) Ω is bounded or (ii) the free field satisfies the generalized radiation condition at infinity. Assumption (ii) can be relaxed so that any free-field u F satisfying Lu F = 0 in Ω (e.g. a plane wave) is permitted. This result rests upon a modified version of the weak formulation ( 7) which reads

where u S = u-u F is the scattered field in Ω -. A derivation along the lines of Appendix (A.1) yields the Lagrangian-derivative form of equation (A.4) as

Identity (A.5) is rearranged by noting that the free field satisfies

(with O as introduced in ( 17)) which, upon carrying out the Lagrangian differentiation, yields

The free field, being insensitive to the inclusion shape and properties, has a Lagrangian derivative given by

Substituting (A.6) into (A.5) and using (A.7) with the divergence theorem, the counterpart of the weak shape sensitivity formulation ( 30) is found as

From that point, the adjoint solution is again defined by [START_REF] Nintcheu Fata | A linear sampling method for near-field inverse problems in elastodynamics[END_REF]. Proceeding as in Sec. 4, one finds that the shape sensitivity formula (30) still holds, in particular because the last integral in (A.5) vanishes by virtue of conditions (32f).

A.3 Material sensitivity of elastodynamic fundamental solution

The full-space elastodynamic Green's tensors Û and T are given, for the inclusion medium, by their components (the minus sign in the exponential being consistent with the assumed implicit time-harmonic factor exp(iωt)). Noting that