Guillaume Bal

Cédric Bellis
email: bellis@lma.cnrs-mrs.fr

Sébastien Imperiale
email: sebastien.imperiale@inria.fr

François Monard
email: fmonard@math.washington.edu

Reconstruction of constitutive parameters in isotropic linear elasticity from noisy full-field measurements

Keywords: Quantitative parameter identification, Internal data, Regularization method, Elastography, Lamé system

Within the framework of linear elasticity we assume the availability of internal full-field measurements of the continuum deformations of a non-homogeneous isotropic solid. The aim is the quantitative reconstruction of the associated moduli. A simple gradient system for the sought constitutive parameters is derived algebraically from the momentum equation, whose coefficients are expressed in terms of the measured displacement fields and their spatial derivatives. Direct integration of this system is discussed to finally demonstrate the inexpediency of such an approach when dealing with noisy data. Upon using polluted measurements, an alternative variational formulation is deployed to invert for the physical parameters. Analysis of this latter inversion procedure provides existence and uniqueness results while the reconstruction stability with respect to the measurements is investigated. As the inversion procedure requires differentiating the measurements twice, a numerical differentiation scheme based on an ad hoc regularization then allows an optimally stable reconstruction of the sought moduli. Numerical results are included to illustrate and assess the performance of the overall approach.

Introduction

The identification of constitutive parameters associated with deformable solids bears relevance to a wide range of applications [START_REF] Bonnet | Inverse problems in elasticity[END_REF] such as structural health monitoring, nondestructive material testing or medical imaging. A static or dynamical excitation applied to an elastic body gives rise to a measurable signature of such quantities of interest. The data, possibly noisy, might be then collected externally, i.e. at the boundary of the domain, or internally in the form of partial or full-field measurements. Classical identification methods aiming at the quantitative reconstruction of parameters revolve around the minimization of a given objective function. Such approaches typically suffer from requiring numerous forward solutions, and from an intrinsic ill-posedness yielding practical instabilities as discussed in e.g. [START_REF] Monard | Taming unstable inverse problems: Mathematical routes toward high-resolution medical imaging modalities[END_REF].

State-of-the-art experimental techniques [START_REF] Greenleaf | Selected methods for imaging elastic properties of biological tissues[END_REF][START_REF] Parker | Imaging the elastic properties of tissue: the 20 year perspective[END_REF][START_REF]Full-Field Measurements and Identification in Solid Mechanics[END_REF][START_REF] Buffière | Imagerie 3D en mécanique des matériaux[END_REF] offer a variety of non-invasive imaging modalities providing internal full-field measurements for (i) biological tissues using, e.g., ultrasound, magnetic resonance imaging or speckle interferometry, and (ii) materials by X-ray, neutron diffraction tomography or digital image correlation. Dedicated identification methods are also flourishing [START_REF] Parker | Imaging the elastic properties of tissue: the 20 year perspective[END_REF][START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF]. The breakthrough is the availability of internal data, which has led to a paradigm shift in imaging: the constitutive parameters entering the model partial differential equation can now be reconstructed given the knowledge of one, or multiple, solutions of that PDE. In other words, the momentum equation is seen as a PDE for the unknown moduli with the measured displacement or strain fields constituting its spatially varying coefficients. This is the perspective adopted in [START_REF] Mclaughlin | Unique identifiability of elastic parameters from time-dependent interior displacement measurement[END_REF] based on time-dependent data. In a similar context, the article [START_REF] Barbone | Elastic modulus imaging: some exact solutions of the compressible elastography inverse problem[END_REF] focuses on the case when one Lamé parameter is known a priori and only one measured displacement field is used to reconstruct the other one. In [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF], a dedicated variational formulation is employed to solve for the elasticity parameters using static displacement field measurements.

