Resolution Guarantees for the Reconstruction of Inclusions in Linear Elasticity Based on Monotonicity Methods

We deal with the reconstruction of inclusions in elastic bodies based on monotonicity methods and construct conditions under which a resolution for a given partition can be achieved. These conditions take into account the background error as well as the measurement noise. As a main result, this shows us that the resolution guarantees depend heavily on the Lam\'e parameter $\mu$ and only marginally on $\lambda$.


Introduction
In this paper we deal with the detection and reconstruction of inclusions in elastic bodies based on monotonicity methods, where the main focus lies on the so-called 'resolution guarantees'. Our results are of special importance, when considering reconstructions based on real data.
Before we introduce the definition of the resolution guarantees we present the setting. Let We base our considerations on the work [17], where the resolution guarantees for the electrical impedance tomography (EIT) problem were considered. Definition 1. An inclusion detection method that yields a reconstruction D R to the true inclusion D is said to fulfill a resolution guarantee w.r.t. a partition (ω s ) s=1,...,N ⊆ Ω if (i) ω s ⊆ D implies ω s ⊆ D R for s ∈ {1, 2, . . . , N} (i.e., every element that is covered by the inclusion will be marked in the reconstruction), (ii) D = ∅ implies D R = ∅ (i.e., if there is no inclusion then no element will be marked in the reconstruction).

Remark 1.
Obviously, a reconstruction guarantee will not hold true for arbitrarily fine partitions. The achievable resolution will depend on the number of applied boundary forces, the inclusion contrast, the background error and the measurement noise.
We aim to construct conditions under which a resolution for a given partition can be achieved. These conditions take into account the background error as well as the measurement noise. Thus, we show that it is generally possible to give rigorous guarantees in linear elasticity.
Before we go into more detail about the resolution guarantees, we would like to give an insight into the theory and methods of inclusion detection considered so far.
The theory of the inverse problem of linear elasticity, i.e., uniqueness results and Lipschitz stability studies, etc were examined, e.g. in the following works: [24,26,30] deal with the 2D case and [19] with two and three dimensions. For uniqueness results in 3D we want to mention [9,31,32] as well as [2,3,26,29,31], where some boundary determination results were proved. In addition, results concerning the anisotropic case can be found, e.g. in [20,22,23]. Finally, [21] discussed the reconstruction of inclusion from boundary measurements.
The following methods, among others, have been used to solve the inverse problem of linear elasticity: A Landweber iteration method was applied in [18,28]. Further on, [25,27] considered regularization approaches. Beside the aforementioned methods, adjoint methods were used in [33][34][35]. Further on, [1,10,36] took a look at reciprocity principles. Finally, we want to mention the monotonicity methods for linear elasticity developed by the authors of this paper in [5].
We focus on the monotonicity methods which are built on the examinations in [37,38]. These methods were first used for EIT (see, e.g. [11-13, 15, 16]) and then on other problems such as elasticity (see, e.g. [5][6][7]). In short, for this method the monotonicity properties of the Neumann-to-Dirichlet operator play an essential role. All in all, this builds the basis for our considerations.
This paper is organized as follows. First, we introduce the problem statement for linear elasticity and the setting, where we distinguish between the continuous and discrete case. Next, we give a summary of the monotonicity methods, i.e., the standard monotonicity tests as well as the linearized monotonicity tests. In section 4, we present the background for the resolution guarantees and introduce the algorithms of the aforementioned monotonicity tests. Further on, we prove the required theorems which build the basis for the algorithms. Finally, we simulate the reconstruction for different settings. As a result, we conclude that the resolution guarantees depend heavily on the Lamé parameter µ and only marginally on λ.

Problem statement and setting
We take a look at the continuous and discrete setting and introduce the corresponding problems, the required assumptions concerning the inclusion, the background and the measurement error.

Continuous case
We start with the introduction of the problems of interest, e.g. the direct as well as the inverse problem of stationary linear elasticity.
For the following, we define Let u : Ω → R d be the displacement vector, µ, λ : Ω → L ∞ + (Ω) the Lamé parameters, ∇u = 1 2 ∇u + (∇u) T the symmetric gradient, where the rows of ∇u consist of the gradient of the components u i , i = 1, . . . , d, of u, n is the normal vector pointing outside of Ω, g ∈ L 2 (Γ N ) d the boundary force and I the d × d-identity matrix. We define the divergence of a Euclidean basis vector and x j a component of a vector from R d .
The boundary value problem of linear elasticity (direct problem) is that u ∈ H 1 (Ω) d solves From a physical point of view, this means that we deal with an elastic test body Ω which is fixed (zero displacement) at Γ D (Dirichlet condition) and apply a force g on Γ N (Neumann condition). This results in the displacement u, which is measured on the boundary Γ N . The equivalent weak formulation of the boundary value problem (1)- We want to remark that for λ, µ ∈ L ∞ + (Ω), the existence and uniqueness of a solution to the variational formulation (4) can be shown by the Lax-Milgram theorem (see e.g. in [4]).

