V F Sciuti
email: vfsciuti@gmail.com

R Vargas

R B Canto

F Hild

Pyramidal Adaptive Meshing for Digital Image Correlation Dealing With Cracks

Keywords:

wedge splitting test MgO hydration crack network computational cost

Introduction

Many materials remain functional with cracks before complete failure. Not only predicting such damage to avoid accidents but also its monitoring is a constant field of development. Since the introduction of Digital Image Correlation (DIC) in the early 1980s [START_REF] Sutton | Determination of displacements using an improved digital correlation method[END_REF][START_REF] Chu | Applications of digital-image-correlation techniques to experimental mechanics[END_REF], it has seen many developments [START_REF] Sutton | Computer vision-based, noncontacting deformation measurements in mechanics: A generational transformation[END_REF]. Among all applications of this technique [START_REF] Sutton | Advances in Two-Dimensional and Three-Dimensional Computer Vision[END_REF][START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF][START_REF] Hild | Digital image correlation[END_REF], the analysis of cracked bodies has received a lot of attention since the early days of DIC [START_REF] Peters | A study of dynamic near-crack-tip fracture parameters by digital image analysis[END_REF][START_REF] Mcneill | Estimation of stress intensity factor by digital image correlation[END_REF]. In the field of fracture mechanics, it was shown that stress intensity factors [START_REF] Mcneill | Estimation of stress intensity factor by digital image correlation[END_REF][START_REF] Roux | Digital image correlation and fracture: An advanced technique for estimating stress intensity factors of 2d and 3d cracks[END_REF], energy release rates [START_REF] Mathieu | Identication of interlaminar fracture properties of a composite laminate using local full-field kinematic measurements and finite element simulations[END_REF], parameters of cohesive zone models [START_REF] Abanto-Bueno | Experimental determination of cohesive failure properties of a photodegradable copolymer[END_REF][START_REF] Fedele | Identification of adhesive properties in glare assemblies by digital image correlation[END_REF][START_REF] Réthoré | Identification of a cohesive zone model from digital images at the micron-scale[END_REF]] could be determined via comparisons of measured displacements with simulations or available closed-form solutions. Some specialized numerical schemes were developed to capture the specific kinematics associated with the presence of cracks. For instance, integrated approaches [START_REF] Roux | Stress intensity factor measurements from digital image correlation: post-processing and integrated approaches[END_REF] were introduced when using Williams' kinematic series [START_REF] Williams | On the stress distribution at the base of a stationary crack[END_REF].

Since the introduction of FE-based DIC [START_REF] Broggiato | Adaptive image correlation technique for full-field strain measurement[END_REF][START_REF] Sun | Finite-element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF] and following on the development of the extended finite element method [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], the same type of kinematic enrichment was considered in extended DIC (X-DIC) [START_REF] Réthoré | From pictures to extended finite elements: Extended digital image correlation (X-DIC)[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF]. It was later generalized to 3D analyses of cracks monitored via computed microtomography [START_REF] Réthoré | Extended three-dimensional digital image correlation (X3D-DIC)[END_REF][START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF]. Another route is 'a priori' node splitting [START_REF] Roux | Mechanical assistance to DIC[END_REF][START_REF] Fagerholt | Investigation of fracture in small-scale SENT tests of a welded X80 pipeline steel using Digital Image Correlation with node splitting[END_REF][START_REF] Fagerholt | Measuring discontinuous displacement fields in cracked specimens using digital image correlation with mesh adaptation and crack-path optimization[END_REF][START_REF] Vargas | On the identification and validation of fracture mechanics models for curved crack in mortar[END_REF]. In both approaches (i.e., X-DIC and node splitting), it was shown that gray level residuals were very useful to adapt the level set function or the mesh to the crack path [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF][START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF][START_REF] Roux | Mechanical assistance to DIC[END_REF][START_REF] Fagerholt | Measuring discontinuous displacement fields in cracked specimens using digital image correlation with mesh adaptation and crack-path optimization[END_REF][START_REF] Vargas | On the identification and validation of fracture mechanics models for curved crack in mortar[END_REF] since they are evaluated pixel-wise.