The approach adopted here is to be linked to the so-called hybrid or multi-wave inverse problems [START_REF] Bal | Hybrid inverse problems and redundant systems of partial differential equations[END_REF][START_REF]Hybrid inverse problems and internal functionals[END_REF], which under the assumption of a coupling between two physical phenomena, come as two-step inverse problems. A typical scenario is that a highresolution imaging modality will provide internal measurements (e.g. the displacement fields, here) for a parameter-reconstruction problem involving an elliptic PDE. While this first step is assumed here, our starting point is the second step of this inversion. For such inverse problems where the forward model is an elliptic PDE, considering inversion from internal measurements (rather than, classically, measurements at the domain's boundary) greatly improve the mathematical ill-posedness, which practically implies great improvements in resolution on reconstructed images. As the problem remains mildly ill-posed, noisy measurements still require being regularized before the inversion step in order not to amplify high-frequency noise [START_REF] Engl | Regularization of Inverse Problems[END_REF], though this regularization is much less stringent than when using boundary measurements, allowing recovery of information at smaller scales.

This study lies at the crossroads of experimental solid mechanics and hybrid inverse problems. It is similar in spirit to the so-called adjoint-weighted variational formulations introduced in [START_REF] Barbone | Adjointweighted variational formulation for the direct solution of plane stress inverse elasticity problems[END_REF][START_REF] Albocher | Adjoint-weighted equation for inverse problems of incompressible plane-stress elasticity[END_REF][START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF]. However, advantage is taken here of an explicit formulation of a gradient system for the sought parameters. This approach stems from the method in [START_REF]Hybrid inverse problems and internal functionals[END_REF][START_REF] Bal | Reconstruction of coefficients in scalar second-order elliptic equations from knowledge of their solutions[END_REF][START_REF] Monard | Inverse diffusion problems with redundant internal information[END_REF] for scalar diffusion equations, which is therefore extended here to the tensorial framework of elasticity. This enables further mathematical characterizations of the reconstruction stability and an inversion from a least squares approach. A regularization-based data differentiation scheme is also proposed to accommodate noisy measurements.

The article outline is as follows: (i) The governing equations of linear elasticity are presented in Section 2. Section 3 introduces the inverse problem of the quantitative reconstruction of the spatiallydependent constitutive moduli, namely the two eigenvalues of the elasticity tensor. Assuming the availability of a set of measured displacement fields, a system of PDEs for the unknown parameters is constructed by algebraic manipulations. Upon satisfying an invertibility condition, these equations are then recast into a simple gradient system. (ii) In Section 4, this gradient system is discussed using a conventional ordinary differential equation-based approach which yields a local uniqueness and stability result. (iii) Alternatively, when the available data is noisy, then a weak formulation of the least squares problem associated with the system at hand is presented in Section 5. Existence and uniqueness of a solution is finally shown using a standard analysis. (iv) Strain and hessian tensors of measured displacement fields that enter the gradient system are clearly strongly detrimental to the inversion in the presence of noise. In this regard, the reconstruction stability in such configurations is analyzed in Section 6. (v) Finally, a numerical differentiation scheme is proposed and analyzed in Section 7. It relies on a regularization of the measurements by L 2 -projections on coarse, yet high-order, finite element spaces. The regularizing operator is constructed to enable an optimally stable reconstruction of the constitutive moduli. A set of numerical results is included in Section 8.

Preliminaries

Let Ω ⊂ R d with d = 2 or 3, denote a regular enough elastic body which undergoes a time-harmonic infinitesimal transformation characterized by the displacement field u and the frequency ω. On noting ∂ j the j-th partial derivative and assuming that |∂ j u i | = o(1) for 1 ≤ i, j ≤ d, the corresponding linearized deformation is quantified by the second-order strain tensor ε = ∇ s u := 1 2 [∇u + (∇u) T ]. In the absence of any body force field, the displacement field satisfies the momentum equation [START_REF] Marsden | Mathematical foundations of elasticity[END_REF]

∇ • (C : ε) + ρ ω 2 u = 0 in Ω, (1)
where the mass density ρ > 0 is assumed to be constant and with reference to the implicit time-harmonic factor e iωt . In the isotropic case considered, the elasticity tensor C is characterized by only two eigenvalues α and β, see [START_REF] Knowles | On the representation of the elasticity tensor for isotropic media[END_REF], such that

C = α d I ⊗ I + β I sym - 1 d I ⊗ I , (2)
where I and I sym are respectively the second-order and symmetric fourth-order identity tensors. In the ensuing analysis it will be seen that inverting for the elastic moduli α > 0 and β > 0, rather than for the Lamé parameters λ = (α-β)/d and µ = β/2, is facilitated by the decomposition and interpretation of the measured strain fields into elementary hydrostatic and deviatoric contributions. The analysis of the proposed reconstruction approach further requires that (α, β) are smooth enough, i.e.

Assumptions 1. The constitutive parameters satisfy:

α, β ∈ L ∞ (Ω) ∩ H 1 (Ω).
The boundedness assumption is a direct consequence of standard thermo-mechanical stability conditions [START_REF] Mehrabadi | Strain energy density bounds for linear anisotropic elastic materials[END_REF]. The requirement α, β ∈ H 1 (Ω) is more specifically associated with the H 1 (Ω)-norm that is employed hereinafter to measure the parameters reconstruction quality. Yet, some regularity on the solution of ( 1) is required, in general, and this is to be achieved with sufficiently smooth coefficients.

To be used in the ensuing developments, for p = 2 or p = ∞ with associated Lebesgue space L p (Ω), and (m 1 , m 2 ) ∈ N 2 , let us introduce a norm associated with matrices

H ∈ L p (Ω) m 1 ×m 2 as H L p := |H| L p , (3)
which is expressed in term of the Frobenius norm |H| = (H : H) 1/2 with double inner product ":". Moreover, given ≥ 0, Sobolev spaces W ,p (Ω) and H (Ω) = W ,2 (Ω), then for all vectors h ∈ W ,p (Ω) m 1 , the definition (3) remains valid with m 2 = 1 when = 0 by replacing the Frobenius norm by the Euclidean norm |h| = (h • h) 1/2 that uses the single inner product. When ≥ 1, given an orthonormal basis {e i } i of R d , one defines

h H = m 1 i=1 h • e i 2 H 1 2 and h W ,∞ = m 1 i=1 h • e i W ,∞ .

Reconstruction of elastic moduli

This study focuses on the reconstruction of the spatially varying elastic moduli α, β from the knowledge of two displacement fields u 1 , u 2 within Ω. The choice of this measurement set is supported by earlier studies [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF][START_REF] Bal | Hybrid inverse problems and redundant systems of partial differential equations[END_REF][START_REF]Hybrid inverse problems and internal functionals[END_REF]. For n = 1, 2, one denotes strain tensors by ε n = ∇ s u n with trace and deviatoric counterpart given by

t n = tr(ε n ), ε D n = ε n - t n d I. (4)

Inversion formula

Plugging the decomposition (2) into the momentum equation (1) one obtains

∇ • (C : ε n ) + ρ ω 2 n u n = ∇ • α t n d I + βε D n + ρ ω 2 n u n = 0, n = 1, 2, (5)
which, expanding using the product rule, yields

t n d ∇α + α∇ t n d + ε D n • ∇β + β∇ • ε D n + ρ ω 2 n u n = 0, n = 1, 2. (6)
On introducing the matrices

A n = tn d I ε D n ∈ R d×2d and B n = ∇ tn d ∇ • ε D n ∈ R d×2 , (7)
then Eqn. (6) reads

A n ∇α ∇β + B n α β = -ρ ω 2 n u n , n = 1, 2.
We then combine these equations into an overdetermined PDE system of the form

A ∇α ∇β +B α β = - ρ ω 2 1 u 1 ρ ω 2 2 u 2
, where

A := A 1 A 2 and B := B 1 B 2 . (8)
We now make more precise when and how to invert the 2d × 2d matrix-valued A.

Lemma 1. Let us define

:= t 1 ε D 2 -t 2 ε D 1 ( 9
)
Suppose that is uniformly invertible, then A defined in ( 8) is uniformly invertible and

A -1 = d ε D 2 -d ε D 1 -t 2 I t 1 I -1 0 0 -1 .
Using Lemma 1 and upon inverting the block matrix A, one arrives at the system

(α, β) = f , where (α, β) := ∇α ∇β + M α β (10)
with

M = a b c d := A -1 B and f = f 1 f 2 := -A -1 ρ ω 2 1 u 1 ρ ω 2 2 u 2 . ( 11
)
When M ∈ L ∞ (Ω) 2d×2 , the operator is linear with values from H 1 (Ω) 2 into L 2 (Ω) 2d . Note that such a system is valid both in 2D and 3D upon substitution of the correct expression for the strain tensors and their traces. Moreover, one is required to verify that the tensor is invertible.

The elasticity parameters can be reconstructed from Eqn. [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF] either by direct integration of this gradient system, which requires some compatibility conditions, or using a least squares based variational formulation. These two approaches will be addressed in sections 4 and 5 respectively. In any case, one can already deduce some regularity requirement so that (10) make sense in a standard functional sense for which M ∈ L ∞ (Ω) 2d×2 and f ∈ L ∞ (Ω) 2d . These hypotheses are summarized below:

Assumptions 2. The displacement field solutions satisfy u n ∈ H 2 (Ω) d ∩ W 1,∞ (Ω) d for n = 1, 2 and inf Ω (| det |) ≥ c 0 > 0.
Remark 1. It is generally assumed that incompressible materials, such as biological soft tissues [START_REF] Fung | Foundations of solid mechanics[END_REF] can be thought of as the limit of the compressible case in the limit α → ∞. In this case, tr(ε) → 0 which implies det → 0 so that Assumptions 2 are violated. Therefore, the compressible and the incompressible cases are very dissimilar, the latter being beyond the scope of this article. Reference to [START_REF] Barbone | Elastic modulus imaging: some exact solutions of the compressible elastography inverse problem[END_REF], Section 8, can be made for a discussion on the discrepancy between these two cases.

Invertibility conditions

As seen in the section above, the reconstruction procedure requires the tensor to be invertible. This section now aims at clarifying this assumption.

Proposition 1. The tensor := t 1 ε D 2 -t 2 ε D 1 is non-invertible if and only if the strain tensors satisfy t 1 ε 2 = t 2 ε 1 when d = 2, or t 1 ε 2 = t 2 ε 1 + γ(φ 1 ⊗ φ 1 -φ 2 ⊗ φ 2 ) when d = 3
for some scalar γ and two unit orthogonal vectors φ 1 , φ 2 .

In particular, is non-invertible if, given strain tensor ε 1 , there exists

κ ∈ R \ {0} such that ε 2 satisfies ε 2 = κ ε 1 (2D) or ε 2 = κ ε 1 + t -1 1 γ(φ 1 ⊗ φ 1 -φ 2 ⊗ φ 2 ) (3D).
Proof of Proposition 1. By construction, the tensor

= t 1 ε D 2 -t 2 ε D 1 = t 1 ε 2 -t 2 ε
1 is symmetric and traceless. In two dimensions, this implies that it takes the form = a b b -a with a, b two real parameters, in which case det = -(a 2 + b 2 ), so det = 0 is equivalent to = 0. In the case d = 3 one has sp( ) = {γ 1 , γ 2 , -(γ 1 + γ 2 )}. with γ 1 , γ 2 ∈ R. The condition det( ) = 0 entails γ i = 0 for i ∈ {1, 2} or γ 2 = -γ 1 so that one can conclude to the existence of an orthonormal basis

{φ i } i of R d such that = γ(φ 1 ⊗ φ 1 -φ 2 ⊗ φ 2 ).
How to achieve condition det = 0 ? Solutions u 1 , u 2 have to be constructed such that det = 0 is achieved locally or globally. For example, a strategy is to generate displacement fields satisfying the properties that (i) t 1 = 0, (ii) ε 1 invertible and (iii) t 2 = 0, since in this case the tensor takes the form = -t 2 ε 1 . From the practical standpoint, the remaining question concerns the control of these solutions from the domain boundary as discussed below.

Constant coefficients: A simple prototypical example in the case of constant coefficients and stationary regime is to choose, in dimension d = 2 or d = 3,

u 1 (x) • e j = i =j x i , 1 ≤ j ≤ d,
and

u 2 (x) = x. (12)
In this case, we have that ε 1 = U -I, where U is an d × d matrix filled with ones, and ε 2 = I. Here we globally have

t 1 = 0, t 2 = d = 0 and det ε 1 = (-1) d-1 (d -1) = 0, so that = t 1 ε 2 -t 2 ε
1 is nowhere singular. Note that leeway is allowed in choosing admissible fields u 1 , u 2 and known analytical solutions of elasticity can be employed.

Non-constant coefficients: In the case of non-constant coefficients, it is no longer clear theoretically that one can achieve the condition det = 0 globally with only two solutions, although numerical results in Section 8 will show that numerically, two solutions, chosen after the example above for instance, are usually enough to obtain satisfactory global reconstructions.

Following ideas initiated by Bal and Uhlmann in [START_REF] Bal | Reconstruction of coefficients in scalar second-order elliptic equations from knowledge of their solutions[END_REF] (see also [START_REF] Bal | Inverse anisotropic conductivity from internal current densities[END_REF][START_REF] Monard | Inverse anisotropic conductivity from power densities in dimension n ≥ 3[END_REF][START_REF] Bal | Boundary control of elliptic solutions to enforce local constraints[END_REF]), let us mention that under some regularity assumptions, it can be proven that there exist boundary conditions such that condition det = 0 holds locally. As seen in the scalar case, such proof will not construct the boundary conditions explicitly. Yet it justifies that the nonvanishing determinant condition can be achieved and provides a way to prove local reconstructibility of parameters from full-field measurements. The main tool for such proofs is the Runge approximation property, which itself follows from unique continuation property, established in e.g. [START_REF] Nakamura | Oscillating-decaying solutions, Runge approximation and its applications to inverse problems[END_REF][START_REF] Lin | Conditional stability for the hexagonal anisotropic elastic dynamical systems[END_REF][START_REF] Lin | Quantitative strong unique continuation for the Lamé system with less regular coefficients[END_REF] in the context of elasticity.

ODE-based approach and stability estimate

In this section we consider a direct integration approach of the gradient system [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF], its main interest (stability estimate) and its main shortcoming (dependence on integration curve). Considering a connected subset Ω 0 ⊆ Ω and two points x 0 , x ∈ Ω 0 , let

γ : [0, 1] t → ξ(t) such that ξ(0) = x 0 and ξ(1) = x,
be a smooth curve with endpoints x 0 and x. Using the chain rule along this curve, the gradient system (10) yields an ODE along the curve γ, of the form

dϕ γ (t) dt +M γ (t)ϕ γ (t) = f γ (t) with ϕ γ (t) := α • γ(t) β • γ(t) and ϕ γ (0) = α(x 0 ) β(x 0 ) (13)
and where one has defined, for t ∈ [0, 1]

M γ (t) := ξ(t) • a • γ(t) ξ(t) • b • γ(t) ξ(t) • c • γ(t) ξ(t) • d • γ(t) and f γ (t) := ξ(t) • f 1 • γ(t) ξ(t) • f 2 • γ(t) .
From the values (α(x 0 ), β(x 0 )), the values (α(x), β(x)) can be computed via direct integration of ODE [START_REF] Engl | Regularization of Inverse Problems[END_REF]. Fixing x 0 and varying x throughout Ω 0 , this induces a reconstruction procedure of (α, β) throughout Ω 0 , since Ω 0 is connected, every x ∈ Ω 0 can be connected to x 0 via a smooth curve. This approach yields a unique and stable reconstruction of the sought moduli in the sense of the following theorem.

Theorem 1. Consider a convex subset Ω 0 ⊆ Ω and two data sets (u 1 , u 2 ), (u 1 , u 2 ) with C 2 -smooth components, associated to two set of parameters (α, β) and (α , β ) and satisfying inf Ω 0 (| det |, | det |) ≥ c 0 > 0.
Then the data sets determine uniquely the parameters (α, β) and (α , β ) over Ω 0 and there exists a positive constant C such that

α -α W 1,∞ (Ω 0 ) + β -β W 1,∞ (Ω 0 ) ≤ C 0 + 2 n=1 u n -u n W 2,∞ (Ω 0 ) , where 0 = |α(x 0 ) -α (x 0 )| + |β(x 0 ) -β (x 0 )| is the error committed at x 0 ∈ Ω 0 .
The proof of Theorem 1 relies on the Picard-Lindelöf theorem and Gronwall's lemma to control the propagation of errors along characteristic curves, similar analysis can be found in coupled-physics contexts in e.g. [START_REF] Monard | Inverse diffusion problems with redundant internal information[END_REF] and is omitted here. While this approach provides us with an explicit reconstruction procedure yielding unique and stable reconstructions, this uniqueness depends on the choice of curve joining x 0 and x. This is because the gradient system (10) is redundant and requires that the measurements fulfill additional integrability conditions in order for the reconstruction to not depend on this choice of integration curve. Noise in the measurements will however violate these conditions, thereby making this reconstruction approach somewhat not as well-adapted to noisy measurements as the variational formulation that follows.

Variational formulation

Hereafter, the regularity requirements in Theorem 1 are relaxed and Assumptions 2 are considered in the ensuing developments.

Noisy data and least squares approach

Following the above discussion, from now on, let {(u 1,δ , u 2,δ )} denote a set of noisy measurements where the real-valued index parameter δ ≥ 0 is intended to represent a measure of the noise. In such a case, direct differentiations of these measurements, as required in [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF], might amplify noise at high spatial frequencies thereby preventing a successful identification. Therefore, the first step towards the formulation of a reconstruction algorithm is to introduce some approximations, in a sense that will be specified later on, of the strain and hessian tensors of the noisy displacements fields as

ε n,δ ∼ ε n and H n,δ ∼ d k=1 ∇ ⊗2 (u n,δ • e k ) ⊗ e k for n = 1, 2. ( 14
)
with the short-hand notation

H k n,δ = H n,δ • e k for k = 1, . . . , d so that (H k n,δ ) ij ∼ ∂ i ∂ j (u n,δ • e k ). (15)
The explicit constructions of the second and third-order tensors ε 1,δ , ε 2,δ and H 1,δ , H 2,δ from the noisy data (u 1,δ , u 2,δ ) is the focus of Section 7. For now, we consider the noisy operator δ , featuring the noisy matrix M δ with components a δ , b δ , c δ and d δ , together with the noisy source term f δ that are formally obtained upon substitutions of

u n ← u n,δ , t n ← t n,δ = tr(ε n,δ ), ε D n ← ε D n,δ = ε n,δ - t n,δ d I, ← δ := t 1,δ ε D 2,δ -t 2,δ ε D 1,δ , ∇t n ← d i,j=1 (H j n,δ ) ij e i , ∇ • ε n ← d i,j=1 d -2 2 d (H j n,δ ) ij + 1 2 (H i n,δ ) jj e i .
for n = 1, 2 into the equations from ( 7) to [START_REF] Bal | Hybrid inverse problems and redundant systems of partial differential equations[END_REF]. With the noisy operator δ at hand, the problem ( 10) is correspondingly recast as the minimization problem inf

(α δ ,β δ ) δ (α δ , β δ ) -f δ 2 L 2 (Ω) , for all δ > 0. ( 16
)
When δ = 0, we define 0 ≡ and f 0 ≡ f with the corresponding solution to ( 16) satisfying (α 0 , β 0 ) ≡ (α, β) and it will be proven next by convergence analysis that this solution coincides with the limit solution to the perturbed problem series [START_REF] Bal | Reconstruction of coefficients in scalar second-order elliptic equations from knowledge of their solutions[END_REF] as δ → 0.

We first use a lifting (α , β ) of the assumed-to-be-known boundary values of (α, β)

α δ = α o δ + α and β δ = β o δ + β , (17)
to pose a problem in (α o δ , β o δ ) ∈ H 1 0 (Ω) 2 . In order to compute the solution of the minimization problem [START_REF] Bal | Reconstruction of coefficients in scalar second-order elliptic equations from knowledge of their solutions[END_REF], see [START_REF] Bochev | Least-Squares Finite Element Methods[END_REF], we consider solving the associated equation * δ

δ (α o δ , β o δ ) = * δ f δ -δ (α , β ) , (18)
where

* δ : L 2 (Ω) 2d → H -1 (Ω) 2 denotes the adjoint of δ : H 1 0 (Ω) 2 → L 2 (Ω) 2d with H -1 (Ω) the dual space to H 1 0 (Ω). Formal manipulations yield * δ h = - ∇ • h 1 ∇ • h 2 + M δ T h 1 h 2 , ∀ h = (h 1 , h 2 ) ∈ L 2 (Ω) 2d .
One can check by application of Green's formula that for all (α o δ , β o δ ) and (α, β) in

H 1 0 (Ω) 2 * δ δ (α o δ , β o δ ), (α, β) L 2 (Ω) = δ (α o δ , β o δ ), δ (α, β) L 2 (Ω) . (19)
Remark 2. The knowledge of the boundary values of α and β over the entire boundary ∂Ω represents redundant data. However, such assumption simplifies the following analysis since only essential boundary conditions, i.e. Dirichlet boundary conditions, can be employed in the variational formulation of the problem. When (α, β) is known only in a subdomain ∂Ω 0 ⊆ ∂Ω then the gradient equation ( 10) can be used to deduce natural boundary conditions, in terms of fluxes ∇α • n and ∇β • n on ∂Ω \ ∂Ω 0 , which would then appear in the integration by parts [START_REF] Knowles | On the representation of the elasticity tensor for isotropic media[END_REF].

Finally, note that an alternative to solving system [START_REF] Marsden | Mathematical foundations of elasticity[END_REF] is the Adjoint-Weighted variational Equation method presented in [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF]. It is based on a variational formulation that features a weighting operator as a substitute for the adjoint operator * δ in (18).

Existence and Uniqueness

In what follows we study existence, uniqueness and continuity estimates for the solutions to Eqn. [START_REF] Marsden | Mathematical foundations of elasticity[END_REF]. Although this analysis is fairly standard when the noise parameter δ is fixed, it remains non trivial to obtain continuity estimates that are independent of δ. This is the focus of the following developments. To this aim, let us express the following hypotheses that describe the assumed uniform boundedness w.r.t. δ of the noisy coefficients (a δ , b δ , c δ , d δ ) featured in Sec. 5.1.

Assumptions 3. There exists B > 0 such that for all δ > 0 one has

a δ , b δ , c δ , d δ ∈ L ∞ B (Ω) d := h ∈ L ∞ (Ω) d , h L ∞ (Ω) ≤ B .
We first prove a uniqueness results associated with the operator * δ δ .

Lemma 2. If Assumptions 3 hold, then for all δ > 0 the operator * δ δ : 2 then multiplying by (α, β), integrating over Ω and integrating by parts yield

H 1 0 (Ω) 2 → H -1 (Ω) 2 is injective. Proof. If * δ δ (α, β) = 0 for some (α, β) ∈ H 1 0 (Ω)
δ (α, β) 2 L 2 (Ω) = Ω |∇α + a δ α + b δ β| 2 + |∇ β + c δ α + d δ β| 2 dx = 0, (20)
thus we have δ (α, β) = 0 throughout Ω, which implies

|∇α| + |∇ β| ≤ C |α| + | β| in Ω,
where C is a constant that depends on B. Then either an approach similar in spirit to establishing uniqueness in Theorem 1, or the unique continuation principle, see for instance [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]Lemma 8.5] adapted to gradient equations and H 1 0 (Ω) solutions, allows to conclude that (α, β) ≡ (0, 0) throughout Ω.