Neumann-to-Dirichlet operator and its monotonicity properties.
Measuring boundary displacements that result from applying forces to Γ N can be modeled by the Neumann-to-Dirichlet operator Λ(λ, µ) defined by This operator is self-adjoint, compact and linear (see corollary 1.1 from [5]). Its associated bilinear form is given by where u g (λ,µ) solves the problem (1)-(3) and u h (λ,µ) the corresponding problem with boundary force h ∈ L 2 (Γ N ) d .

Discrete case
Next, we go over to the discrete case and take a look at the bounded domain Ω ⊂ R d with piecewise smooth boundary representing the elastic object. Further on, let λ, µ : Ω → R + be the real valued Lamé parameter distribution inside Ω.
We apply the forces g l on the Neumann boundary of the object, where the location of their support is denoted by Γ (l) N ⊆ Γ N , l = 1, . . . , M. We assume that the patches are disjoint. Thus, the discrete boundary value problem is given by The resulting displacement measurements are represented by the discrete version of Λ(λ, µ): with and u (k) solves the boundary value problem (7)-(10) for the boundary load g k .

Assumptions regarding the inclusion, the background as well as the measurement error
In the following, we introduce our assumptions and definitions concerning the Lamé parameters for the inclusion and background including their error considerations.
(a) Distribution of Lamé parameter (λ(x), µ(x)): where D denotes the unknown inclusion and B the background.
We distinguish between the following two cases where the lower bounds c λ and c µ are known.

Summary of the monotonicity methods
First, we state the monotonicity estimates for the Neumann-to-Dirichlet operator Λ(λ, µ) and denote by u g (λ,µ) the solution of problem (1)-(3) for the boundary load g and the Lamé parameters λ and µ.

be an applied boundary force, and let u
Then

be an applied boundary force, and let u
We give a short overview concerning the monotonicity methods, where we restrict ourselves to the case λ 1 ⩾ λ 0 , µ 1 ⩾ µ 0 . In the following, let D be the unknown inclusion and χ D the characteristic function w.r.t. D. In addition, we deal with 'noisy difference measurements', i.e., distance measurements between u g (λ,µ) and u g (λ0,µ0) affected by noise, which stem from the corresponding system (1)-(3).
We define the outer support in correspondence to [5] as follows: let ϕ = (ϕ 1 , ϕ 2 ) : Ω → R 2 be a measurable function, the outer support out ∂Ω supp(ϕ) is the complement (in Ω) of the union of those relatively open U ⊆ Ω that are connected to ∂Ω and for which ϕ| U = 0.

Standard monotonicity tests
We start our consideration with the standard monotonicity tests and take a look at the case for exact as well as noisy data. Here, we denote the material without inclusion by (λ 0 , µ 0 ) and the Lamé parameters of the inclusion by (λ 1 , µ 1 ).

Linearized monotonicity tests
We also introduce the linearized monotonicity tests as a modification of the standard methods. Similar as before, we deal with the exact as well as perturbed problem.

Tests for exact and noisy data.
Corollary 4. Linearized monotonicity test (Corollary 2.7 from [5]) Corollary 5. Linearized monotonicity test for noisy data (Corollary 2.9 from [5]) (17) with noise level δ > 0. Then for every open set ω ⊆ Ω there exists a noise level δ 0 > 0, such that for all 0 < δ < δ 0 , ω is correctly detected as inside or not inside the inclusion D by the following monotonicity test

Resolution guarantees
In this section we formulate the algorithms for the monotonicity tests, i.e., the standard monotonicity tests as well as the linearized tests and follow the considerations in [17], where resolution guarantees for EIT were analyzed.

Algorithms
Before we take a look at the algorithms for the reconstruction, we define the corresponding notations which we will use in the following. We set where the quantities are given in section 2.3 assumptions (a)-(d).

Algorithms for standard monotonicity tests.
We now formulate the algorithms for the standard monotonicity tests. We start with the case