The design of the mesh is crucial for FE-based methods, e.g., numerical simulation or DIC, especially in the presence of localized phenomena. The first FE-based DIC algorithms were pyramidal to enable for better convergence of the minimization scheme [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF][START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Fedele | Global 2D digital image correlation for motion estimation in a finite element framework: a variational formulation and a regularized, pyramidal, multi-grid implementation[END_REF]. Such h-refinement procedure has to be dealt with care since DIC analyses are inverse problems, namely, the finer the discretization, the higher the measurement uncertainty [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF][START_REF] Bornert | Assessment of digital image correlation measurement errors: Methodology and results[END_REF], yet the interpolation error is lower [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF]. Heterogeneous h-refinement is appealing when dealing with localized phenomena such as strained bands [START_REF] Wang | Mesh-Based Digital Image Correlation Method Using Non-Uniform Elements for Measuring Displacement Fields with High Gradient[END_REF] or cracks [START_REF] Baldi | Assessment of h-refinement procedure for global digital image correlation[END_REF][START_REF] Panin | Effect of the mesh size of the vector displacement field on the strain estimate in the digital image correlation method[END_REF] to better capture high strain gradients [START_REF] Hild | Damage measurements via DIC[END_REF]. As before, one critical aspect is related to the criterion for selecting elements to be refined. Several choices are possible, among them strains [START_REF] Baldi | Assessment of h-refinement procedure for global digital image correlation[END_REF][START_REF] Panin | Effect of the mesh size of the vector displacement field on the strain estimate in the digital image correlation method[END_REF] and gray level residuals (GLR) [START_REF] Wang | Mesh-Based Digital Image Correlation Method Using Non-Uniform Elements for Measuring Displacement Fields with High Gradient[END_REF][START_REF] Baldi | Assessment of h-refinement procedure for global digital image correlation[END_REF][START_REF] Vargas | On the identification and validation of fracture mechanics models for curved crack in mortar[END_REF].

Another route is given by p-refinement, which is generally easier to implement since it does not affect the mesh connectivity. This procedure has already been used in DIC [START_REF] Ma | Mesh-based digital image correlation method using higher order isoparametric elements[END_REF][START_REF] Li | Local Hermite (LH) Method: An accurate and robust smooth technique for high-gradient strain reconstruction in digital image correlation[END_REF], increasing the resolution without reducing the element size. Such implementations can be self-adaptative [START_REF] Wittevrongel | A self adaptive global digital image correlation algorithm[END_REF][START_REF] Wittevrongel | C -continuity in Digital Image Correlation: Implementation and Validation of C -1 , C 0 and C 1 Algorithms[END_REF]. Kleinendorst et al. [START_REF] Kleinendorst | On the use of adaptive refinement in isogeometric digitalimage correlation[END_REF] implemented p-refined DIC using Non Uniform Rational Basis Splines (NURBS).

Hereafter, the numerical procedures are first presented, i.e., DIC for the measurements, the Lagrange multipliers in the DIC scheme to deal with the hanging nodes, and the proposed adaptive meshing procedure. Then, three cases are shown to illustrate the Adaptative Meshing (AM) procedure. For the artificial case (#1.1 to #1.3), three images were created to compare the gains of applying AM with different crack densities. Then two experiments are analyzed, namely, a Wedge Splitting Test (WST, Case #2) dealing with a dominant crack, and curing and drying of an MgO refractory castable (Case #3) leading to the formation of a crack network, where the effects of AM over the DIC results are highlighted. For all cases, the locally refined meshes are compared with their analog uniform mesh (UM) with the smallest element size.

Numerical procedures 2.1. Digital Image Correlation

DIC is a full-field measurement technique based on the optical flow, i.e., the gray level of any pixel location in the reference configuration is equal to the gray level in the deformed state image at location + ( ) instead of .

The displacement field is determined by minimizing the gray level residual Φ 2 over the selected domain Ω

Φ 2 = ∑ Ω [ ( ) -( + ( ))] 2 (1)
Searching for pixel-wise is an ill-posed problem, and thus it is regularized by searching in a subspace made of trial fields

Ψ ( ) = ∑ ( ) (2)
where are nodal displacements gathered in the column vector { }. The minimization of Φ 2 ({ }) is carried out with an iterative Gauss-Newton scheme that leads to linear systems

[ ]{ } = { } (3)
where { } gathers the updates to the nodal displacement for the current iteration, and the other terms are the Hessian matrix

= ∑ Ω ( ⋅ ∇ ∇ ∇ )( )( ⋅ ∇ ∇ ∇ )( ) (4)
and the residual vector

ℎ = ∑ Ω [ ( ) -( + ̃ ( ))]( ⋅ ∇ ∇ ∇ )( ) (5)
where ̃ is the current estimate of the nodal displacement component .

Depending on the choice of domain Ω different correlation procedures are obtained. When Ω reduces to a small interrogation window, the previous setting corresponds to that of local DIC [START_REF] Sutton | Image correlation for shape, motion and deformation measurements: Basic Concepts, Theory and Applications[END_REF][START_REF] Sutton | Computer vision-based, noncontacting deformation measurements in mechanics: A generational transformation[END_REF]. Conversely, Ω may cover the whole region of interest (ROI). In that case case, global approaches are obtained [START_REF] Wagne | Spectral approach to displacement evaluation from image analysis[END_REF]. Among many choices for kinematic bases, finite element descriptions are very versatile [START_REF] Sun | Finite-element formulation for a digital image correlation method[END_REF][START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]. In the present case, three noded (T3) triangles were considered (i.e., T3-DIC) [START_REF] Leclerc | Integrated digital image correlation for the identification of mechanical properties[END_REF]).

In the present studies, very fine meshes will be constructed. At one stage, the discretization will be too fine for standard FE-based DIC to converge. It is however possible to supplement DIC with mechanical input, namely, via mechanical regularization [START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis of fractured samples: The equilibrium gap method as a mechanical filter[END_REF][START_REF] Tomičevć | Mechanics-aided digital image correlation[END_REF] based upon the equilibrium gap method [START_REF] Claire | A finite element formulation to identify damage fields: The equilibrium gap method[END_REF]. In linear elasticity, the discretized equilibrium equations read

[ ]{ } = { }, (6)
where [K] is the stiffness matrix, and { } the vector of nodal forces. If the displacement field does not satisfy equilibrium, force residuals { } arise

{ } = [ ]{ } -{ } (7)
In the absence of body forces, the nodal forces of inner nodes must vanish, which makes the goal of the equilibrium gap approach to minimize [44]

Φ 2 = ‖{ }‖ 2 = { } ⊤ [ ] ⊤ [ ]{ } ( 8
)
where ⊤ is the transpose operator, and Φ 2 the summation of the L2-norm of all equilibrium gaps for inner nodes and traction-free boundaries (i.e., Neumann nodes). For Dirichlet boundary nodes, the discrete Laplace-Beltrami operator is considered instead [START_REF] Mendoza | Complete mechanical regularization applied to digital image and volume correlation[END_REF]. In the present case, the previous regularization was applied to all nodes.

To simultaneously reduce the correlation residuals Φ 2 and the equilibrium gap Φ 2 , the total cost function Φ 2

(1 + )Φ 2 = Φ2 + Φ2 (9)
is minimized, where is the weight that defines the scale associated with the equilibrium gap, Φ2 and Φ2 normalized cost functions. The importance of normalization is to convert the residuals into dimensionless quantities [START_REF] Leclerc | Voxel-scale digital volume correlation[END_REF][START_REF] Mendoza | Complete mechanical regularization applied to digital image and volume correlation[END_REF]. It is carried out by considering a trial displacement field in the form of a plane wave ( ) = exp( ⋅ ), where is the amplitude and the wave vector. The normalized cost functions become

Φ2 = Φ 2 { } ⊤ [ ]{ } , Φ2 = 2 { } ⊤ [ ] ⊤ [ ]{ } . ( 10
)
The wavelength dependence of { } ⊤ [ ] ⊤ [ ]{ } is of the fourth-order wrt. ‖ ‖, and the quantity { } ⊤ [ ]{ } is independent of the wave number ‖ ‖. The weight reads

= (‖ ‖ ) 4 (11)
where is the so-called regularization length (for Φ 2 ). Therefore, the weight put on the cost function is proportional to raised to the power 4. The value chosen for should be small enough to fit in small ROIs but large enough to avoid divergence of the minimization scheme [START_REF] Leclerc | Voxel-scale digital volume correlation[END_REF]. For this work, whose objective is to monitor cracks, the value of needs to be as small as possible, otherwise displacement jumps in the crack path will be spread over adjacent elements.

Instead of solving Systems (3), the new systems become

[ ]{ } = { } (12)
with

[ ] = [ ] + { } ⊤ [ ]{ } { } ⊤ [ ] ⊤ [ ]{ } [ ] ⊤ [ ] (13)
and

{ } = { } - { } ⊤ [ ]{ } { } ⊤ [ ] ⊤ [ ]{ } [ ] ⊤ [ ]{̃ } (14)
The iterations are stopped when the L2-norm of the vector containing the root mean square (RMS) of the horizontal and vertical displacement corrections was less than 10 -3 px.

Adaptive meshing

Refinement strategy

In the present work, the mesh size will be lowered within FE-based DIC procedures (i.e., heterogeneous h-refinement).

During a local h-refinement, some nodes located in the middle of an edge between two elements (e.g., one subdivided and others not), have different numbers of connectivities. To illustrate this phenomenon, let us consider a structured mesh with two triangular elements (Figure 1 are not hanging nodes because they are on the outer edges. For localized h-refinement, at least one node is created on the edge that separates the subdivided domain and the intact one. The presence of hanging nodes makes the mesh non-conformal, which can induce numerical artifacts. A natural alternative to deal with hanging nodes is to divide the elements creating a new node at the centroid of a parent element and connect it to the element vertices (Figure 1(d)). However, this approach generates slender elements that may compromise the calculations. To overcome this drawback, Baldi and Bertolino [START_REF] Baldi | Assessment of h-refinement procedure for global digital image correlation[END_REF] reconstructed the mesh using the new nodes provided by local refinement. An alternative solution to remeshing strained regions (in the present case cracks) is proposed herein, which uses Lagrange multipliers to link the displacement of hanging nodes based on those of the adjacent (also called parent) nodes (nodes 1 and 4 in Figure 1) and the shape function of the considered element.

If an element edge is divided at its middle point, the displacements of the hanging node is the average of the displacement of its parent nodes. If the edges are not split in their middle, weights related to the distance from the hanging node to each parent node should be considered and the shape functions. Such an approach can be easily implemented in a DIC algorithm and does not require the DIC Hessian to be changed, because only an auxiliary (Lagrange) matrix [ ] is added. The [ ] matrix is pre-multiplied by the ratio of Frobenius norms of [ ] and [ ],

, to avoid degrading the full system conditioning. Once all the desired constraints are accounted for, the augmented

Hessian [ ] becomes [ ] = ⎡ ⎢ ⎢ ⎣ [ ] [ ] ⊤ [ ] [ ] ⎤ ⎥ ⎥ ⎦ (15)
The vector of unknowns { } now gathers all nodal displacements { } and Lagrange multipliers { }, and the augmented residual vector { } reads

{ } ⊤ = {{ } ⊤ { } ⊤ } (16)
resulting in the augmented system to solve as

[ ]{ } = { } (17)
In experimental cases, defining which elements need to be divided is not straightforward. There are several choices, among them the maximum eigen strains [START_REF] Hild | Damage measurements via DIC[END_REF][START_REF] Baldi | Assessment of h-refinement procedure for global digital image correlation[END_REF][START_REF] Panin | Effect of the mesh size of the vector displacement field on the strain estimate in the digital image correlation method[END_REF] and the gray level residuals [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF][START_REF] Roux | Mechanical assistance to DIC[END_REF][START_REF] Fagerholt | Measuring discontinuous displacement fields in cracked specimens using digital image correlation with mesh adaptation and crack-path optimization[END_REF][START_REF] Wang | Mesh-Based Digital Image Correlation Method Using Non-Uniform Elements for Measuring Displacement Fields with High Gradient[END_REF][START_REF] Hild | Damage measurements via DIC[END_REF][START_REF] Vargas | On the identification and validation of fracture mechanics models for curved crack in mortar[END_REF]. It is worth noting that the eigen strains are computed as the average over each element. Conversely, the gray level residuals

( ) = ( ) -( + ( ))
, where are the measured nodal displacements (i.e., at convergence of the DIC code), are defined at the pixel scale. In both cases, a user-defined threshold will be needed.

Refinement Criterion

Since cracks are to be detected and quantified, it is proposed to introduce a refinement criterion explicitly accounting for displacement jumps and their associated measurement uncertainty. In the presence of cracks, the mean strain tensor in each element is expressed as [START_REF] Vakulenko | Continuum theory of medium with cracks, Isv[END_REF] = 1

|Ω | ∫ Ω ( )d + 1 |Ω | ∫ Γ ( ) ⊠ ( ) d (18)
where denotes the displacement jump, the normal to the cracked surface Γ , ⊠ the symmetrized tensorial product, the strain tensor of the uncracked matrix, and the curvilinear abscissa along the crack path. In the following, it is assumed that the elastic strain levels are negligible with respect to the singular contribution provided by cracks.

Consequently, the mean strain tensor reduces to

≈ 1 |Ω | ∫ Γ ( ) ⊠ ( ) d (19)
The mean crack opening displacement per element is defined as

⊠ = 1 |Γ | ∫ Γ ⊠ d (20)
where denotes the mean crack normal. Further, it is assumed that the main opening contribution is associated with the mode I regime

= ( 21
)
such that then corresponds to the eigen direction associated with the average maximum eigen strain 1 . The latter is then related to the mean crack opening displacement by

|Ω | 1 = |Γ | (22)
It is worth noting that |Ω |∕|Γ | scales with the element size , provided only one crack is present in each element.

This hypothesis is likely to be true for small element sizes. Consequently, the mean crack opening displacement is approximated by

= 1 (23)
The mean crack opening displacement threshold * will be set to classify elements as damaged (i.e., ≥ * ).

It is proposed to relate * to its measurement uncertainty (i.e., typically 5 times or more its standard uncertainty ( )). According to Equation ( 23), the uncertainty ( ) is related to the standard strain uncertainty ( 1 ) by ( ) = ( 1 ). Further, the strain uncertainty ( 1 ) is dependent on the element size

( 1 ) = (24)
and characterizes the compromise to be carried out between the uncertainty and the discretization level (i.e., spatial resolution). Typical values of the power range between -1 and -2 (i.e., when the measurement uncertainties are random and controlled by white Gaussian noise [START_REF] Hild | Comparison of local and global approaches to digital image correlation[END_REF]). The constant is dependent on the contrast in the picture and the noise level [START_REF] Besnard | Finite-element" displacement fields analysis from digital images: Application to Portevin-Le Chatelier bands[END_REF]. The standard uncertainty ( 1 ) (and therefore ( )) is obtained by performing an uncertainty quantification prior to each analyzed test. Consequently, the threshold

* becomes * ( ) = +1 (25)
where is a proportionality constant between the standard uncertainties. Thus, the threshold * explicitly depends on the discretization level. This observation means that the threshold * has to be changed for each refinement step to account for the reduction of element size (or equivalently to capture elements for which the crack opening displacement is greater than * ). During the analysis the elements were divided by two and, consequently, the threshold was multiplied by two in each refinement step.

By using Equation ( 23), a strain threshold * ( ) = * ( )∕ can also be defined. With the previous proposition, it is concluded that * is also length dependent. This proposition departs from the use of a unique strain threshold [START_REF] Wang | Mesh-Based Digital Image Correlation Method Using Non-Uniform Elements for Measuring Displacement Fields with High Gradient[END_REF][START_REF] Baldi | Assessment of h-refinement procedure for global digital image correlation[END_REF][START_REF] Panin | Effect of the mesh size of the vector displacement field on the strain estimate in the digital image correlation method[END_REF] when multiple refinement steps are to be performed. Both thresholds (i.e., * 1 and * ) introduced herein relate to the measurement uncertainties of the associated quantities (i.e., 1 and

). Throughout this work the discussion will be focused on the use of * and its dependence on the parameter .

Refinement Steps

It is proposed to start off with a coarse mesh and to use the previous criterion to refine the 'damaged' elements (i.e., ≥ * ( )), which will better capture the kinematics and morphology of the crack locations while possibly remaining with low computational cost. The complete procedure consists of the following steps:

1. perform a DIC analysis using a coarse uniform mesh (UM);

2. choose an image for crack evaluation based on the displacement jump ( ) ;

3. list the damaged elements for this image pair and considered discretization;

4. execute one subdivision in the damaged elements (i.e., each selected element generates four new ones);

5. identify the hanging nodes and respective parent nodes;

6. construct the Lagrange matrix and perform DIC analyses for this image pair using the locally refined mesh;

7. repeat items 3 to 6 until the smallest elements achieve the size desired by the user;

8. perform a DIC analysis for all images using the optimized mesh.

This approach was implemented within the Correli 3.0 framework [49] and is described in the flowchart of Figure 2.

It is worth noting that, in the present implementation, the end of the adaptive procedure is a user decision, utilizing the size of the smallest element after division as a stop criterion. The smallest element size may be guided by the highest measured Mean Crack Opening Displacement (MCOD) using UMs. If even smaller elements were envisaged (i.e., down to one pixel size), damage could be used along with mechanical regularization [START_REF] Hild | Damage measurements via DIC[END_REF].

Artificial Cases (#1.1 to #1.3)

This statement is specific to the artificial images .

The artificial images of 200 × 200 pixel definition shown in Figure 3 are analyzed to illustrate the benefits of the adaptive meshing proposal. Such pictures may correspond to the gray level residuals at the end of a DIC run [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF][START_REF] Réthoré | Extended digital image correlation with crack shape optimization[END_REF][START_REF] Rannou | Three dimensional experimental and numerical multiscale analysis of a fatigue crack[END_REF][START_REF] Roux | Mechanical assistance to DIC[END_REF][START_REF] Fagerholt | Measuring discontinuous displacement fields in cracked specimens using digital image correlation with mesh adaptation and crack-path optimization[END_REF][START_REF] Vargas | On the identification and validation of fracture mechanics models for curved crack in mortar[END_REF]. In the present cases these images were built to demonstrate the method (i.e., the threshold for binarization is known), the binary images remove the user-dependence to identify the cracked elements. To represent practical situations, a single crack (Case #1.1, Figure 3(a)) is assumed with a sinuous front. An additional case consists of a major macrocrack with some branches (Case #1.2, Figure 3 These images are first binarized, in which one represents a cracked state, and zero otherwise. Then, starting from the simple mesh shown in Figure 4(a), if any pixel that represents a crack is contained in an element, it is divided into four smaller elements (see Figure 1(c)), and this procedure is iteratively repeated. The starting element is of the image size, with a mean element size 0 = 141 px, obtained as the square root of the element area ( √ 0 ). Hereafter, the nomenclature it will be used for the mean element size of iteration it for the reference case where all elements are divided, with it=0 the initialization. It is worth noting that it is also related to the smallest element size in the AM methodology proposed herein. Seven iterations were chosen for the smallest element to remain greater than 1 px (i.e., In the case of the reference meshes where all elements were subdivided, the DIC Hessian matrix [ ] was also computed for comparison purposes. For three refined meshes, the augmented Hessians [ ] were assembled to also consider the new constraints using Lagrange multipliers. The sizes of the Hessian matrices are shown in Table 1 for all three cases. In virtually all cases, there is a clear benefit of the AM procedure, especially for very fine discretizations. The relative gains in terms of size of the sparse matrices are reported in Table 2 when compared to uniformly refined solution, calculated as one minus each Hessian size divided by its reference for a given iteration. There were negative values in some iterations, only for the crack network (i.e., because more hanging nodes were created than edges that had not to be split). The gains of this refinement methodology are inversely proportional to the crack density. In cases with few cracks, it is interesting to have fewer elements where no localized phenomena are taking place, for a more efficient calculation and no loss of accuracy. However, for very dense crack networks, this approach may even be worse in some cases than a full mesh refinement (e.g., when ≈ 10 px). If very small elements are required (i.e., about the pixel size), it is likely that this method allows for gains in most cases. One experimental setup that allows for stable crack propagation is the so-called Wedge Splitting Test (WST) [START_REF] Tschegg | Prüfeinrichtung zur Ermittlung von bruchmechanishen Kennwerten sowie hiefür geeignete, Prüfkörper, Austrian Pat[END_REF][START_REF] Harmuth | Stability of crack propagation associated with fracture energy determined by wedge splitting specimen[END_REF]. It is used hereafter to exemplify the benefits of the method proposed herein when not too many cracks are suspected to occur. In the WST, a vertical force is transmitted to a pre-notched sample using a wedge and cylinders to prescribe a horizontal (mode I) splitting load. DIC was used to analyze WST experiments via subset-based DIC [START_REF] Belrhiti | Investigation of the impact of micro-cracks on fracture behavior of magnesia products using wedge splitting test and digital image correlation[END_REF][START_REF] Dai | Determination of the fracture behaviour of MgO-refractories using multi-cycle wedge splitting test and digital image correlation[END_REF][START_REF] Dupré | Detection of cracks in refractory materials by an enhanced digital image correlation technique[END_REF][START_REF] Khlifi | Investigation of microstructure-property relantionships of magnesia-hercynite refractory composites by a refined digital image correlation technique[END_REF] or using T3-DIC [START_REF] Vargas | Analysis of wedge splitting test on refractory castable via integrated DIC[END_REF][START_REF] Vargas | Comparison of two full-field identification methods for the wedge splitting test on a refractory[END_REF] based on structured meshes. Further, Ref. [START_REF] Vargas | Analysis of a castable refractory using the wedge splitting test and cohesive zone model[END_REF] discussed how to calibrate a cohesive zone model for such experiment, which required a mesh to be adapted in the crack propagation zone.

(b)). Last, a crack network is considered (Case #1.3,
The studied material was a commercial castable refractory DD40, made by IBAR, that can be utilized in catalytic cracking units of oil industries [START_REF] Stafuzza | Application of the digital image correlation technique to the fracture energy method for refractory concretes[END_REF][START_REF] Miyaji | The coke effect on the fracture energy of a refractory castable for the petrochemical industry[END_REF][START_REF] Pereira | A study about the contribution of the -phase transition of quartz to thermal cycle damage of a refractory used in fluidized catalytic cracking units[END_REF]. Wedge Splitting Tests (WSTs), which were performed along the guidelines proposed in Ref. [START_REF] Ribeiro | Projeto, adaptação, instalação e testes preliminares para um sistema de medida de energia de fratura de materiais cerâmicos pelo método da cunha[END_REF], were first analyzed via integrated DIC [START_REF] Vargas | Analysis of wedge splitting test on refractory castable via integrated DIC[END_REF]. The specimen used in the test reported hereafter was similar with the sole difference being the absence of lateral grooves to prescribe straight propagation. A parallelepipedic sample of size 100 × 100 × 75 mm 3 was molded with a pre-notch to insert the loading parts. The reference image and the initial mesh (in yellow) used in the procedure developed herein are shown in Figure 6(a). In the same figure, the wedge, both cylinders, and blocks that were used to transmit the load horizontally, are also highlighted in white dashed lines. A random pattern was created by spraying black and white paints to increase the image contrast.

In the present experiment, there was no groove to guide the crack so that a tortuous crack path may be observed with possibly some bifurcations. A pre-load of 40 N was applied before starting the test, to adjust the loading parts (Figure 6(b)). The experiment was controlled with a constant actuator velocity of 80 µm/min and was stopped at 5 % of the maximum load. A 0.1 Hz image acquisition rate was used (i.e., a total of 354 pictures was available). Only the odd images ranging from 5 to 257 (i.e., down to 10% of the peak load) were analyzed to avoid lengthy discussions on the last part of the WST, resulting in a total of 126 images. The 30th analyzed frame (blue asterisk in Figure 6(b)) will be used throughout the discussion.

The DIC parameters used hereafter are reported in Table 3. The five images acquired during the pre-loading step were used to compute the measurement uncertainties for UMs. The temporal variances of the maximum eigen strain were computed in each element, and their root mean is shown in Figure 7(a) as a function of the inverse element size -1 . The expected power law dependence described by Equation ( 24) is observed with = -1.66. This value is close but not equal to -2 (i.e., some other sources of uncertainties occurred in this experimental configuration). The standard MCOD uncertainty is reported in Figure 7(b).

A similar agreement is observed in comparison to Equation [START_REF] Fagerholt | Investigation of fracture in small-scale SENT tests of a welded X80 pipeline steel using Digital Image Correlation with node splitting[END_REF]. The black stars in Figure 7 refer to the uncertainties computed for the optimal AM ( = 3), and the procedure to achieve this adapted mesh is described hereafter. The black stars mark the uncertainties for the optimal AM with = 3

For each frame, the standard deviation of the gray level residuals (GLR) defined as

( ) = ( ) -( + ( )) (26)
was computed and normalized by the dynamic range of the reference picture (i.e., Δ = max ROI ( ) -min ROI ( )),

and the effect of the element size using UMs are shown in Figure 8. In the beginning of the test, the gray level residuals were independent of the element size. Their overall level was related to acquisition noise, and the fluctuations due to small lighting variations. As the crack initiated and propagated, the residuals became more dependent on the element size (i.e., the larger the element size, the higher the residuals). This feature is explained by the lack of degrees of freedom to properly describe the kinematics of the crack, which is improved by finer meshes. This effect corresponds to interpolation errors [START_REF] Réthoré | Shear-band capturing using a multiscale extended digital image correlation technique[END_REF]. Further, for small element sizes, the gain in discretization is less pronounced.

This phenomenon is due to the regularization scheme used herein, which dampens out sharp displacement gradients over the area defined by the regularization length. GLR fields are shown in Figure 9 for the frame #30. The GLR field of 64 px UM was set as a reference and the GLR fields of the other UMs were subtracted from it to highlight their differences. As expected, the crack is better accounted for as the element size decreases. However, using a very fine mesh implies not only higher computational cost, but also more difficult convergence of the analysis due to the reduced number of pixels inside each element. This issue can be addressed via regularization strategies (e.g., mechanical regularization [START_REF] Réthoré | An extended and integrated digital image correlation technique applied to the analysis fractured samples[END_REF][START_REF] Tomičević | Mechanics-aided digital image correlation[END_REF] as used in this work). Last, let us note that even with very fine meshes, there are still high gray level residuals in the cracked region, which were induced by the hypothesis of continuous displacement fields since no node splitting strategy was used herein [START_REF] Roux | Mechanical assistance to DIC[END_REF][START_REF] Fagerholt | Investigation of fracture in small-scale SENT tests of a welded X80 pipeline steel using Digital Image Correlation with node splitting[END_REF][START_REF] Fagerholt | Measuring discontinuous displacement fields in cracked specimens using digital image correlation with mesh adaptation and crack-path optimization[END_REF][START_REF] Vargas | On the identification and validation of fracture mechanics models for curved crack in mortar[END_REF].

One can notice gray level residuals almost spread over the entire ROI, which does not change with the element size.

The residuals are due to the acquisition noise and to fluctuations due to small lighting variations.

The mean crack opening displacement (MCOD) fields for frame #30 computed using different element sizes are shown in Figure 10. MCOD fields were less affected by the acquisition noise in comparison to the grainy GLR fields, pointing MCOD as a better metric for crack detection. Further, crack details (e.g., bifurcations) are better captured by small element sizes. However, it is essential to stress that uniform meshes with small element sizes require a higher computational cost. This fact motivates the use of meshes whose refinement is "MCOD-guided" as discussed hereafter.

The refinement from 64-px UM down to 8-px UM showed a visible benefit to describe the crack and its ramification.

However, the 4-px UM does not present significant improvements when compared to 8 px UM (even with a four times Figure 9: Normalized gray level residual field and differences for the frame #30 of the analysis using UMs with different element sizes. For the sake of visualization the results from finer meshes were subtracted from the coarser mesh result (64 px). The color bar range was reduced to highlight the patterns in the field larger Hessian, see Table 4). This observation is due to the regularization length utilized herein. Considering all the 228 previous points, the best compromise for the tested UMs is provided by the 8-px discretization. The following analyses are devoted to adapted meshes. In Figure 11, the map of mesh sizes are reported at the end of the adaption steps for 5 different levels of (i.e., the ratio between the measurement uncertainty and the threshold used to detect damaged elements). When the threshold is very high (i.e., = 12, 6), only the zones surrounding the main crack are refined. Conversely, when is too low (e.g., = 1.5), the mesh is refined in areas that are not cracked.

From these first analyses, it is concluded that = 3 leads to the best compromise. This observation is confirmed when analyzing the maps of MCOD reported in Figure 12. For the first two levels 235 (i.e., = 12, 6), the presence of secondary cracks is smeared over some coarse elements. Conversely, when is too 236 small (i.e., = 1.5), the additional refinement does not bring any new insights into the crack pattern. Again = 3 is 237 a good compromise between capturing the details around the main crack and the related computational cost. The acquisition noise in the GLR fields for the analyses using uniform meshes (Figure 9) is also observed when using adaptive meshes (Figure 13). This result was expected because acquisition noise does not depend on the element size. It is possible to depict the full developed crack, even on frame #30, when the crack has not developed that far (see Figure 9), due to residuals caused by mesh refinement. This result makes the use of GLR field difficult to analyze crack propagation when using the proposed adaptive mesh, and reinforces the use of MCOD fields for this application.

However, the residuals are concentrated around the cracked regions, and the localized refinement does not generate fake cracks in the boundaries of domains with different element sizes (Figure 11). It is interesting to note that the crack branches seen in the MCOD field are not visible in any of the GLR fields, thereby indicating that the corresponding displacement levels remained very small (i.e., well below the pixel level). The normalized standard deviation GLR for each frame of the different adapted meshes are shown in Figure 14 and compared to the uniform mesh (designated as UM in the sequel) results from Figure 8. The standard deviations of the GLR at the beginning of the experiment are virtually identical, showing that all adapted meshes perform the same way, which is a further indication that no crack initiated. For frames closer to the end of the experiment, where the residuals due to crack openings are higher than those due to acquisition noise and lighting fluctuations, the AM with a threshold = 3 leads to GLR close to the UM 16 px discretization. The = 1.5 AM results are close to the 8 px UM, even with this AM containing varying element sizes (Figure 11(c-d)). The previous observations mean that the same quality of results is achieved using less degrees of freedom (Table 4).

The smallest standard deviation (STD) of the GLR for the last frame was found for the 4-px UM (i.e., it is equal to 2.95%). Some 4 px elements are also present in AM for ranging from 6 down to 1.5 (Figure 11). The STD of GLR levels are close to those of the 8-px UM, mainly for the AM with = 1.5. The size of the Hessians related to ≤ 3 are 6 and 4 times smaller than those of 4-px UM, respectively. This result indicates that the main source of GLR increase is the lack of degrees of freedom in the cracked region. It demonstrates that a localized refinement can considerably improve the description of cracks without huge increase in computational cost. Further, the computational cost is also represented by the computation time (CT) for each analysis, which was obtained by summing the elapsed time during system solving in each iteration. Both CT and average iteration per image analyzed are shown in Table 4. The AM with = 3 is the best compromise between a fine discretization and a low computational cost since the refinement is concentrated in the whole cracked region (Figure 11 and Table 4). Higher values (i.e., = 6, 12) do not refine enough the cracked region and do not contain 4-px elements that help to describe the crack path. Conversely, when = 1.5, spots where there is no crack are refined. In terms of MCOD maps, the AM with = 3 leads to the same type of resolution (Figure 12) as the 8-px UM (Figure 10) but with a smaller Hessian (see Table 4). Interestingly, even if the conditioning of the system is degraded by the use of Lagrange multipliers, clearly visible in the CT of similar Hessian sizes, the problem needed less iterations to converge due to the description of the crack region. Repeating the uncertainty analyses using only the initial images and the AM with = 3 led to standard uncertainties for the maximum eigen strain equal to 7.5 × 10 -5 , and for the MCOD to 5 × 10 -2 µm (see Figure 7).

One last question arises about subsequent crack propagation after frame #30 since the mesh does not account for the fully propagated stage. Figure 15 summarizes the results after adapting the mesh for the very last frame using = 3. The element size distribution now has many 4 px elements concentrated along the crack path. Figure 15(b-c) reports the MCOD maps using this new AM, for the 30-th and the last frames, respectively. The fine discretization in the whole crack path leads to noisier MCODs in earlier frames in regions that are probably not cracked yet (i.e., below the mid-height for frame #30, see Figures 10 and12). Figure 16 shows the GLR ratio between the first AM with = 3 and the one defined with the last frame. Both residuals have the same levels until frame #30, and from that point on their ratio rises up to 17% since the first mesh did not account for crack propagation after this point. The AM at the end of the test shows noisier MCODs for earlier frames but with lower GLR. Such observations call for gradual mesh adaptions throughout the frames. However, more meshes (and consequently Hessian matrices) would need to be stored, and the benefits on the computational cost would need to be tested in future works.

Damage due to MgO hydration in refractory castable (Case #3)

Magnesia (MgO) is used in some refractory castable formulations [START_REF] Salomão | The role of hydraulic binders on magnesia containing refractory castables: calcium aluminate cement and hydratable alumina[END_REF]. MgO hydrates during curing and drying, which may damage the material. In previous works, a climatic chamber was developed to evaluate such mechanism using DIC [START_REF] Sciuti | Digital Image Correlation applied to in situ evaluation of surface cracks upon curing of MgO-containing refractory castables[END_REF]. A cubic specimen was kept in the climatic chamber at 50°C and 50% relative air humidity ( The maximum eigen strain ( 1 ) and MCOD uncertainties were calculated using the ten images acquired at the beginning of the experiment, before crack initiation. The temporal variance of 1 and MCOD were computed for each element, and their root means are shown in Figure 18 as functions of the inverse element size -1 . The power is now equal to -1.47, which is smaller than the first experiment. It shows that uncertainty quantification should be performed for each test as the thresholds * and * may vary from one experimental configuration to another one. The black stars in Figure 18 refer to the uncertainties computed for the optimal AM ( = 1.5) Figure 19 shows the standard deviation of the GLR when normalized by the dynamic range of the reference picture.

The effect of the element size for a uniform mesh (UM) is first investigated. At the beginning of the test, the overall GLR level is mainly related to acquisition noise. The fluctuations are also due to lighting variations and the environment of the climatic chamber. They are more pronounced in this second experiment, as highlighted by the spikes. Corrections to the gray level fluctuations were suggested for such cases [START_REF] Sciuti | On the benefits of correcting brightness and contrast in global digital image correlation: Monitoring cracks during curing and drying of a refractory castable[END_REF]. However, the correction was not used herein to evidence the increase in noise sensitivity with a decrease of the element size (i.e., increase in spike amplitudes for small elements). Further, the GLR for each UM are very close until frame #25, for which cracks start to show measurable MCOD. For the following frames, the GLR curves are distinguishable. It is caused by the lack of degrees of freedom (DOFs) to describe the crack network with coarse discretizations. This contribution to the GLR is less pronounced as the element size decreases (and the number of DOFs increase). The GLR fields for frame #38 are shown in Figure 20. The field for the 64-px UM was set as a reference and the GLR fields of the other UMs were subtracted from it to highlight their differences as in Case #2 (see Figure 9).

Acquisition noise in the undamaged parts is identical for all UMs, and the different meshes sizes affect the regions around the cracks as expected from Case #2 results. The GLR caused by the crack network are thinner for UMs with 8-or 4-px elements. However, the 4-px UM was very sensitive to lighting changes (Figure 19), which hindered cracks to be properly detected and quantified. Such observation makes the 8-px mesh a better choice. The corresponding GLR fields are shown in Figure 24 for the studied thresholds . The acquisition noise remains present in all cases since it is not dependent on the element size. Similarly to the previous results, the crack pattern is more visible in the GLR maps of higher thresholds due to interpolation errors related to fewer DOFs.

V. F. The normalized standard deviation of GLR comparing UMs and AMs is shown in Figure 25 for the whole experiment. As before, the smallest residuals are observed for the 4 px UM, and AM with = 3 and 1.5, which are close to 16 px and 8 px UMs, respectively. The similarity of the curves before frame #38 further evidences that no (or very little) cracks propagated until this point. From frame #38 onward, higher residuals are related to higher gray levels induced by the newly created cracks, and consequently the normalized GLR follows the amount of DOFs to allow for such complex kinematics associated with the crack network to be captured. An additional contribution is given by the fact that mesh adaption was not performed again. The Hessian size for each case shown in Figure 25 is reported in Table 5 along with the normalized GLR for the last frame, as well as the computation time (CT), and the average iteration number. As pointed out for the numerical example, the benefits of using AMs are inversely proportional to the crack density. Further, the local refinement of AMs did not reduce the average iteration number in Case #3 as was noted in Case #2. The main benefit comes from using the AM with = 1.5, with similar GLR and Hessian size to the 8 px UM but containing few 4 px elements. The uncertainty analyses using only the initial images and AM with = 1.5 led to standard uncertainties for the maximum eigen strain and MCOD equal to 2 × 10 -4 and 0.9 µm, respectively (see Figure 18). As in the previous case, one last point of interest arises from the frame used for mesh adaption. Figure 26 shows the results of an AM with = 1.5 for the very last frame of the test. The crack path is finely described with very small elements (Figure 26(a)), with considerably more 4-px elements than for the 38-th frame (Figure 22). The corresponding MCOD fields are compared in Figure 26(b-c). A similar trend (to the WST case) of noisier results for the earlier frame, and a cleaner description for the last frame. The fluctuations seen for frame #38 are due the fine discretization related to crack openings and not its initiation. Although the crack network was already visible in earlier frames, its openings were very small. Last, the ratio between GLR for the AM for the 38-th frame and the last picture is shown in Figure 27. For earlier frames, a coarser discretization (provided by the 38-th frame) was slightly better, and as cracks subsequently opened, the second mesh provided little gain. This observation is further evidence that the crack network was already formed in early frames, which is different from the WST case in which propagation occurred during the whole experiment.

Conclusion

A procedure to locally refine meshes for Digital Image Correlation (DIC) analyses was proposed. It is referred to as adaptive meshing as it is guided by auxiliary fields, i.e., gray level residuals (artificial case) or maximum eigen strain (experimental cases). At each refinement iteration, the selected elements were divided into four new ones by splitting the edges into two new ones, which created hanging nodes (i.e., nodes with a lack of connectivity). Lagrange multipliers reconnected these nodes to the mesh, ensuring that the displacement of each hanging node be the mean of its parent nodes. This methodology was developed to improve the meshing process for DIC analyses of cracked materials with efficient computational cost, and with at least similar accuracy compared to fine uniform meshes. Although the discussion mainly focused on cracks, other auxiliary fields related to localized phenomena (e.g., phase transition, plasticity) could be used.

The artificial test case illustrated the relationship between the number of cracks and the gain in applying AM procedures. The refinement criterion was based on the gray levels at the pixel scale. It was demonstrated that the gains of AM were inversely proportional to the surface crack density (i.e., higher for single crack analysis and decayed as the number and crack extensions increased). In other words, the studied phenomenon (i.e., cracks) passed from localized to diffuse states.

In the investigated experimental cases, AM was first applied to a Wedge Splitting Test (WST), since a low crack density was expected and highlighted the benefits of the proposed procedure. Five different Uniform Meshes (UMs) (i.e., structured meshes with uniform element size) were tested and then compared to the AM procedure. It was demonstrated that AMs reduced memory storage by decreasing the DIC Hessian size while keeping similar descriptions of the crack kinematics.

The second experimental case was devoted to crack network formation due to the hydration of an MgO castable, with a considerably higher crack density. This case was also related to higher uncertainty levels since the experiment was performed inside a climatic chamber with additional artifacts (i.e., glass windows and humidity/temperature of the environment). Even in this challenging case, the AM procedure showed benefits for reducing the memory load with a smaller Hessian matrix. An additional benefit for this case was that the map for the element size after adapting the mesh already showed the overall crack network. Even if the artificial example pointed out that the benefits would be smaller for high crack densities, the experimental cases proved that it was indeed more efficient for smaller crack densities but it could also be efficient for dense crack networks when a very fine description was sought.

The criterion to divide the elements was based on the standard uncertainty of the Mean Crack Opening Displacement (MCOD) per element. This metric was estimated by analyzing the initial frames of each experiment, where no crack initiation was expected, and using uniform meshes with the same element sizes as those that will make part of the adapted meshes. It was shown that the uncertainties were different for both experiments (because of differences in experimental environment), which further validated the choice of uncertainty quantifications. The final threshold was proportional to the MCOD uncertainty to only account for cracks (i.e., a multiplicative factor > 1 was considered).

Different values of were chosen for the two investigated cases, which indicates that may depend on the surface crack density. Additional investigation are needed to confirm such observation.

Although the numerical test case underscored that the methodology would be best suited for single cracks, it was not straightforward for the experimental cases. Even though both cases benefited from an adapted mesh, the crack network had one additional advantage since crack initiations showed the full crack path from early frames on, meaning that the mesh could be adapted for different frames with similar results. For the WST, with one macrocrack and few branches, the crack path changed throughout the experiment and further studies should be conducted to find a suitable temporal mesh adaptation strategy to describe gradual crack propagation (e.g., using developing meshes, damaging elements in the cracked region or using node splitting techniques).

Last, the AM procedure has a high potential of reducing the computational cost of Digital Volume Correlation analysis, where data are usually considerably bigger than DIC, especially for cases with localized phenomena. Implementing a refinement termination based on the gray level residuals, and also, real-time applications of this procedure during the DIC analyses of many images, are other future improvements of this new methodology.

  (a)). If each element is divided into four new elements of the same size, there is no lack of connectivity at any node, as shown in Figure1(b). However, if only the element 1-4-3 is subdivided into four new ones, it will generate one hanging node (i.e., node 6 in Figure1(c) is a hanging node in element 1-2-4), which does not have the correct number of connectivities (see node 6 in Figure1(b)). Nodes 5 and 7 in Figure1(c)

Figure 1 :

 1 Figure 1: Examples of various localized h-refinements. (a) Original mesh. (b) Global h-refinement. (c) Local refinement (non conformal mesh). (d) Local refinement with conformal mesh

Figure 2 :

 2 Figure 2: Flowchart of the AM procedure applied to DIC, using MCOD as division criterion

Figure 3 (Figure 3 :

 33 Figure3(c)). The corresponding damaged areas (i.e., ratio of cracked pixels wrt. total number of pixels) is equal to 1, 3, and 15%, respectively.

Figure 4 :

 4 Figure 4: Reference and first uniform meshes used in the artificial cases. (a) Starting mesh ( 0 =141 px). (b) First iteration ( 1 =70 px). (c) Second iteration ( 2 =35 px)

Figure 5 :

 5 Figure 5: Meshes for the three artificial cases after 7 refinements (Figure 3). (a) Case #1.1 (b) Case #1.2 (c) Case #1.3

Figure 6 :

 6 Figure 6: (a) Patterned specimen and coarse mesh in yellow (mean element size: ≈ 64 px or 3 mm) used for preliminary analyses. The loading parts are depicted in white dashed lines. The applied vertical force and the horizontally transmitted forces ( ≈ 5.7 ) to the sample are shown as blue arrows. (b) Load displacement curve for the studied WST down to 10% of the peak load. The red circles depict image acquisitions and the blue asterisk depicts the frame used in the discussion.

Figure 7 :

 7 Figure 7: Standard uncertainties of the maximum eigen strain (a) and MCOD (b) as functions of the inverse element size -1 (UMs) for the initial images of the WST. The dashed lines depict power law fits according to Equations (24) and (25).

Figure 8 :

 8 Figure 8: Normalized standard deviation of the gray level residuals (GLR) fields in each frame (% of the dynamic range) for UMs with different element sizes. The vertical dashed line depicts frame #30

Figure 10 :

 10 Figure 10: Mean crack opening displacement (MCOD) fields for the frame #30 of the analysis using UMs with different element sizes (1 px ≡ 48.5 µm)

Figure 11 :

 11 Figure 11: Element size field for different initial thresholds ( * = ( )) for the frame #30 of the WST case. (a) = 12, (b) = 6, (c) = 3, and (d) = 1.5

Figure 12 :

 12 Figure 12: MCOD fields for different initial thresholds ( * = ( )) for the frame #30 of the WST case. (a) = 12, (b) = 6, (c) = 3, and (d) = 1.5

Figure 13 :

 13 Figure 13: Gray level residual field for different initial thresholds ( * = ( )) for the frame # 30 of the WST case. (a) = 12, and the difference between two cases to depict the GLR caused by the refinement (b) 12 -6, (c) 12 -3, and (d) 12 -1.5

5 Figure 14 :

 514 Figure 14: Comparison of the standard deviation of normalized gray level residuals field for the different discretizations (i.e., UM -solid lines, and AM -dash dotted lines). The vertical dashed black line marks frame #30

Figure 15 :

 15 Figure 15: Results using a mesh adapted for the last frame. (a) Element size map. MCOD fields for the 30-th (b) and the last (c) frames

30 Figure 16 :

 3016 Figure 16: Ratio between the GLR for the AM for frame #30-th and for the last one

Figure 17 :

 17 Figure 17: Schematic drawing of the climatic chamber with the refractory castable cube in the center (a), adapted from [66]. Reference image (b) and frame #38 (c) for the MgO hydration experiment. The image sizes are ≈ 67 × 67 mm 2 (≈ 1340 × 1340 px) and the ROI (yellow square) size is ≈ 60 × 60 mm 2 (≈ 1200 × 1200 px)

Figure 18 :

 18 Figure 18: Standard uncertainties of the maximum eigen strain (a) and MCOD (b) as functions of the inverse of the element size for the initial images of the curing and drying experiment. The dashed lines depict power law fits according to Equations (24) and[START_REF] Fagerholt | Investigation of fracture in small-scale SENT tests of a welded X80 pipeline steel using Digital Image Correlation with node splitting[END_REF]. The black stars mark the uncertainties for the optimal AM ( = 1.5)

Figure 19 :

 19 Figure 19: Normalized standard deviation of the gray level residuals (GLR) in each frame (as % of the dynamic range) for uniform meshes with different element sizes. The vertical dashed line depicts frame #38

Figure 20 :

 20 Figure 20: Normalized gray level residual fields for frame #38 using UMs with different element sizes. For the sake of visualization the results from finer meshes were subtracted from the coarser mesh result (64 px). The color bar range was reduced to highlight the patterns in the fields

Figure 21 :Figure 22 :

 2122 Figure 21: Mean crack opening displacement (MCOD) fields for frame #38 when UMs with different element sizes are used (1 px ≡ 50 µm)

Figure 23 :

 23 Figure 23: MCOD fields for picture #38 and different initial thresholds ( * = ( )) of the MgO hydration case. (a) = 12, (b) = 6, (c) = 3, and (d) = 1.5

Figure 24 :

 24 Figure 24: Gray level residual field for different initial thresholds ( * = ( )) for the frame # 38 of the WST case. (a) = 12, and difference between two of them to depict the GLR caused by the refinement (b) 12 -6, (c) 12 -3, and (d) 12 -1.5

5 Figure 25 :

 525 Figure 25: Standard deviation of GLR for the different discretizations (i.e., UM: solid lines, and AM: dash dotted lines). The vertical dashed line marks frame #38

Figure 26 :

 26 Figure 26: Results using a mesh adapted for the last frame. (a) Element size map. MCOD fields for the 38-th (b) and the last (c) frames

38 Figure 27 :

 3827 Figure 27: Ratio between the GLR for the AM for frame #38-th and for the last one

Table 1

 1 Size of Hessian matrices as a function of smallest element sizes it for the artificial cases iteration it [px] reference Case #1.1 Case #1.2 Case #1.3

	1	70	27	27	27	27
	2	35	75	63	72	75
	3	18	243	147	204	243
	4	8.8	867	303	558	882
	5	4.4	3,267	639	1,416	3,552
	6	2.2	12,675	1,335	3,219	11,661
	7	1.1	49,923	2,856	7,278	31,008

Table 2

 2 Percent gains on size of the Hessian matrices in comparison to the reference case.

	iteration	it [px] Case #1.1 Case #1.2 Case #1.3
	1	70	0	0	0
	2	35	16	4	0
	3	18	40	16	0
	4	8.8	65	36	-2
	5	4.4	80	57	-9
	6	2.2	89	75	8
	7	1.1	94	85	38

Table 3

 3

	DIC analysis parameters
	DIC software	Correli 3.0 [49]
	Image filtering	none
	Shape functions	linear (T3)
	Matching criterion	Regularized sum of squared differences
	Interpolant	cubic
	Strain calculation	derivative of shape functions
	Regularization length 20 px

Table 4

 4 Sizes of the Hessian matrices for uniform (UM) and adapted (AM) meshes, computation time, average iterations per image pair, standard deviation of normalized gray level residuals for the series of 126 frames analyzed, and gain in AM Hessian size in comparison to UM with the same smallest element size. The column Feature contains the average element size ( )

	and number of 4-px elements in brackets for AMs
	Type	Feature	Hessian size CT [s]	Average iteration # per image pair	STD of GLR [%] Gain [%]
	UM	= 64 px	672	0.86	6	4.96	-
	UM	= 32 px	2,511	2.79	6	695	-
	UM	= 16 px	9,360	11.5	6	3.37	-
	UM	= 8 px	37,440	58	6	3.08	-
	UM	= 4 px	149,760	445.4	8	2.95	-
	AM ( = 12)	= 12 px [2,592]	10,002	108	7	4.04	93
	AM ( = 6)	= 11 px [3,124]	13,230	131	6	3.52	91
	AM ( = )	= 10 px [3,952]	23,817	188.6	5	3.34	84
	AM ( = 1.5)	= 8 px [6,348]	37,668	371.5	5	3.16	75

Table 5

 5 Sizes of the Hessian matrices for uniform (UM) and adapted (AM) meshes, computation time, average iterations per image pair, standard deviation of normalized gray level residuals for the series of 200 frames analyzed , and gain in AM Hessian size in comparison to UM with the same smallest element size. The column Feature contains the average element size ( )

	and number of 4-px elements in brackets for AMs
	Type	Feature	Hessian size CT [s]	Average iteration # per image pair	STD of GLR [%] Gain [%]
	UM	= 64 px	588	0.9	4.6	3.85	-
	UM	= 32 px	2,028	5.18	4.7	3.10	-
	UM	= 16 px	7,500	15.3	4.8	2.58	-
	UM	= 8 px	30,000	64.6	5.0	2.26	-
	UM	= 4 px	120,000	1677	23.5	2.11	-
	AM ( = 12)	= 34 px [0]	3,078	12.72	4.6	3.34	59
	AM ( = 6)	= 21 px [0]	7,803	34	4.6	2.83	74
	AM ( = 3)	= 14 px [0]	16,845	87.6	4.7	2.53	44
	AM ( = . )	= 10 px [312]	31,422	356.6	4.7	2.29	74

=1.1 px in the present case). Figure4shows the starting mesh along with the first two iterations of the reference mesh, i.e., the case where all elements were subdivided.

Acknowledgments

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001, and the support by the PDSE grant #8881.188511/2018-01 used during the internship of VS at LMT, CAPES (Brazil) and #2010/20920-9, #2018/02801-4, #2018/15266-0 and #2018/23081-0, São Paulo Research Foundation (FAPESP). The authors would like to thank D.B.D. Stafuzza for his help in performing the WST experiment.

Credit authorship statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


A Olmos
email: aolmos@ivia.es

N Boonham

T Candresse

P Gentit

B Giovani

D Kutnjak

L Liefting

H J Maree

A Minafra

A Moreira

M K Nakhla

F Petter

M Ravnikar

B Rodoni

J W Roenhorst

M Rott

A B Ruiz-Garc Ia

J Santala

G Stancanelli

R Van Der Vlugt

C Varveri

M Westenberg

T Wetzel

H Ziebell

S Massart

High-throughput sequencing technologies for plant pest diagnosis: challenges and opportunities

High-throughput sequencing (HTS) technologies have revolutionized plant pest research and are now raising interest for plant pest diagnostics, with plant virus diagnostics at the forefront of development. However, the application of HTS in plant pest diagnostics raises important challenges that plant health regulators will have to address. Adapted infrastructures, technical guidelines and training are pivotal for further use and adoption of the HTS technologies in the phytosanitary framework.

Introduction

Early diagnosis and a rapid response are crucial to reduce the risk of entry and spread of plant pests into an area. In 2016, the Commission on Phytosanitary Measures adopted a recommendation recognizing that 'pest diagnosis is a cross-cutting issue that underpins most International Plant Protection Convention (IPPC) activities. In order to take action against a pest, it must be accurately identified. To enable safe trade, pest diagnosis must further be completed quickly and to a high level of confidence' (FAO, 2016). National Plant Protection Organizations (NPPOs) routinely perform pest diagnosis to support export certification, import inspections, pest surveillance and eradication programmes. Plant pests can be managed most effectively when control measures are implemented at an early stage of infestation. Plant pest diagnostics is based on the use of a range of methods underpinned by different biological principles (e.g. bioassay, biochemical, isolation/ extraction methods, molecular methods, morphological and morphometric and serological methods), some of them being highly specific and others more generic. The ability to detect plant pests varies with the sensitivity and specificity of the detection tools used. The recent development of high-throughput sequencing (HTS) technologies, also called next-generation sequencing (NGS) or deep sequencing, has revolutionized the research on plant-associated organisms. These techniques are beginning to be used for pest diagnostics, in particular for viruses. In this paper, the opportunities offered by these new technologies for the diagnosis of plant pests are discussed. The challenges they will pose and the issues that need to be addressed for appropriate use of HTS within phytosanitary frameworks are highlighted, with an emphasis on plant viruses for which the adoption of HTS technology is more advanced.

What is HTS?

HTS technologies enable the simultaneous sequencing, and thus detection, of any organism present in a sample (see Fig. 1). These technologies are broad-spectrum and generic in nature and can potentially be used to replace a wide range of conventional methods, especially in the screening stage. For known pests, including regulated ones, such unbiased approaches allow the detection of very distant isolates or strains of the same pest. Moreover, HTS technologies have the ability to reveal in any plant sample the presence of unexpected as well as unknown organisms that might be potential threats.

HTS technologies are widely applied in research, accelerating the discovery of new potential pests from symptomatic (resolving the aetiology of disease) and asymptomatic plants and substrates (soil, water, wood, etc.). In the field of plant virology, HTS is revolutionizing scientists' understanding of the ubiquitous nature of viruses and clearly demonstrates that we are just beginning to touch on a large and previously uncharacterized viral diversity [START_REF] Simmonds | Consensus statement: virus taxonomy in the age of metagenomics[END_REF][START_REF] Zheng | VirusDetect: an automated pipeline for efficient virus discovery using deep sequencing of small RNAs[END_REF][START_REF] Yutin | Discovery of an expansive bacteriophage family that includes the most abundant viruses from the human gut[END_REF]. However, these discoveries are often only based on the identification of a partial genomic sequence and the biology of such new species is most often very poorly characterized or not characterized at all. Thus, the increased use of HTS technologies for plant pest diagnosis will accelerate the pace of discovery of new plant pests for which challenging and quick decisions will have to be taken by NPPOs on the basis of very limited information and imprecisely evaluated potential phytosanitary risks.

Opportunities for the use of HTS in pest diagnosis

HTS technologies open new possibilities and opportunities in routine diagnostics for (a) understanding the status of a pest in a region through surveillance programmes, (b) certifying nuclear stock and plant propagation material, (c) (post-entry) quarantine testing to prevent the introduction of pests into a country or area, and (d) monitoring of imported commodities for new potential risks. HTS offers important benefits for each of these applications [START_REF] Rwahnih | Comparison of next-generation sequencing versus biological indexing for the optimal detection of viral pathogens in grapevine[END_REF][START_REF] Hadidi | Next-generation sequencing and genome editing in plant virology[END_REF][START_REF] Rott | Application of next generation sequencing for diagnostic testing of tree fruit viruses and viroids[END_REF].

Knowledge of pest status in a territory is key to phytosanitary regulation and justification of phytosanitary measures. HTS technologies are powerful tools for both targeted and non-targeted analysis, and their use in surveillance programmes has a lot of potential. Pre-designed automated bioinformatics workflows will allow an efficient analysis of datasets, providing information on the presence or absence of all targeted pest(s) in the selected crop(s). Moreover, national large-scale non-targeted analysis of viruses, for example virome scanning for viruses present in cultivated and wild plants of a botanical family or genus, can be performed with these technologies. Countries that know the potential pests already present in their territory will have an advantage in managing the associated phytosanitary risks.

The certification of nuclear stocks and propagation materials also requires screening of specific targeted pests. As such, HTS analysis can be applied in a similar way as described for surveillance programmes. HTS analyses have the potential to eventually replace lengthy and tedious resource-intensive biological indexing for regulated pests, shortening and streamlining indexing cycles; increasing the pace of flow of novel varieties to growers is an additional advantage. In this respect, in several field applications, HTS has been found to be superior to standard bioassays for the detection of low-titre viruses [START_REF] Rwahnih | Comparison of next-generation sequencing versus biological indexing for the optimal detection of viral pathogens in grapevine[END_REF].

Post-entry quarantine testing, as required, for example, for imports of potato germplasm, prescribes specific testing for a number of pests and a broad screen for other pests by visual inspection and bioassays. For screening of specific pests, targeted HTS analysis can be applied. However, HTS will also allow a non-targeted analysis of the generated datasets expanding the post-entry testing to unexpected and unknown pests present in imported plant material.

Other advantages of HTS technologies are related to an increase in the amount of genomic data obtained compared with conventional molecular techniques. This can improve researchers' understanding of the genomic diversity within a species (or at a lower taxonomic level), resulting in a better taxonomic assignment for unambiguously identified pests, a better appraisal of the impact of a pest's genetic diversity and better design of targeted molecular tests.

HTS offers the possibility to store the generated sequencing libraries and datasets for future analyses. As databases expand with information on new pests and bioinformatics tools improve, historical datasets will become available for re-analysis, enabling NPPOs to perform post hoc bioinformatic surveys. This is a clear benefit compared with conventional molecular methods which require new laboratory analyses on properly stored original samples or isolated nucleic acids.

Challenges of the application of HTS within phytosanitary frameworks

Adopting HTS technologies for routine diagnostic purposes raises challenges that may have significant implications within a phytosanitary framework.

Laboratory and bioinformatics challenges

As for conventional methods, interpretation of the test results derived from sequence analysis requires a high level of expertise in plant pathology and taxonomy. For laboratories, validation and routine application of HTS technologies demand new bioinformatics skills not previously required in a routine testing laboratory and investment in infrastructure to enable rapid analysis and data storage. Whilst similar principles of validation are required (analytical sensitivity, analytical specificity, repeatability and reproducibility), HTS experiments produce data on a large scale (essentially data for many tests performed all at the same time), making analysis more complex and necessitating novel skills in manipulation and storage of the large numbers of sequence data thus generated.

Interpretation of results, biological significance of new findings and pest risk analyses

The use of HTS-based diagnosis by regulators will be a major challenge due to the increased number of new pests being discovered for which a biological context may be missing. For plant viruses, [START_REF] Massart | A framework for the evaluation of biosecurity, commercial, regulatory and scientific impacts of plant viruses and viroids identified by NGS technologies[END_REF] proposed a framework within which relevant information can be generated for new findings. This framework will be useful for any novel organism detected by HTS, describing the characterization that should take place when starting from the sole genomic sequence. It is essential that when making regulatory decisions consideration continues to be given to the potential risks associated with the discovery of a new organism.

The complexity of disease aetiology in field conditions may pose additional challenges, for example in the case of mixed viral infections. In fact, plants are commonly infected with more than one virus species, especially in woody, vegetatively propagated plants such as grapevine [START_REF] Jooste | Identification and distribution of multiple virus infections in Grapevine leafroll diseased vineyards[END_REF]. A possible synergism between coinfecting viral species might lead to unpredictable variations in symptoms, infectivity, accumulation and/or vector transmissibility. In fact, previously undescribed viruses not associated with deleterious symptoms found in some hosts may provide benefits to the host and may have been positively selected for during cultivation [START_REF] Roossinck | Move over bacteria! Viruses make their mark as mutualistic symbionts[END_REF].

The characterization of the biological properties of any organism identified through HTS technologies will need to be performed within a short timeframe in order to provide the basis for assessing the risk it poses. An additional challenge is the allocation of appropriate funding to perform these biological characterizations and risk assessments [START_REF] Massart | A framework for the evaluation of biosecurity, commercial, regulatory and scientific impacts of plant viruses and viroids identified by NGS technologies[END_REF]. Moreover, it means that pest risk analysis (PRA) processesfrom initiation, through to the pest risk assessment and finally the pest risk managementwill need to follow the same pace as discovery of new potential pests via HTS technologies.

Data sharing

Whilst it is standard practice among academics to share HTS datasets generated for research purposes, the sharing of such 'raw' datasets is less likely in a regulatory context. There could be benefits in sharing these datasets to facilitate movement of plants (e.g. by sharing information on the health status of propagative material). However, it is recognized that release of such data may pose a risk of the misuse of new findings to hamper the international movement of plants and plant products. Therefore, HTS data management and accessibility is another important decision to be taken by plant health regulators, for which a cost-benefit evaluation will be required. There may be additional complexities for regulators with datasets generated for research purposes that concern endemic pests or even healthy hosts explored by HTS technologies. These datasets may be publicly accessible and inadvertently reveal the presence of regulated pests in previously unknown locations or hosts.

What are the main issues to address for plant health regulators in the HTS era?

Implementing HTS technologies for plant pest diagnosis will raise some questions that need to be addressed:

Technical guidelines for routine analysis by laboratories

HTS technologies, as is the case for other techniques, face similar challenges with regard to proper sampling procedures, nucleic acid template preparation, appropriate controls, contamination and validation issues. HTS also requires confirmation of the results by additional techniques. Contamination is more problematic than with other methods, especially when high-titre infections are encountered and despite additional measures being taken (e.g. rolling programmes of tag use, increased instrument decontamination). Well-defined controls/reference materials and confirmatory testing will solve this problem and help to unambiguously identify infested samples. Thus, technical and bioinformatics guidelines for the use of HTS in plant pest diagnostics are needed to help address these issues. Guidelines will also be needed for reporting the results and for documenting and safeguarding the HTS data and the results of the bioinformatic analyses.

Infrastructure to analyse, host and share HTS data

Implementation of HTS in plant health will require significant investments in information technology and bioinformatics. Large files up to a few gigabytes per sample are generated and need to be stored and properly backed-up on servers for long periods of time. Transferring files from servers to data analysis computers requires a fast network, with machines that need high computational power to handle large data files.

Education/training to adapt to the new paradigm

Expertise in taxonomy of pests, classical diagnostics and epidemiology will continue to be essential. Additionally, analysts will require a more in-depth knowledge on bioinformatics for the appropriate interpretation of the data. Risk assessors, risk managers and plant health regulators should be trained on the background of HTS-based diagnosis to enable them to make fast and appropriate decisions.

A framework for biological characterization of micro-organisms to support decision-making by the regulators

A framework for the evaluation of biosecurity, commercial, regulatory and scientific impacts of plant viruses and viroids identified by HTS technologies has been proposed by a European consortium of researchers [START_REF] Massart | A framework for the evaluation of biosecurity, commercial, regulatory and scientific impacts of plant viruses and viroids identified by NGS technologies[END_REF]. It provides useful guidance for strategies for generating the biological and other relevant data needed for PRA. Appropriate communication channels between regulatory authorities (such as NPPOs) and diagnostic laboratories need to be enhanced.

Reinforcing baseline surveys

The application of HTS technologies in research and their future large-scale application in pest diagnostics will identify poorly characterized or new potential pests. Ultimately for a NPPO, the key question will be whether this novel finding from HTS is truly a new incursion rather than a pest that has been present in an area or imported commodity but previously unreported due to limitations in existing diagnostic technologies. This will put a far greater emphasis on baseline survey activity to demonstrate what pests are currently present in an area/imported commodity. For example, some countries (such as Belgium) have already launched national-based surveys using HTS to detect any virus present in cultivated and wild plants of selected botanical families or genera within their respective territories.

To summarize

In summary, HTS technologies are powerful tools that are rapidly evolving and open up unprecedented possibilities. By overcoming the current challenges, they have the potential to transform the field of plant diagnostics by enabling quick, efficient and broad-spectrum testing using generic procedures applied in a range of diverse disciplines, thus providing a major contribution to the progress of the management of plant pest control programmes.

Technologies de s equenc ßage g en etique a haut d ebit pour le diagnostic des organismes nuisibles aux v eg etaux: d efis et perspectives Les technologies de s equenc ßage g en etique a haut d ebit ont r evolutionn e la recherche des organismes nuisibles aux v eg etaux et suscitent maintenant de l'int erêt dans le contexte de leur diagnostic, avec au premier plan celui des virus des v eg etaux. Cependant, l'application de ces technologies au diagnostic des organismes nuisibles aux v eg etaux fait emerger de nombreux challenges auxquels les personnes d eveloppant la r eglementation en sant e des v eg etaux devront faire face. Des infrastructures adapt ees, des lignes directrices techniques ainsi que la formation sont des el ements clefs pour l'utilisation future et l'adoption de ces technologies dans le contexte phytosanitaire.

Fig. 1

 1 Fig. 1 High throughput sequencing technologies and traditional detection and identification methods in plant pest diagnosis. Examples shown in this figure are for viruses, fungi and bacteria but HTS can be used other pests such as arthropods and nematodes.

  Пpимeнeниe тexнoлoгий выcoкoпpoизвoдитeльнoгo ceквeниpoвaния для диaгнocтики вpeдныx для pacтeний opгaнизмoв: cyщecтвyющиe пpoблeмы и oткpывaющиecя вoзмoжнocти Texнoлoгии пиpoceквeниpoвaния (HTS) пoдняли нa кaчecтвeннo нoвый ypoвeнь иccлeдoвaния вpeдныx для pacтeний opгaнизмoв, ивн acтoящee вpeмя oни пoвышaют интepec кв oпpocaмд и aгнocтики вpeдныx opгaнизмoв. Пpиэ т oмд и aгнocтикa виpycoв pacтeний явлeтcя oднoйи зп pиopитeтныx зaдaч.O днaкo пpимeнeниe HTS для диaгнocтики вpeдныx для pacтeний opгaнизмoв coпpяжeнocцeлым pядoмв aжныx пpoблeм, нa кoтopыe дoлжнo быть oбpaщeнo внимaниe нopмaтивныxo p гaнoв фитocaнитapии.B o тп oчeмy для дaльнeйшeгo иcпoльзoвaния и пpинятия нa вoopyжeниe тexнoлoгий пиpoceквeниpoвaния (HTS) в фитocaнитapныxc тpyктypax нeoбxoдимa иxa дaптaция, coглacoвaнныep y кoвoдящиe пpинципы и тpeнинг.

Acknowledgements

The authors would like to thank Ms Sophie Petter for preparing the figure.