Next, we show an inequality of Poincaré-Friedrichs type which is of key importance for the ensuing developments and numerical schemes.

Proposition 2. Provided that Assumptions 3 are satisfied, then there exits C > 0 depending on B but not on δ such that for all δ > 0 one has

α 2 H 1 0 (Ω) + β 2 H 1 0 (Ω) ≤ C δ (α, β) 2 L 2 (Ω) , ∀(α, β) ∈ H 1 0 (Ω) 2 . ( 21
)
Proof. By contradiction, assume that, for any δ ≥ 0 the inequality [START_REF] Fung | Foundations of solid mechanics[END_REF] does not hold, i.e. one can construct a sequence

(α δ(n) , βδ(n) , a δ(n) , b δ(n) , c δ(n) , d δ(n) ) n≥1 ∈ H 1 0 (Ω) 2 × L ∞ B (Ω) d 4 , such that αδ(n) 2 H 1 0 (Ω) + βδ(n) 2 H 1 0 (Ω) > C n δ(n) (α δ(n) , βδ(n) ) 2 L 2 (Ω) .
After re-normalization and without loss of generality we can assume that

(i) αδ(n) 2 H 1 0 (Ω) + βδ(n) 2 H 1 0 (Ω) = 1, (ii) δ(n) αδ(n) , βδ(n) 2 L 2 (Ω) < 1 n .
From (i), the sequence (α δ(n) , βδ(n) ) is bounded in the Hilbert space H 1 0 (Ω) 2 and by definition of the Banach space L ∞ B (Ω), there exist