Then the reconstruction D R is given by the union of the marked resolution elements.
Further on, we take a look at the case for 'smaller' Lamé parameter inclusions and assume such that

Algorithms for linearized monotonicity tests. Replacing the monotonicity test for the case
with their linearized approximations yields the linearized monotonicity test where κ λ , κ µ ∈ R are suitable contrast levels defined in the following algorithm. Further, we assume the λ max and µ max are global bounds with for all x ∈ Ω.

Algorithm 3. Mark each resolution element ω s for which
Then the reconstruction D R is given by the union of the marked resolution elements.
As for the standard monotonicity test, we formulate the linearized test for inclusions with smaller Lamé parameter which fulfill (λ Dmax , µ Dmax ) < (λ B min , µ B min ).

Algorithm 4. Mark each resolution element ω s for which
for s ∈ {1, 2, . . . , N}. Then the reconstruction D R is given by the union of the marked resolution elements.

Formulation of theorems
We will analyze the algorithms in more detail and take a look at the required theorems.

Theorem 1. The reconstruction of algorithm 1 fulfills the resolution guarantees if
where 'eig' stands for the set of eigenvalues of the input matrix.
The knowledge, that from (λ 1 , µ 1 ) ⩽ (λ 2 , µ 2 ) it follows that implies that Hence, so that ω s will be marked by the algorithm. This shows that part (i) of the resolution guarantee is satisfied. To prove part (ii) of the resolution guarantee, assume that D = ∅ and D R = ∅. Then there must be an index s ∈ {1, 2, . . . , N} with Again, with the monotonicity relation (23), we obtain and thus ν ⩾ −2δ, which is a contradiction to ν < −2δ ⩽ 0.
All in all, this theorem gives a rigorous yet conceptually simple criterion to check whether a given resolution guarantee is valid or not. by solving the boundary value problem (7)- (10). If this yields a negative value for ν, then the resolution guarantee holds true up to a measurement error of 0 ⩽ δ < − ν 2 .
Next, we formulate the corresponding theorem for case (18).

Theorem 2. The reconstruction of algorithm 2 fulfills the resolution guarantee if
Proof. The proof of part (i) of the resolution guarantee is analogous to the proof of part (i) in the theorem before. To show part (ii) of the resolution guarantee, assume that D = ∅ and D R = ∅. Then there must be an index s ∈ {1, 2 . . . , N} with Using the results from before, we obtain and thus ν ⩽ 2δ, which is a contradiction.

Theorem 3. The reconstruction of algorithm 3 fulfills the resolution guarantee if
Proof. First, let ω s ⊆ D and let a ∈ R M . In a body with interior Lamé parameters (λ B min , µ B min ), let ug be the corresponding displacements resulting from applying the boundary loadg. Based on the discrete version of the Neumann-to-Dirichlet operator (11), the variational formulation (4) as well as the associated bilinear form (5), we obtain where the last inequality holds due to lemma 2. Since ω s ⊆ D implies λ − λ B min ⩾ (c λ + λ 0 ϵ λ )χ ωs and µ − µ B min ⩾ (c µ + µ 0 ϵ µ )χ ωs , it follows in an analogous way that Hence, we obtain that For the proof of (ii), the reader is referred to the corresponding proof of theorem 1.
Finally, we present the theorem for the case (18). Proof. First, let ω s ⊆ D and let a ∈ R M . In a body with interior Lamé parameters (λ Bmax , µ Bmax ), let ug be the corresponding displacements resulting from applying the boundary loadg. As in the proof of the theorem before, we obtain This yields Hence, ω s will be marked, which shows that part (i) of the resolution guarantee. The second part is analogue to the proof of part (ii) from the theorem before.

Numerical simulations
We examine an elastic body (Makrolon) with possible inclusions (aluminium), where the corresponding Lamé parameters are given in table 1. We consider two different settings of test cubes (5 × 5 × 5 and 10 × 10 × 10) as well as two configurations of Neumann patches. Specifically, we apply boundary forces on five faces of the elastic body with either 5 × 5 or 10 × 10 Neumann patches on each face. Figure 1 shows exemplary the setting with 5 × 5 × 5 testcubes and 125 Neumann patches.
The forward problem is solved with COMSOL Multiphysics with LiveLink for MATLAB, where finite elements of degree 2 and tetrahedrons are used. In addition, we want to mention [14], which can be used as a tutorial for implementing finite element methods for inverse coefficient problems in elliptic partial differential equations.
Our simulations are based on noisy data. We assume that we are given a noise level η ⩾ 0 and set δ = η · Λ(λ, µ) F .