α * , β * in H 1 0 (Ω) and a * , b * , c * , d * ∈ L ∞ B (Ω) d such that up to subsequences one has (α δ(n) , βδ(n) ) H 1 0 (α * , β * ) and (α δ(n) , βδ(n) ) → L 2 (α * , β * ), (a δ(n) , b δ(n) , c δ(n) , d δ(n) ) * L ∞ (a * , b * , c * , d * ). (22)
Next, one seeks to assess the weak convergence limit of δ(n) αδ(n) , βδ(n) , i.e. whether

δ(n) αδ(n) , βδ(n) , (v, w) L 2 (Ω) -→ n→+∞ * α * , β * , (v, w) L 2 (Ω) , ∀(v, w) ∈ L 2 (Ω) 2d .
where the operator * is defined as in [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF] in terms of a * , b * , c * and d * . Given that

δ(n) αδ(n) , βδ(n) - * α * , β * , (v, w) L 2 (Ω) = + Ω ∇(α δ(n) -α * ) • v + ∇( βδ(n) -β * ) • w dx + Ω (a δ(n) -a * ) • v + (c δ(n) -c * ) • w α * + (b δ(n) -b * ) • v + (d δ(n) -d * ) • w β * dx + Ω (a δ(n) • v + c δ(n) • w)(α δ(n) -α * ) + (b δ(n) • v + d δ(n) • w)( βδ(n) -β * ) dx.
Owing to the respective convergences [START_REF] Bal | Inverse anisotropic conductivity from internal current densities[END_REF] of the right-hand side terms in the previous equation, i.e. weak convergence in

H 1 0 , weak- * in L ∞ and strong in L 2 , one can conclude that δ(n) αδ(n) , βδ(n) converges weakly in L 2 to * α * , β * . This implies * α * , β * L 2 (Ω) ≤ lim inf n→+∞ δ(n) αδ(n) , βδ(n) L 2 (Ω) = 0.
Finally, one has * α * , β * = 0 which using Lemma 2 yields α * = β * = 0. Remark that

1 n ≥ δ(n) αδ(n) , βδ(n) 2 L 2 (Ω) ≥ 1 2 αδ(n) 2 H 1 0 (Ω) + βδ(n) 2 H 1 0 (Ω) - a δ(n) αδ(n) + b δ(n) βδ(n) 2 L 2 (Ω) + c δ(n) αδ(n) + d δ(n) βδ(n) 2 L 2 (Ω) ,
and so, using the strong convergence in L 2 :

4 B 2 αδ(n) 2 L 2 (Ω) + βδ(n) 2 L 2 (Ω) ≥ 1 2 - 1 n ⇒ α * 2 L 2 (Ω) + β * 2 L 2 (Ω) ≥ 1 4 B 2 , which contradicts α * = β * = 0.
With Proposition 2 at hand, we are now in position to prove existence and uniqueness of a solution to Eqn. [START_REF] Marsden | Mathematical foundations of elasticity[END_REF]. Given data satisfying Assumptions 3, finding (α o δ , β o δ ) ∈ H 1 0 (Ω) 2 solving ( 18) is equivalent to solving the weak form

δ (α o δ , β o δ ), δ (α, β) L 2 (Ω) = f δ -δ (α , β ), δ (α, β) L 2 (Ω) , ∀(α, β) ∈ H 1 0 (Ω) 2 . ( 23
)
From the Cauchy-Schwarz inequality, [START_REF] Mehrabadi | Strain energy density bounds for linear anisotropic elastic materials[END_REF] and the definition of the space L ∞ B (Ω) then there is a constant C > 0 such that

δ (α o δ , β o δ ), δ (α, β) L 2 (Ω) ≤ C α o δ 2 H 1 0 (Ω) + β o δ 2 H 1 0 (Ω) 1 2 α 2 H 1 0 (Ω) + β 2 H 1 0 (Ω) 1 2 . ( 24
)
Moreover, Proposition 2 entails

1 C α o δ 2 H 1 0 (Ω) + β o δ 2 H 1 0 (Ω) ≤ δ (α o δ , β o δ ) 2 L 2 (Ω) . (25)
With the boundedness [START_REF] Bal | Boundary control of elliptic solutions to enforce local constraints[END_REF] and coercivity (25) now verified, the existence of a unique solution to the variational problem [START_REF] Monard | Inverse anisotropic conductivity from power densities in dimension n ≥ 3[END_REF] follows directly from the Lax-Milgram theorem which ensures that

α o δ 2 H 1 0 (Ω) + β o δ 2 H 1 0 (Ω) 1 2 ≤ C f δ -δ (α , β ) L 2 (Ω) . 6
(α o , β o ) = f -(α , β ) and δ (α o δ , β o δ ) = f δ -δ (α , β ) + f δ
where the subscript is omitted in noise-free quantities (α, β) ≡ (α 0 , β 0 ), ≡ 0 and f ≡ f 0 , see Section 5.1. Accounting for the fact that f δ no longer belongs to the range of the operator δ , the additional term f δ satisfies * δ f δ = 0. Defining the operator

E δ (α, β) := (M -M δ ) α β , then the previous equations entail * δ δ (α o δ -α o , β o δ -β o ) = * δ f δ -f + E δ (α, β)
where α o δα o = β o δβ o = 0 on the boundary ∂Ω. Therefore, upon taking the L 2 (Ω)inner product in the above equation with the vector (α o δα o , β o δβ o ), one obtains

δ (α o δ -α o , β o δ -β o ) 2 L 2 (Ω) = f δ -f + E δ (α, β), δ (α o δ -α o , β o δ -β o ) L 2 (Ω) , so that δ (α o δ -α o , β o δ -β o ) L 2 (Ω) ≤ f δ -f + E δ (α, β) L 2 (Ω)
. Assumptions 1 together with Proposition 2 yield the fundamental result of this section: Proposition 3. For all δ ≥ 0, if Assumptions 3 hold and f δ ∈ L 2 (Ω) 2d then there exists C > 0 which does not depend on δ such that

α o δ -α o H 1 0 (Ω) + β o δ -β o H 1 0 (Ω) ≤ C f δ -f L 2 (Ω) + E δ . where E δ 2 := v=a,b,c,d v -v δ 2 L 2 (Ω) .

Solution stability w.r.t. noisy displacement measurements

We now clarify the meaning of ( 14) and deduce a stability result for the reconstructed parameters in terms of the approximation quality associated with the noisy data.

Theorem 2. If assumptions 2 and 3 hold and there exists an increasing function η(δ) such that, for all δ > 0, the functions (u 1,δ , u 2,δ ) and

(ε 1,δ , ε 2,δ , H 1,δ , H 2,δ ) satisfy 2 n=1 u n,δ -u n L ∞ (Ω) + ε n,δ -ε n L ∞ (Ω) + d k=1 H k n,δ -∇ ⊗2 (u n •e k ) L 2 (Ω) ≤ η(δ), (26)
then, for δ small enough, the solution to [START_REF] Marsden | Mathematical foundations of elasticity[END_REF] with ( δ , f δ ) defined in Section 5.1 satisfy

α o δ -α o H 1 0 (Ω) + β o δ -β o H 1 0 (Ω) ≤ C η(δ)
, where C is a positive scalar which does not depend on δ.

Proof. The proof uses standard estimates as well as Assumptions 2. A direct application of Proposition 3 is used to conclude. Remark 3. As the proof of Theorem 2 relies on bounding the term E δ , it requires to estimate either the strain tensor ε n,δ in the L ∞ -norm and the hessian tensor H k n,δ in the L 2 -norm, or conversely. The choice adopted in ( 26) is actually the more practical. Remark 4. Theorem 2 is consistent with the derivation of Eqn. [START_REF] Marsden | Mathematical foundations of elasticity[END_REF] in that no error has been considered in the knowledge of the boundary values of α and β on ∂Ω.

Numerical differentiation scheme

We now address the question of constructing explicitly the quantities {ε 1,δ , ε 2,δ } and {H 1,δ , H 2,δ } from noisy measurements u 1,δ , u 2,δ . In [START_REF] Kohn | A variational method for parameter identification[END_REF] the discrepancy between noisy and exact measurements is assumed to be small in H 1 -norm. Yet, in practice the functions u n,δ may not have bounded gradients or second-order derivatives. Therefore, we aim at defining a differentiation operator

D : (u 1,δ , u 2,δ ) ∈ L ∞ (Ω) d 2 → (ε 1,δ , ε 2,δ , H 1,δ , H 2,δ ) ∈ L ∞ (Ω) d×d 2 × L 2 (Ω) d×d×d 2 ,
such that, for δ small enough, Eqn. ( 26) holds with a function η(δ) to be determined. Therefore, we are required to regularize the measurements to define the operator D. While various methods have been discussed in the literature, see e.g. [START_REF] Engl | Regularization of Inverse Problems[END_REF], we discuss hereinafter a regularization by L 2 -projections onto coarse finite element spaces and show in which sense this approach is optimal.

Approach overview

Consider a partition of the domain Ω ⊂ R d using a number of finite elements K e that are characterized by a mesh size h with 0 < h < 1 and such that

with C independent of h. This partition is associated with the finite element spaces

L r h = {v h ∈ L 2 (Ω) | (v h )| Ke ∈ Q r e ∀ K e ⊂ Ω }
with Q r e a local finite element space on K e that includes at least the polynomials of order r. Let V h ⊂ H 1 0 denote a finite element space approaching H 1 0 (Ω) in the limit h → 0.

The approach considered is as follows: Noisy displacement fields measurements u n,δ are assumed to be available on a fine grid. They are interpreted as functions whose components belong to a low-order finite element space, such as L 0 h 0 , on a fine mesh. Note that this choice is arbitrary and L 1 h 0 or a standard conforming finite element space can be used for interpolating displacement data. Next, we define regularization and differentiation schemes to construct the strain and hessian tensors ε n,δ , H n,δ of these quantities. They are associated to coarser meshes but higher order finite elements, i.e. to functional spaces L r h 1 and L r h 2 respectively, with h 0