Example 1.
For our simulations we calculate the maximal noise η perturbing Λ(λ, µ) for different background error parameters ϵ λ and ϵ µ (see figure 2), based on theorem 3. In more detail, we take a look at the given partition shown in figure 1 and consider the inclusion contrast c λ = 0, c µ = 2 × 10 10 . It should be noted that the algorithm performs better by choosing c λ and c µ as close to the difference (λ D (x), µ D (x)) − (λ 0 , µ 0 ) as possible. Hence, selecting c λ = 0 makes the reconstruction worse in theory, however we will see in the results that the Lamé parameter λ and thereby the choice of c λ only marginally affects the results of the resolution guarantee. Further on, we apply the boundary loads in the normal direction on the Neumann patches. For solving the problem, we use the linearized monotonicity test in the form of algorithm 3, since λ < λ max and µ < µ max as denoted in table 1. All in all, our simulations will show us if and for which noise levels we obtain a resolution guarantee.
Note that patches in the above figures were only drawn, if a resolution guarantee exists for the tuple (ϵ λ , ϵ µ , η). Figure 2 tells us that the maximal η of approximately 1.413% is reached for ϵ λ = 0 = ϵ µ . The background error ϵ λ does not show much impact. Even for ϵ λ = 100%, we obtain a resolution guarantee. The maximal background error w.r.t. µ with ϵ λ = 0% is ϵ µ ≈ 7.692% at η = 0%. Remark 3. All in all, we conclude that the resolution guarantees depend heavily on the Lamé parameter µ and only marginally on λ. This is in accordance with the results in other papers, e.g. in [6].

Example 2.
Based on the result of Example 1, we change our configuration and set ϵ λ = 0% for a better comparability. The results are shown in figures 3-5, where we analyze the relation of ϵ µ (x-axis) and η (y-axis) with both values given in %. The considered numbers of  testcubes and Neumann patches are given in the caption of the figure. As expected, the smaller the background error ϵ µ can be estimated, the more noise on the data can be handled.
In figure 3, we deal with 5 × 5 × 5 = 125 testcubes and 125 Neumann patches as shown in figure 1. We can observe an approximately linear connection between ϵ µ and η showing that a resolution guarantee is given for all pairs (ϵ µ , η) on the black line and the gray area below for ϵ λ = 0%.
In figure 4, we change our setting and increase the number of testcubes to 10 × 10 × 10 = 1000, while simulating the reconstruction for the same 125 Neumann patches.
If we now compare figures 3 and 4, we see that for more testcubes, our method is less stable w.r.t. both ϵ µ and η. This behavior is expected since smaller pixels are to be reconstructed with the same amount of data from the Neumann patches. Nevertheless, we achieve a resolution guarantee, if the pair η, ϵ µ is located on the black line or the gray area below. The maximal noise on the data is given by η ≈ 0.200% for ϵ µ = ϵ λ = 0% and the maximal background noise for µ is given by ϵ µ ≈ 0.927% for ϵ λ = η = 0%. Increasing the resolution by using more Neumann patches is also possible as indicated in figure 5. This figure shows the set-up with 1000 testcubes, the same as in figure 4, but with 500 Neumann patches instead of 125. This increases both the stability regarding η as well es ϵ µ , however, the improvement is small. In fact, the maximal noise on the data is given by η ≈ 0.213% for ϵ µ = ϵ λ = 0% and the maximal background noise for µ is given by ϵ µ ≈ 0.942% for ϵ λ = η = 0%. For a better resolution guarantee, even more Neumann patches have to be used, but the numerical effort to do that will increases heavily.

Conclusion and outlook
Our main focus was the construction of conditions under which a resolution for a given partition can be achieved. Thus, our formulation takes both the background error as well as the measurement noise into account. The numerical simulations showed that for more testcubes our method is less stable w.r.t. ϵ λ , ϵ µ and η. This behavior is expected since more as well as smaller pixels are to be reconstructed with the same amount of data from the Neumann patches. As a result, the resolution guarantees depend heavily on the Lamé parameter µ and only marginally on λ. Finally, we want to remark that the algorithm is more stable w.r.t. ϵ λ , ϵ µ as w.r.t. η. All in all, our results are of special importance, when considering simulations based on real data, e.g. in [8] or in the framework of monotonicity-based regularization (see, e.g. [6]).

Data availability statement
Source data for figures are provided within the paper. The data that support the findings of this study are available upon reasonable request from the authors.