h 1 ≤ h 2 . Reconstruction of moduli (α o δ , β o δ
) is finally achieved upon solving [START_REF] Marsden | Mathematical foundations of elasticity[END_REF] in V h 0 on the fine discretization.

Table 1: Summary of the approach where r ≥ 2 and

h 2 ≥ h 1 h 0 . Data Gradient Hessian Reconstruction Quantities u n,δ ε n,δ H n,δ (α o δ , β o δ ) Fct. spaces L 0 h 0 L r h 1 L r h 2 V h 0
Note that in the ensuing analysis, we do not account for the numerical errors associated with (i) the numerical quadratures required for the computation of ε n,δ and H n,δ , and (ii) the computation of (α o δ , β o δ ) using the discretized version of the variational formulation discussed above.

The L 2 -projection operator

Let P r h denote the scalar orthogonal projection operator on the finite element space of order r considered

P r h : u ∈ L 2 (Ω) -→ u h,r ∈ L r h
where u h,r is defined as the unique function of L r h such that

(u h,r , v h,r ) L 2 (Ω) = (u, v h,r ) L 2 (Ω) ∀ v h,r ∈ L r h .
As the spaces L r h are not conforming this equation is equivalent to have, for all K e ⊂ Ω,

(u h,r , v h,r ) L 2 (Ke) = (u, v h,r ) L 2 (Ke) ∀ v h,r ∈ Q r e , (28)
which defines the local orthogonal projection operator P r h,e that is associated with element K e and satisfies P r h,e (u| Ke ) = u h,r | Ke . Standard projection properties are satisfied locally and globally, i.e.

P r h,e (u h,r | Ke ) = u h,r | Ke and P r h (u h,r ) = u h,r ∀ u h,r ∈ L r h P r h,e (u| Ke ) L 2 (Ke) ≤ C u L 2 (Ke) and P r h (u) L 2 (Ω) ≤ C u L 2 (Ω) ∀u ∈ L 2 (Ω).
We now state an approximation result, which can be proven using standard results of the theory of finite element method [START_REF] Brenner | The mathematical theory of finite element methods[END_REF][START_REF] Heimsund | Superconvergence for the gradient of finite element approximations by L2 projections[END_REF]. This result quantifies how the local L 2projection onto high-order finite element spaces yields a smooth local approximation of a given function u from a L 2 approximation of this function. In these inequalities, r denote the order of the featured finite element space while the positive integers and m stand for the Sobolev indices respectively associated with the a priori regularity of the function u considered and the required regularity of the approximation.

Proposition 4. For any u ∈ H (K e ), v ∈ L 2 (K e ) and element K e ⊂ Ω the following inequalities hold -If m < ≤ r + 1, we have

u -P r h,e (v) H m (Ke) ≤ C h -m u -v L 2 (Ke) + C h -m |u| H (Ke) . (29)
-If m + d/2 < ≤ r + 1 then u -P r h,e (v) W m,∞ (Ke) ≤ C h -m-d/2 u -v L 2 (Ke) + C h -m-d/2 |u| H (Ke) . (30)
The interpretation of Proposition 4 is as follows: Within the framework of finite element spaces, one can approximate locally in element K e , the derivatives up to order m of a given function u ∈ H (K e ) using a function v ∈ L 2 (K e ) that approximates u in the L 2norm only. Yet, the approximation quality is penalized by the detrimental term h -m , with h being the mesh size. Therefore, to obtain an optimal approximation in the sense of Theorem 2, the parameter h has to be chosen in such a way that, roughly speaking, δh -m 1. This issue is discussed in the next section.

Construction of noisy operator δ by regularization

Consider noisy displacement measurements u 1,δ , u 2,δ ∈ L ∞ (Ω) d satisfying for all δ > 0

u 1 -u 1,δ L ∞ (Ω) + u 2 -u 2,δ L ∞ (Ω) ≤ C δ. (31)
The differentiation operator D is now given by

D(u 1,δ , u 2,δ ) = (ε 1,δ , ε 2,δ , H 1,δ , H 2,δ ) (32)
where the strain and hessian tensors are defined element-wise in all elements K e ⊂ Ω by

ε n,δ | Ke = ∇ s d k=1 P r h 1 ,e (u n,δ | Ke • e k )e k , H n,δ | Ke = d k=1 ∇ ⊗2 P r h 2 ,e (u n,δ | Ke • e k ) ⊗ e k .
For the sake of generality, two different meshes of sizes h 1 and h 2 are used to project the data (u 1,δ , u 2,δ ) to compute first and second-order derivatives respectively. For simplicity, we assume that the function η(δ) in ( 26) is given by a power law. Now, the aim is to tune the parameters h 1 and h 2 as functions of δ to obtain the best possible estimate in [START_REF] Lin | Conditional stability for the hexagonal anisotropic elastic dynamical systems[END_REF]. For the sake of argument, let ≥ 0 ≥ 2 such that u n , u n,δ ∈ H (Ω) ∩ W -1,∞ (Ω) for all δ > 0 and n = 1, 2. In this case one can show using classic interpolation results, see [START_REF] Adams | Sobolev spaces[END_REF]Theorem 4.17], that a direct differentiation of the measurements u 1,δ , u 2,δ yields the convergence rate η(δ) = δ -2 . Therefore this is the optimal bound that constrains the final estimate on the reconstruction.

Theorem 3. Let Assumptions 2 be satisfied and assume that the solutions to (5) satisfy u 1 , u 2 ∈ H (Ω) d where is an integer such that 3 ≤ ≤ r + 1, while the noisy measurements satisfy [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]. For δ sufficiently small, when choosing

h 1 = δ 1 -d/2
and h 2 = δ 1 and constructing the operator δ using D defined by [START_REF] Heimsund | Superconvergence for the gradient of finite element approximations by L2 projections[END_REF], then the reconstructed parameters (α o δ , β o δ ) obtained solving [START_REF] Marsden | Mathematical foundations of elasticity[END_REF] satisfy

α o δ -α o H 1 0 (Ω) + β o δ -β o H 1 0 (Ω) ≤ C δ -2 , ( 33
)
where the constant C > 0 does not depend on δ.

Proof. From Definition [START_REF] Albocher | Adjoint-weighted equation for inverse problems of incompressible plane-stress elasticity[END_REF] and Property [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF] we have

H k n,δ -∇ ⊗2 (u n • e k ) L 2 (Ke) ≤ C P r h 2 ,e (u n,δ ) -u n H 2 (Ke) ≤ C h 2 -2 u n -u n,δ L 2 (Ke) + C h 2 -2 |u n | H (Ke) .
Squaring the above, summing over elements K e and taking the square root, entail

H k n,δ -∇ ⊗2 (u n • e k ) L 2 (Ω) ≤ C h 2 -2 u n -u n,δ L 2 (Ω) + C h 2 -2 Ke |u n | 2 H (Ke) 1/2 .
As u 1 , u 2 ∈ H (Ω) d , using Eqn. [START_REF] Brenner | The mathematical theory of finite element methods[END_REF] and the boundedness of Ω finally yield

H k n,δ -∇ ⊗2 (u n • e k ) L 2 (Ω) ≤ C δ h 2 2 + h 2 -2 . (34)
Moreover, owing to Property [START_REF] Kohn | A variational method for parameter identification[END_REF], one has

ε n,δ -ε n L ∞ (Ω) ≤ sup Ke P r h 1 (u n,δ ) -u n W 1,∞ (Ke) ≤ C sup Ke h -1-d/2 1 u n -u n,δ L 2 (Ke) + h -1-d/2 1 |u n | H (Ke) .
Using Eqn. [START_REF] Brenner | The mathematical theory of finite element methods[END_REF], the first term can be bounded uniformly with respect to K e ⊂ Ω as

h -1-d/2 1 u n -u n,δ L 2 (Ke) ≤ Ch -1 1 u n -u n,δ L ∞ (Ke) ≤ C h -1
1 δ, since the mesh is assumed to be uniform for all h 1 owing to [START_REF] Lin | Quantitative strong unique continuation for the Lamé system with less regular coefficients[END_REF]. The above result and the inequality sup

Ke |u n | H (Ke) ≤ C|u n | H (Ω) yield ε n,δ -ε n L ∞ (Ω) ≤ C δ h 1 + h -1-d/2 1 . (35)
Combining [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF] and (35) finally entails

2 n=1 ε n,δ -ε n L ∞ (Ω) + d k=1 H k n,δ -∇ ⊗2 (u n •e k ) L 2 (Ω) ≤ C δ h 1 + h -1-d/2 1 + δ h 2 2 + h -2 2 .
To conclude the proof, let h 1 = δ r 1 and h 2 = δ r 2 with the positive real parameters r 1 and r 2 chosen to maximize the quantities min(r 1 ( -1-d/2), 1-r 1 ) and min(r 2 ( -2), 1-2r 2 ).

Since ≥ 3, one can show that these maxima are achieved when r 1 = 1/(d/2) and r 2 = 1/ . Then Theorem 2 is used to finally obtain Estimate [START_REF] Adams | Sobolev spaces[END_REF].

Remark 5. When d = 2, the choice h 1 = h 2 = δ 1/ is actually sufficient to obtain [START_REF] Adams | Sobolev spaces[END_REF]. However, in this case the featured constant C would be suboptimal compared to this obtained using Theorem 3.

Numerical results

A set of numerical results is presented in this section to assess the performances of the proposed approach. These numerical examples correspond to solving Eqn. [START_REF] Marsden | Mathematical foundations of elasticity[END_REF] in dimension d = 2 with the noisy operators constructed using the method presented in Section 7.3 with h 1 = h 2 ≡ h. Given a distribution of smooth enough constitutive parameters (α, β), the solutions (u 1 , u 2 ) are computed on a sufficiently fine mesh from Eqn. ( 1) augmented with non-homogeneous Dirichlet boundary conditions

u n = g n on ∂Ω, n = 1, 2. ( 36
)
These solutions are then polluted by a noise, which, for the sake of reproducibility, is parameterized and therefore deterministic, as an increasing function in L ∞ (Ω)-norm of a parameter δ. In what follows, we assess the behavior of the proposed numerical algorithm with respect to the parameters h and δ. Let Ω = [0, 1] 2 and the parameters α(x) and β(x) by defined as

α(x) = α 0 + Nα i=1 α i c(|x -x α i | ; a - i , a + i ), β(x) = β 0 + N β i=1 β i c(|x -x β i | ; b - i , b + i ), ( 37
)
where the function c(r ; r -, r + ) equals 1 if r < r -, 0 if r > r + and c(r ; r -, r

+ ) = 1 -(r -r -)/(r + -r -) 2 1 + 2(r -r -)/(r + -r -) r -≤ r ≤ r + .
The radial function c(|x| ; r -, r + ) is C 1 (Ω) while all its second order derivatives are L ∞ (Ω). This property guarantees that the solutions u 1 , u n to (1) and (36) are at least in H 3 (Ω) for smooth enough boundary conditions g n . In turn, it ensures that Assumptions 2 are satisfied owing to the injection of H 3 (Ω) into W 1,∞ (Ω). Therefore, the algorithm can at least be defined at the continuous level and Theorem 3 grants us theoretical convergence of the algorithm as δ → 0. Numerical computations are performed using continuous fifth-order nodal finite elements on a square mesh based on Gauss-Lobatto points as described in [START_REF] Cohen | Higher-order numerical methods for transient wave equations[END_REF]. We use fifth-order quadrature formulae based on Gauss points for the local projection [START_REF] Bochev | Least-Squares Finite Element Methods[END_REF] and Gauss-Lobatto points for the computation of the finite element matrices. The synthetic measurements u 1 and u 2 are computed using a conjugate gradient (CG) technique, on a reference fine grid of size h 0 = 1/120, i.e. L 0 h 0 = V h 0 and dim(V h 0 ) = 601 × 601. Then, deterministic noise is added to the solution as

u n,δ (x) = u n (x) + δ M m=-M |m| M χ |m| M x √ δ ,
with χ(x) = cos(2π x • e 1 ) cos(2π x • e 2 )(e 1 + e 2 ) and M = 20. At the discrete level, the noise is directly interpolated on the nodal functions generating L 0 h 0 . Moreover, Assumption [START_REF] Brenner | The mathematical theory of finite element methods[END_REF] is satisfied by construction, yet u n,δ (x)u n (x) H 2 (Ω) = O(1), which entails that a non-regularized direct data differentiation approach with dim V h 0 = +∞ should not converge as δ → 0.

Finally, it is considered that the solution to the elliptic problem ( 18) is computed sufficiently accurately with a CG method, i.e. the space V h 0 in Table 1 is large enough and the relative stopping criterion of the CG method is good enough.

Simple static case. The parameters α and β are given by formulae (37) with

α 0 = 22, β 0 = 2, N α = N β = 1, α 1 = β 1 = 18, a - 1 = b - 1 = 0.1, a + 1 = b + 1 = 0.2, x α 1 = x β 1 = [1/2, 1/2]
T , see Fig. 3a. The static case corresponds to ω 1 = ω 2 = 0 in (5). As discussed in Section 3.2, the invertibility condition in Assumptions 2 can be satisfied by choosing appropriate Dirichlet boundary conditions g 1 and g 2 in (36), such as Computed elasticity solutions u 1 and u 2 are represented Figure 1, together with the term det , defined by ( 9), which does not vanish so that Assumptions 2 are satisfied.

g 1 (x) = 1 + (x • e 2 ) e 1 + (x • e 1 ) e 2 , g 2 (x) = 1 + (x • e 1 ) e 1 + (x • e 2 ) e 2 . ( 38
) (a) u 1 (x) • e 1 (b) u 1 (x) • e 2 (c) u 2 (x) • e 1 (d) u 2 (x) • e 2
Convergence in h. The reconstruction behavior is investigated by solving Eqn. ( 18) for different mesh sizes h in the definition of δ from Section 7.3. This parameter is varied as h = k h 0 using an integer k and given the size h 0 of the reference mesh associated with the noisy data. Figure 2 shows the relative H 1 (Ω)-error of the reconstruction for noise values δ = 10 -5 and δ = 10 -7 . In the light of Eqns. ( 34) and (35), Figure 3 shows that the computation of the data derivatives is penalizing if the mesh employed is too fine, as negative powers of h are involved. Alternatively, if the mesh is too coarse then the quality of the solutions approximation is too deteriorated, as they involves positive powers of h. Hence, there exists an optimal value h corresponding, in each case, to the minimum of the associated curve in Figure 2.

Convergence in δ. The convergence of the reconstruction error with respect to the noise level δ is now investigated. In Figure 4, the relative reconstruction errors in H 1 (Ω)-norm are compared for different values h. Two expected trends are highlighted. Firstly, at a given noise value δ there exists an optimal mesh size parameter h, for the projection and differentiation steps, for which the reconstruction is the best in H 1 (Ω)- Secondly, for a fixed h, when the noise level decreases in L (Ω)-norm the reconstruction quality reaches a plateau. This is due to a loss of resolution associated with the projection on a coarse mesh. These numerical results essentially show that the proposed algorithm achieves the standard trade-off between regularization and resolution.

Frequency dependent case. In light of formulae ( 4) and ( 9), it is clear that using oscillating elasticity solutions, with locally vanishing gradients, might be detrimental to the invertibility condition of Assumptions 2. This implies that low frequency configurations should be preferred so that we set ω 1 = 1 and ω 2 = 0. The boundary conditions are prescribed as

g 1 (x) = 1 + (x • e 1 ) e 1 + (x • e 2 ) e 2 , g 2 (x) = (1 + x • e 1 + x • e 2 )(e 1 -e 2 ).
The exact moduli distributions are computed using (37) and displayed in Fig. 6a. The modulus of the corresponding elasticity solutions and det are plotted in Figure 5.

Figure 6 shows three reconstructions. Fig. 6b corresponds to a noise-free configuration (δ = 0) and no regularization (h = h 0 = 1/120). The associated relative reconstruction error in H 1 (Ω)-norm is 0.0033. Figure 6c corresponds to δ = 10 -7 and mesh size h = 1/120, i.e. a reconstruction without regularization. The corresponding relative error is 0.83 and the reconstruction of α(x), β(x) is poor. Finally, for this noise value but with regularization, i.e. h = 1/24, then the relative error decreases to 0.67 and the reconstruction is qualitatively improved in terms of identification of the number of heterogeneities, their locations and relative sizes, see Fig. 6d.

Random parameters distributions This section is concluded with the investigation of a case of randomly varying moduli. Their studied exact distributions in Fig. 7a are constructed from (37), using N α = N β = 1000 and random parameters. They are characterized by a correlation length much smaller than the domain size in order to give rise to oscillating strain field solutions. The settings are these of the previous paragraph but ω 1 = ω 2 = 0. The reconstruction in Fig. 7b from noise-free measurements is fairly good, with a relative error of 0.0067. A noise value δ = 10 -6 makes the non-regularized reconstruction ineffective: Fig. 7c, relative error 0.87. Alternatively, the recovery of the unknowns is remarkably improved using the proposed regularization method as shown in Fig. 7d with a relative reconstruction error of 0.20.

Conclusion

An algorithm for the quantitative reconstruction of constitutive elasticity parameters has been investigated. It is based on the construction and inversion of a linear operator computed from full-field internal measurements, possibly noisy, of elasticity solutions. The highlights of this study are: (i) Construction of a linear operator from available full-field measurements which in turn are required to satisfy some invertibility and compatibility conditions. (ii) Characterization of this symmetric definite positive linear operator and derivation of an associated variational formulation. (iii) Investigation of a numerical differentiation scheme based on regularizations by projections on coarse meshes but within a high-order finite element formulation. (iv) Theoretical results proving convergence of the algorithm when noise in the data vanishes in the L ∞ (Ω)-norm yet not necessarily in the H 2 (Ω)-norm. (v) The proposed algorithm is theoretically valid at any frequency. Invertibility and compatibility conditions make it applicable at least to low-frequency configurations. (vi) Numerical results highlighting the method's potential for practical applications.

Looking forward, it would be relevant to address the cases where the invertibility or compatibility conditions might locally be not satisfied, e.g. when dealing with highly oscillating measurements. A possible strategy is to work with a larger set of data so that inversion formula [START_REF] Barbone | Adjoint-weighted variational formulation for the direct solution of inverse problems of general linear elasticity with full interior data[END_REF] can be constructed pointwise using the pair of solutions maximizing | det | locally. Moreover, the proposed approach might be extended to the case of non-complete data, where either internal measurements or boundary conditions are partially available. Finally, while the approach finds direct applications in the field of non-destructive material testing, its extension to the model of nearly incompressible solids used in medical imaging is the subject of ongoing research.

K

  e with K e ∩ K e = ∅ if e = e and Ke 1 dx ≤ Ch d ,

Figure 1 :

 1 Figure 1: Computed elasticity solutions in the static case.

  Fig.3b

Fig. 3c

  Fig. 3c

Fig. 3d

  Fig. 3d

Figure 2 :

 2 Figure 2: Relative error in H 1 (Ω)-norm on reconstructed parameters (α, β) w.r.t. mesh size h = k h 0 . The doubly-circled dots correspond to the reconstructions in Figure 3.

Figure 3 :

 3 Figure 3: Reconstruction β δ for different values of h = k h 0 with δ = 10 -7 .

  Fig. 3c

Fig. 3d

  Fig. 3d

Fig. 3b

  Fig.3b

Figure 4 :

 4 Figure 4: Relative error in H 1 (Ω)-norm on parameters (α, β) w.r.t. noise value δ for different values of h = k h 0 . The doubly-circled dots correspond to Figure 3.

Figure 5 :

 5 Figure 5: Computed elasticity solutions in the frequency dependent case

Figure 6 :

 6 Figure 6: (a) Exact values of α(x) (top) and β(x) (bottom); Corresponding reconstructions with (b) no noise nor regularization, (c) with noise but no regularization, (d) with noise and regularization.

Figure 7 :

 7 Figure 7: (a) Exact values of α(x) (top) and β(x) (bottom); Corresponding reconstructions with (b) no noise nor regularization, (c) with noise but no regularization, (d) with noise and regularization.

  . Stability of reconstruction formula with noisy data 6.1. Solution stability w.r.t. noisy operator δ

	Looking forward to the proposition of a regularization scheme in Section 7 we now derive
	a stability result in the presence of noise. Considering Assumptions 1 and Eqn. (17), let
	(α o , β o ) and (α o δ , β o δ ) denote the solutions in H 1 0 (Ω) 2 to Eqn. (18) that are respectively
	associated with exact measurements and noisy data satisfying Assumptions 3. Then
	one has


Fatima Zahra

Frédéric Jamin ¹ , Lkhadir Lakhal ² El Korchi¹²

Mohamed Elomari²
email: elomari@uca.ma

Moulay Said

El Youssoufi

Effect of water on granular matter mechanics, macroscopic scale: Collapse of granular media during wetting

Keywords: Collapse, wetting, Granular media, Glass bead

In soil mechanics, it is well known that the behaviour of the granular soils, as the sand, in the dry state is similar to the saturated state. Between these two extremes states exist various domains of water content often named: hygroscopic, pendular and funicular.

In the funicular domain, the partially saturated soils can undergo a collapse due to the only

variations of water content. This phenomenon was studied under a wetting paths loading by [START_REF] Barden | The collapse mechanism in partly saturated soil[END_REF]. In the literature, several models describe the collapse behaviour as [START_REF] Alonso | A constitutive model for partially saturated soils[END_REF] with a phenomenological model based on hydro-mechanical approach, or as [START_REF] Liu | Numerical study of soil collapse behaviour by discrete element modeling[END_REF] with a numerical model based on the micromechanical approach by Discrete Elements Method.

In the pendular domain, the granular media are more sensitive to the water content, which can be at the origin of an apparent cohesion by capillarity; this behaviour has been shown by simple compression tests [START_REF] Soulié | Capillary cohesion and mechanical strength of polydisperse granular materials[END_REF] and by direct shear tests [START_REF] Richefeu | Shear strength properties of wet granular materials[END_REF].

The aims of this experimental study is to identify and analyse the collapse phenomenon on a granular material constituted by a polydisperse assembly of glass bead during a wetting path from pendular domain to saturation with two different grain sizes. All samples are prepared with initial gravimetric water content w 0 = 5%. This initial water content corresponds to pendular domain and thus creates capillary bridges between grains witch cause apparent cohesion.

After the positioning of the sample in the triaxial cell, we arrange the axial and radial local sensors

ALERT Workshop -Poster session

Aussois 2014 on the sample. Local sensors are used for monitoring the evolution of axial and radial displacements at different points of the sample during the wetting path, without external stress. The wetting path is done by applying a gradient of water pressure ( u w = 5 kPa at the bottom of the sample by a pressure-volume controller).

The experimental results during the wetting path show that the collapse of the granular assembly is purely due to hydraulic effects, caused by the progressive fusion of the capillary bridges. The scope of the collapse is even more significant for smaller grain size (Figure 1). This phenomenon can be explained by a large number of contacts (initial number of capillary bridges per volume). Our results show that the collapse depends on the density of the granular assembly and on the grain size.

Aussois 2014

Illustrations

Glass bead Φ < 50 μm Glass bead Φ 100 -200 μm