Mehdi Ennassir

Keywords: Digital volume correlation, Image-based models, X-ray micro-tomography, Finite element method, Fictitious domain methods, Cellular materials Corrélation d'images volumiques, Modèles basés images, Micro-tomographie par rayons X, Méthode des éléments finis, Méthodes

Il est temps pour moi de remercier toutes les personnes qui ont rendu l'achèvement de cette thèse réalisable ainsi que tous ceux et celles que j'ai pu rencontrer à Toulouse pendant ces six dernières années.

Ma gratitude s'adresse tout d'abord à mes encadrants de thèse Jean-Noël Périé, Robin Bouclier et Jean-Charles Passieux. Merci de m'avoir donné l'opportunité de faire cette thèse et de l'avoir encadrée dans les meilleurs conditions. Merci pour votre bienveillance, votre rigueur scientifique et votre écoute. Je vous en serai toujours reconnaissant. J'ai beaucoup apprécié le fait de vous côtoyer aussi hors le labo, tout en haut dans les airs de Luchon ou au bord de la mer à Giens. Je voudrais chaleureusement remercier Pascal Doumalin, sans qui le passage de mes travaux de thèse de la 2D à la 3D n'aurait pas vu le jour. Merci de nous avoir fourni les données liées aux essais in-situ sur la mousse polyuréthane ainsi que pour toutes les discussions que nous avons eues.

Je souhaite également remercier chaleureusement tous les membres de mon jury de thèse, en particulier Johan Hoefnagels et Pierre Kerfriden pour avoir accepté de rapporter mes travaux ainsi que pour leurs remarques et questions. Merci également à Annalisa Buffa, Éric Maire et Julien Réthoré pour avoir également examiné et contribué à la discussion. Merci à Ludovic Barrière et Michel Fournié pour avoir aussi participé à cette dernière. Encore, un grand Merci à Michel pour m'avoir intégré dans l'équipe d'enseignement de l'IUT GMP de Toulouse ainsi que pour tous tes conseils.

Un grand Merci à Guillaume Colantonio, Morgane Chapelier et Raphaël Fouque pour l'ambiance très sympathique qui régnait dans

Abstract

Measuring displacement and strain fields at low observable scales in complex microstructures still remains a challenge in experimental mechanics often because of the combination of low definition images with poor texture at this scale. The problem is particularly acute in the case of cellular materials, when imaged by conventional micro-tomographs, for which complex highly non-linear local phenomena can occur. As the validation of numerical models and the identification of mechanical properties of materials must rely on accurate measurements of displacement and strain fields, the design and implementation of robust and faithful image correlation algorithms must be conducted.

With cellular materials, the use of digital volume correlation (DVC) faces a paradox: in the absence of markings of exploitable texture on/or in the struts or cell walls, the available speckle will be formed by the material architecture itself. This leads to the inability of classical DVC codes to measure kinematics at the cellular and a fortiori sub-cellular scales, precisely because the interpolation basis of the displacement field cannot account for the complexity of the underlying kinematics, especially when bending or buckling of beams or walls occurs.

The objective of the thesis is to develop a DVC technique for the measurement of displacement fields in cellular materials at the scale of their architecture. The proposed solution consists in assisting DVC by a weak elastic regularization using an automatic image-based mechanical model. The proposed method introduces a separation of scales above which DVC is dominant and below which it is assisted by image-based modeling.

First, a numerical investigation and comparison of different techniques for building automatically a geometric and mechanical model from tomographic images is conducted. Two particular methods are considered: the finite element method (FEM) and the finite-cell method (FCM). The FCM is a fictitious domain method that consists in immersing the complex geometry in a high order structured grid and does not require meshing. In this context, various discretization parameters appear delicate to choose. In this work, these parameters are adjusted to obtain (a) the best possible accuracy (bounded by pixelation errors) while (b) ensuring minimal complexity.

Concerning the ability of the mechanical image-based models to regularize DIC, several virtual experimentations are performed in two-dimensions in order to finely analyze the influence of the introduced regularization lengths for different input mechanical behaviors (elastic, elasto-plastic and geometrically non-linear) and in comparison with ground truth. We show that the method can estimate complex local displacement and strain fields with speckle-free low definition images, even in non-linear regimes such as local buckling.

Finally a three-dimensional generalization is performed through the development of a DVC framework. It takes as an input the reconstructed volumes at the different deformation states of the material and constructs automatically the cellular micro-architeture geometry. It considers either an immersed structured B-spline grid of arbitrary order or a finite-element mesh. An experimental evidence is performed by measuring the complex kinematics of a polyurethane foam under compression during an in situ test.

Résumé

La mesure de champs de déplacement et de déformation aux petites échelles dans des microstructures complexes représente encore un défi majeur dans le monde de la mécanique expérimentale. Ceci est en partie dû aux acquisitions d'images et à la pauvreté de la texture à ces échelles. C'est notamment le cas pour les matériaux cellulaires lorsqu'ils sont imagés avec des micro-tomographes conventionnels et qu'ils peuvent être sujets à des mécanismes de déformation complexes. Comme la validation de modèles numériques et l'identification des propriétés mécaniques de matériaux se base sur des mesures précises de déplacements et de déformations, la conception et l'implémentation d'algorithmes robustes et fiables de corrélation d'images semble nécessaire.

Lorsque l'on s'intéresse à l'utilisation de la corrélation d'images volumiques (DVC) pour les matériaux cellulaires, on est confronté à un paradoxe: l'absence de texture à l'échelle du constituant conduit à considérer l'architecture comme marqueur pour la corrélation. Ceci conduit à l'échec des techniques ordinaires de DVC à mesurer des cinématiques aux échelles subcellulaires en lien avec des comportements mécaniques locaux complexes tels que la flexion ou le flambement de travées.

L'objectif de cette thèse est la conception d'une technique de DVC pour la mesure de champs de déplacement dans des matériaux cellulaires à l'échelle de leurs architectures. Cette technique assiste la corrélation d'images par une régularisation élastique faible en utilisant un modèle mécanique généré automatiquement et basé sur les images. La méthode suggérée introduit une séparation d'échelles au dessus desquelles la DVC est dominante et en dessous desquelles elle est assistée par le modèle mécanique basé sur l'image.

Une première étude numérique consistant à comparer différentes techniques de construction de modèles mécaniques basés sur les images est conduite. L'accent est mis sur deux méthodes de calcul particulières: la méthode des éléments finis (FEM) et la méthode des cellules finies (FCM) qui consiste à immerger la géométrie complexe dans une grille régulière de haut ordre sans utiliser de mailleurs. Si la FCM évite une première phase délicate de discrétisation, plusieurs paramètres restent néanmoins délicats à fixer. Dans ce travail, ces paramètres sont ajustés afin d'obtenir (a) la meilleure précision (bornée par les erreurs de pixellisation) tout en (b) assurant une complexité minimale.

Pour l'aspect mesure par corrélation d'images régularisée, plusieurs expérimentations virtuelles à partir de différentes simulations numériques (en élasticité, en plasticité et en non-linéarité géométrique) sont d'abord réalisées afin d'analyser l'influence des paramètres de régularisation introduits. Les erreurs de mesures peuvent dans ce cas être quantifiées à l'aide des solutions de référence éléments finis. La capacité de la méthode à mesurer des cinématiques complexes en absence de texture est démontrée pour des régimes non-linéaires tels que le flambement.

Finalement, le travail proposé est généralisé à la corrélation volumique des différents états de déformation du matériau et à la construction automatique de la micro-architecture cellulaire en utilisant soit une grille B-spline d'ordre arbitraire (FCM) soit un maillage éléments finis (FEM). Une mise en évidence expérimentale de l'efficacité et de la justesse de l'approche proposée est effectuée à travers de la mesure de cinématiques complexes dans une mousse polyuréthane sollicitée en compression lors d'un essai in situ.

General introduction

Cellular materials exist everywhere in nature and have been used for centuries by mankind. Throughout the human history, the technological development was based on the innovative design and use of complex hierarchical cellular materials such as wood [Gibson 2003]. More recently, in multiple industries, complex materials such as foams, honeycombs or wood are widely used to their exceptional specific performances [Xiong et al. 2019,Feng et al. 2020,Miranda et al. 2021]. From an environmental perspective, a large interest is now re-directed towards this kind of materials given their low environmental impact. Characterizing the mechanical behavior of materials with complex microstructures has therefore become necessary to design low material consumption and efficient eco-structures. However, the need is not limited to industry specific solutions [START_REF] Benedetti | [END_REF]. In biology applications for instance, an accurate characterization of the mechanical behavior of bones could contribute to a better design of bone grafting components. Especially, with the increase of interest towards additive manufacturing, it has become possible to design complex biomimetic cellular architectures [Wegst et al. 2015, Wang et al. 2016, Robles-Linares et al. 2019].

Natural materials usually present multi-scale architectures. As an example, the human body bone structure is subdivided into (i) the macro-structure which defines the largest observed scale with the naked eye (skeleton bones ranging from the smallest phalanx bone ≈ 5 mm to the femoral bone ≈ 480 mm), (ii) the micro-structure which contains multiple complex components such as the spongy bone architecture (with bone cells ≈ 0.1 mm) and blood vessels, (iii) the nano-strucutre or ultrastructure which defines what is under the cell-strut scale to finally achieve (iiii) the collagen/crystal molecular structure. The top of Fig. 1 illustrates the different scales of bone. Similarly, this hierarchical structure can also be observed in plants such as bamboo, see the bottom of Fig. 1. This figure clearly illustrates the complex and organized high heterogeneity of natural materials.

Multiple questions can come to our mind when observing the architectures of Fig. 1: can we find models that describe their mechanical behavior ? Does it make sense to try to model blindly a geometry at the finest scale using finite element structural analysis ? We can all agree that the answer to the last question is no as no existing computer can allow to do such computations. Therefore, the main challenge for characterizing architected materials is to find a way of making an accurate dialogue between the different scales and with minimal computational cost. For the particular case of the hierarchical bamboo architecture, a modeling approach based on homogenization along with scanning electron microscopy (SEM), allowed to predict the stiffness and strengh based on the microstructure [Gangwar & Schillinger 2019].

In this general context of the charaterization of materials, computed micro-tomography (µ-CT) has become one of the most efficient non-invasive imaging tools for representing and inspecting the physical reality of cellular materials [Baruchel et al. 2000, Maire et al. 2003]. It currently makes it possible to reveal the internal architecture of certain materials at a micrometric scale [Amani et al. 2018], or even information on the microstructure of metallic materials [START_REF] Ludwig | [END_REF], Plumb et al. 2018]. By using in situ testing machines [Buffiere et al. 2010], it is possible to assess the effects of loading on internal deformation at various scales [START_REF] Gustafsson | [END_REF], Amani et al. 2020, Wu et al. 2022] or damage [Amani et al. 2018]. In this context, digital volume correlation (DVC) is now commonly used to obtain a three-dimensional displacement field from a sequence of absorption contrast tomographic Figure 1: The fascinating complex hierarchical structure of bone and bamboo. Image taken from [Wegst et al. 2015].

images [START_REF] Bay | [END_REF]. DVC consists in essence in finding the unkown transformation that performs the best possible correlation (or registration) of the images recorded at rest and once the specimen is loaded. DVC has been successfully applied to a wide range of natural and artificial materials [START_REF] Bay | [END_REF], Verhulp et al. 2004, Lenoir et al. 2007, Forsberg et al. 2010] presenting, once imaged, favorable textures at the macro-scale.

This work focuses on open-cell foam like materials with the strong hypothesis of considering only two scales: the macro scale (i.e. volumes containing multiple cells) and the micro-scale which focuses on each individual cell strut. One of the long term aims of this work is to develop procedures to identify the local mechanical properties of cellular materials using a rigorous dialogue between full field measurements and numerical simulations based on µ-CT images. This kind of dialogue can Contents be practically inscribed in the Finite Element Method Updating (FEMU) class of algorithms [Avril et al. 2008,Roux & Hild 2020]. On the one hand, it is necessary to simulate the mechanical response of the material (direct problem). Given the µ-CT images, we can try to build a faithful mechanical imagebased model. On the other hand, such an experimentation/simulation coupling requires minimally invasive measurements that can be obtained using DVC.

However, such dialogues are usually conducted at low spatial resolution and in the case of elastic behaviour [START_REF] Zauel | [END_REF]]. One of the challenges in the field of experimental mechanics is indeed to perform such DVC measurements at the architecture scale [Dall'Ara et al. 2014, Xu 2018], i.e. to see how each individual cell-strut deforms.

As a first step towards this long-term goal, this present contribution focuses on the development of adequat numerical tools that allow to measure accuratly displacement fields below the cellular scale. This method can be applied to a wide spectrum of materials including metallic/polymeric foams, biological tissues (trabecular bones), cellular woods, or additive manufacturing materials such as lattice structures. As an example, an image of a Rohacell-51 polymetacrylimid closed cell foam microstructure obtained using X-ray micro-tomography is given in Fig. 2.

Figure 2: Image of a Rohacell-51 polymetacrylimid closed cell foam microstructure obtained using X-ray micro-tomography. The voxel size is equal to 16 µm and the cell-struts are defined by only 2 to 3 pixels along the thickness direction. The acquisition was peformed using a conventional X-ray tomograph (EasyTom 130) available in Clément Ader Institut.

A two-dimensional analysis (using DIC or Stereo-DIC) is possible by artificially adding a high frequency speckle pattern to the observed surfaces [Maraghechi et al. 2020]. Numerous techniques exist that allow textures to be deposited over a wide range of scales. However, in volume analysis (DVC), depending on the imaging modality, the variations in grey-levels that generate a DIC suitable texture are associated with the micro architecture and/or the heterogeneity of the constituents. For instance in Fig. 2, the acquisition parameters and the size of the sample were such that the resolution of 16 µm per voxel allows for only 3 voxels on average in the strut thickness. We can see that the struts are not textured at all. Anyway, with such a resolution, one would not even be able to see a sub-cellular speckle, even if it exists (absorption contrast provided by the architecture). With such microstructures, we are confronted with a paradox: the scale of the constituents is merged with that of the texture, whereas the texture should be defined at a lower scale. This problem has led DVC users to consider elements (global DVC) or subsets (local DVC) of very large size compared to the micro architecture [START_REF] Bay | [END_REF],Gomes Perini et al. 2014,Xu 2018]. The strain fields obtained with such choices are therefore associated with a meso (or macro) scale which is homogenized with respect to the material architecture scale. The lack of texture at a smaller scale precludes the consideration of smaller elements or subsets, and therefore to access to more resolved measurements. Of course, there have been attempts to deposit texture in volume, especially in manufactured materials using Barium Sulfate [Passieux et al. 2013] or copper particles [Brault et al. 2013] as contrast agent for instance. But, apart from the fact that it is not easy to guarantee a homogeneous and isotropic texture and that it cannot be generalised to any material (especially biological ones), this invasive technique may have effects on the behaviour of the material we want to characterize. This technical barrier which prevents performing strain measurement under the cellular scale represents one of today's most challenging issues in DVC.

This work proposes a method that breaks this barrier and increases the resolution despite the absence of texture at the cell scale. In order to represent the complex cellular kinematics, a voxel (or sub-voxel measurement resolution) is needed. This is equivalent to solving the original optical flow problem [Horn & Schunck 1981]. In order to achieve this goal, the proposed methodology makes use of two main ingredients: (a) a kinematic basis that is adapted to the sub-cellular scale is defined using image-based mechanical modeling and (b) mechanical regularization which consists in adding a smoothness constraint to the DVC measurement. More precisely, this smoothness constraint is given by a penalization term (w.r.t the DVC optimization problem) based on the elastic equilibrium gap in the bulk of the volume and a stabilization term acting on the non-free boundaries [Réthoré Part I

Image-based mechanical modeling

Image-based mechanical modeling consists in generating automatically a geometric representation of an object from image data for numerical simulation. Depending on the studied material, this technique can be relatively complex as it depends on multiple factors imposed by the acquisition process. Resolution, noise and the sample size are usually imposed by the experimental environement and the type of the studied material. For cellular materials, the geometry characterization can be considered relatively simple as we can suppose that the material is constituted from void and a solid, therefore the geometry can be fully characterized by a unique level-set function, which can be viewed as a binary geometrical representation which differentiates between void and material. However, when dealing with more complex materials such as composites or metals, the segmentation step in addition to the assembly of all the geometric parts into one unique object is not immediate. As we are particularly interested in the construction of a mechanical model of materials with complex cellular mircro-architectures (such as foams, lattice structures, bones or cell woods for example), this first part of the thesis presents some techniques that allow to build geometric and mechanical models of cellular materials automatically from grey-level images in addition to our choices and contributions in this subject. In this part, we present our strategy for generating a mechanical image-based model based on two-dimensional fronto-parallel imaging and three-dimensional reconstructed tomographic images. We trust the reconstruction process of the volumes and no investigation on the physics behind imaging is proposed. We rely on two methods: the finite element (FE) approach based on meshing and the finite-cell method (FCM) which is a fictitious domain method based on level-sets. The methodology is general and can be applied to any type of physical problem. The focus on linear elasticity is justified by the fact that the regularization model for the image correlation problem suggested in part II relies on the equilibrium gap using an elastic model. A numerical implementation of both methods in the context of linear elasticity is performed. As the final goal of this thesis is to perform DVC measurements on real cellular architectures, we insist more on the practical aspect and no clear mathematical investigation of the methods is considered. This first part details two major aspects in image-based modeling: the geometry part and the mechanical analysis part and it is organized as follows:

• Chapter 1:

This chapter deals with the extraction of contours from grey-level data. Even though this subject is very large and rich, we only present a restricted approach that applies to cellular architectures with a reasonable level of noise and with a minimal sufficient resolution. All the practical and difficult aspects related to imaging artefacts are omitted. We present how we can practically construct a continuous and smooth representation from voxels. This has two main interests for this work. First the extraction of regular contours (implicit surfaces and closed polyhedric surfaces) and the sub-voxel evaluation of grey-levels which is a fundamental aspect in image correlation (as it will be underlined in Part II). This chapter details the practical implementation of a unified image representation framework based on B-splines that is higly inspired from the works by [Unser 1999]. The focus is set on finding the "best" image representation that allows to extract smooth and regular boundaries so that geometric errors are reduced in the mechanical analysis and correlation step.

• Chapter 2:

We present some image-based mechanical models that exist in the litterature starting from finite element methods (FEM) and finishing with immersed methods (more precisely the finitecell method (FCM)). We show different two and three dimensional examples in addition to a verification of the classical numerical convergence evidence for the FCM method on an exact analytical problem. Finally, a numerical development is adapted for mechanical image-based modeling. As the FCM method induces more parameters than the classical FEM method, their tuning can lead to heavy unnecessary computations especially when the geometric error of the model is high. We present a numerical investigation of the image-based FCM in order to analyze its convergence behavior when geometries are erroneous (i.e. when they come from an imaging reconstruction process).

Chapter 1

Boundary extraction from grey-level images

The goal of this chapter is to construct a continuous representation of images for both smooth boundary extraction and efficient sub-voxel grey-level evaluation in DIC/DVC.

Before proceeding to that, let us recall what is a level-set and how it is determined. Boundary or level-set extraction from grey-level images is a vast and general subject in image-processing. It can also be named as contour detection and is classified as image segmentation problems. It consists in finding the contours of elements of different shapes and topologies within images. We emphasize again that in the context of this thesis, only cellular architectures with two phases (void and solid) are studied, therefore the domain of interest is fully characterized by only one signed distance field.

In this case, the signed distance scalar function denoted φ and that completely represents the domain Ω is defined as follows:      φ(x) < 0 x ∈ int(Ω) (interior: solid domain) φ(x) = 0 x ∈ ∂Ω (boundary) φ(x) > 0 elsewhere (void domain)

.

(1.1)

Numerically, the level-set function is discretized (or sampled) with a finite number of degrees of freedom. Therefore, from a discrete point of view, a level-set function can be viewed as a matrix of distances. A large class of level-sets named evolution level-sets propagate an initial curve using the convection-diffusion equation in order to fit to regions in the image [START_REF] Osher | [END_REF], Chan & Vese 2001,Bernard et al. 2009]. A second manner of identifying regular contours of images is to find a smooth physical representation with minimal oscillations and that conserves the local variations of the target images. A regular contour can therefore be obtained by taking an iso-value of this continuous representation. This approach is easy to implement but again is not an explicit parametrization of the geometry (by the level-set definition). Once again, these smooth image representations can also be used for efficient sub-voxel evaluation for image correlation (see chapter 3). In the following, we denote the studied image I which is a piece-wise constant function from R d to R. In the same time, we represent it with the tensor that represents the reconstructed voxel grid (I ∈ R n i ×n j ×n k in three-dimensions).

1.1 Evolution level-sets based on the transport equation, a brief state of the art

In two dimensions, the idea of evolution level-sets is to make evolve an initial curve or contour according to a transport equation. Assuming that the temporal variation of the contour curve is guided by the normal speed, it was shown in [START_REF] Osher | [END_REF]] that the level-set function φ is the solution of a convection diffusion equation with an initial condition (φ(t = 0, x) = φ 0 (x)). φ 0 is the initial level set function. Depending on the chosen speed, a large class of evolution equations can be derived (for example flame propagation models or mean-curvature flows). The propagation speed contains a term of edge detection or attraction given by a function that vanishes when the gradient is very large. In addition to this function, the curvature of the level-set can also be added in order to add a smoothness constraint. The level-set method is general and has been applied to different problems such as two-phase fluid flows [Sethian & Smereka 2003], structural optimization [Sethian & Wiegmann 2000], crack propagation in mechanics [START_REF] Moës | [END_REF] and the representation of complex microstructures [Legrain et al. 2011]. In image processing, it is known as active contour models [Chan & Vese 2001] (Chan-Vese algorithm). For this algorithm, the image contours are defined as a level-set function which minimizes an energy functional that takes into account different fitting parameters.

The Chan-Vese model mimizes the sum of four criterias: the total length of the contour, the total area bounded, the grey-level dispersion in the first phase (foreground) and the second phase (background).

When differentiating the Chan-Vese function with respect to the level set φ, a transport equation is obtained as it was defined in [START_REF] Osher | [END_REF]] and can be solved numerically. The active contours approach is usually used to target a specific region in an image. The algorithm usually starts with an initial spline (set by a user) and makes it evolve in order to fit to the desired boundary (like a snake, that is why some variants of these models are also called snake models). When using these methods, the user usually observes the image and fixes a number of initial closed curves (typically circles) and launches the algorithm so that each curve performs its segmentation. This approach seems to be not suitable for our applications where the volumes of cellular architectures contain an unkown and large number of cells all connected. We however find that these methods can be applied in cases where the classical binary segmentation fails (see [Getreuer 2012] for more details about the Chan-Vese algorithm and some two-dimensional examples). In the field of Digital Image Correlation, other techniques based on the correlation between the image to segment and virtual images had also been proposed in [START_REF] Réthoré | [END_REF].

Binary thresholding based on smooth representations of images

The second manner of identifying regular contours is to find a smooth physical representation of the images with minimal oscillations and that conserves the local variations of the target images.

To construct the level-set function, we consider a binary thresholding of this smooth representation denoted Ĩ. By introducing the grey-level thresholding value γ, the domain can be entirely described by φ(x) = γ -Ĩ(x) ∀x ∈ Ω.

(1.2)

This section details how we can efficiently construct the image Ĩ using essentially B-spline functions.

1.2. Binary thresholding based on smooth representations of images

B-spline representation of images

One of the tools for representing images smoothly are splines. Splines are piece-wise polynomial functions with a chosen condition of regularity on the boundaries of each interval. Given a set of data values, multiple approaches can be used for fitting the values using either interpolation splines, least squares fitting, etc. Generally, a continuous function Ĩ can be viewed as a linear combination of basis functions such as B-splines (or any other continuous basis) as it will be detailed below.

Definition 1.2.1 (The spline space). A function s belongs to the space of splines of degree p (denoted S p ) if it satisfies the following two conditions [Schoenberg 1973]:

• s is p-1 times differentiable and all its derivatives are continuous to the order p-1 (s ∈ C p-1 (R)).

• s is a piecewise polynomial function on each knot span

[k, k + 1[ with k ∈ Z (s ∈ P p ([k, k + 1[)).
The indices k are called spans in the spline terminology.

Definition 1.2.2 (B-spline of degree p). A B-spline function of degree p is defined as the following recursive convolution

β p + = β p-1 + * β 0 + with β 0 + = 1 [0,1] (1.3)
where * is the convolution operation and

1 [0,1] is the indicator function of [0, 1].
B-splines are symmetric, positive and have a compact support of length p + 1. A degree p basis function can be seen as a convolution or a regularization of the step function. Eq. (1.4) shows the direct analytic formula of quadratic and cubic B-splines. For p = 1, B-spline basis functions are equal to the linear Lagrange basis functions.

β 2 + (x) =                        0 x ≤ 0 1 2 x 2 x ∈ [0, 1] 1 2 + x -1 -(x -1) 2 x ∈ [1, 2] 1 2 -(x -2) + 1 2 (x -2) 2 x ∈ [2, 3] 0 x ≥ 3 , β 3 + (x) =                                0 x ≤ 0 1 6 x 3 x ∈ [0, 1] - 1 2 x 3 + 2x 2 -2x + 2 3 x ∈ [1, 2] 1 2 x 3 -4x 2 + 10x - 22 3 x ∈ [2, 3] - 1 6 x 3 + 2x 2 -8x + 32 3 x ∈ [3, 4] 0 x ≥ 4
.

(1.4)

Fig. 1.1 displays the B-spline basis functions for p ∈ {0, 1, 2, 3}.

Figure 1.1: B-spline basis functions from p = 0 to p = 3.

Binary thresholding based on smooth representations of images 13

Definition 1.2.3 (Shifted B-spline). The symmetric B-spline β p is defined as the translated B-spline

β p + : β p (x) = β p + x + p + 1 2 (1.5)
Theorem 1.2.1 (Representation of splines). Every spline s ∈ S p has a unique and stable representation in terms of its B-spline expansion [START_REF] Schoenberg | [END_REF], Unser 1999]:

s(x) = k∈Z c k β p (x -k) (1.6)
Remark 1.2.1 (Properties). The B-spline functions (β(xk)) k∈Z represent a basis for the space of polynomial splines. They are linearly independant and all sum to 1 (partition of unity property).

When considering a unit spacing and an infinite number of knots, the elements of the spline space are called cardinal splines and are well suited for image processing [Unser et al. 1991, Unser 1999]. A spline function s can be viewed as a continuous representation of a signal of coefficients (c k ) k∈Z . This means that by making the hypothesis that a signal is piece-wise polyomial implies by theorem 1.2.1 that there exists coefficients c k , the coefficients of the B-spline basis, that allow to approximate it. Given a discrete signal I, a spline representation of I consists in finding the optimal coefficients c k (optimal relative to a problem dependant criteria) that verify:

Ĩ(x) = k∈Z c k β p (x -k).
(1.7)

From now on, let us consider the application of formula (1.7) in order to construct the continuous representation of a voxel grid I ∈ R n i ×n j ×n k (the definition is the same for one and two-dimensional images). If the origin of the image is denoted by (x min , y min , z min ) and the voxel size is (∆x, ∆y, ∆z) then we define the continuous representation of I as follows:

• First, as the image is of finite size and non periodic, the domain of approximation of Ĩ is restricted to:

Ω = [x min + p∆x 2 , x max - p∆x 2 ] × [y min + p∆y 2 , y max - p∆y 2 ] × [z min + p∆z 2 , z max - p∆z 2 ]. (1.8)
This domain restriction is perfomed so that for any degree p, the chosen number of basis functions is equal to the number of voxels. For a chosen B-spline basis of degree p, p half voxels are removed from the boundaries as illustrated in Figs. 1.2-1.3-1.4. The red lines in these figures define the domain Ω. We justify this choice by the fact that we are usually interested in regions a little bit far from the image boundaries. The image is prolongated outside the boundary of Ω with constant values equal to the boundary values.

• Then, each B-spline function is positionned at the center of each voxel. Each basis function corresponding to a voxel (i, j, k) has the support in x direction equal to:

[x min + i ∆x 2 - p + 1 2 ∆x, x min + i ∆x 2 + p + 1 2 ∆x].
(1.9)

The support in y and z directions are defined similarly.

Binary thresholding based on smooth representations of images

• Finally, the general form of the B-spline representation Ĩ is defined as follows:

∀(x, y, z) ∈ Ω, Ĩ(x, y, z) = kx+p i=kx ky+p j=ky kz+p k=kz c i,j,k β p x -x i ∆x β p y -y i ∆y β p z -z i ∆z (1.10) with kx = x -(xmin + ∆x/2) ∆x - p + 1 2 , ky = x -(ymin + ∆y/2) ∆y - p + 1 2 , kz z -(zmin + ∆z/2) ∆z - p + 1 2 (1.11)
and x i , y i , z i are the voxel coordinate centers defined by x i = x min + (i + 1)∆x/2, y i = y min + (i + 1)∆y/2, z i = z min + (i + 1)∆z/2 for i ≥ 0. We set in practice (x min , y min , z min ) = (0, 0, 0) and ∆x = ∆y = ∆z = 1. When considering images, it is rare that we look for interpolatory conditions at the boundaries of the image in the opposit of physical simulation where boundary conditions need to be imposed. In this case, the Cox-deBoor formalism can allow to obtain these conditions by allowing spans of empty measures (see appendix A for this definition). We also note the the generalization of the precedent spline formalism in order to obtain a regularity smaller than C p-1 is possible using the Cox-deBoor methodology. It is this approach, wich is very used in the computational mechanics community 1.2. Binary thresholding based on smooth representations of images (i.e. the Isogeometric Analysis (IGA) formalism). Using the Cox-deBoor definition (see Eq. (A.4) in appendix A), we can simply define the closed knot vector defined by:

ξ ξ ξ = x min -p ∆x 2 , ..., x max + p ∆x 2
(1.12) subdivided with N + p + 1 uniform knots (or N + p uniform elements with N the number of voxels in one direction) of size equal to the voxel size ∆x. This closed knot vector (i.e. with not repeated first and last knots) allows us to keep a number of basis functions equal to the number of voxels to obtain exactly the same representation as Eq. (1.10) but with the IGA formalism.

We summarise finally Eq. (1.10) by writting it in a global form in the finite element fashion:

Ĩ(x, y, z) = l i=1 m j=1 n k=1 c i,j,k N i,p (x)N j,p (y)N k,p (z) (1.13)
where N i,p , N j,p and N k,p are the total l, m, n univariate B-spline functions of degree p (this set is including zero and non-zero functions at (x, y, z)).

The practical question is: how can we choose the coefficients c i,j,k so that the new representation of the image still represents the data, and is regular enough so that it reproduces a smooth boundary that is similar to the "real" boundary of the cellular architecture ? We present in the following some approaches for setting the grey-level coefficients c i,j,k so that we obtain in the same time a smooth and edge preserving continuous representation of an image. We will first illustrate for each method how a discrete sigmoid function is reconstructed. We consider sampling the function:

S(x) = I min + I max -I min 1 + e -k(x-x 0 )
(1.14) using n points. When imaging a continuous interface located at position x 0 with low resolution, the grey-level variation between the void and solid regions can be modeled with a sigmoid function that is sampled using few voxels (or measurement points). In order to test the different fitting strategies, we set for example, I min = 0, I max = 1 and k = 10 (corresponding to a high slope so that we approximate a step function). In addition, we perturbate this functions using a Gaussian noise with a standard variation equivalent to 10 grey-levels when the signal amplitude is rescaled to [0, 255].

Coefficients setting as the input signal

As a first solution, one can set directly the coefficients c i,j,k = I i,j,k . When setting p = 1, trilinear interpolation is obtained. The data fitting is exact but robustenss to noise is very poor and the function is not smooth enough to represent the edges. In practice, when using trilinear interpolation, images are usually pre-filtered with an adequat noise removal filter. When p > 1, we obtain a smoothing function that does not interpolate the data but acts as a local filter. In the case where images have two characteristic gray-level values (void and solid), this type of approximation might be relevant as it does preserve the dynamic of the input image (thanks to the partition of unity property of B-splines). Fig. 1.5 illustrates this fitting function. We see that the function Ĩ approximates smoothly the noisy image. It is therefore not possbile to preserve the initial image dynamic.

Smoothing splines using the regularized L 2 projection

In order to fit the data or to get a quasi interpolatory function, the L 2 projection can be considered in its regularized form [Unser 1999]:

Ĩ = 1 2 arg min Ĩ∈Sp || Ĩ -I|| L 2 (Ω) + λ 2 R( Ĩ).
(1.15)

The optimality condition of Eq. (1.15) is defined as follows:

(M + λT)c = b (1.16)
where

M i,j = Ω N i,p (x)N j,p (x)dx, b j = Ω N j (x)I(x)dx, (1.17)
T is the Tikhonov operator relative to the regularization energy R and c is the column vector corresponding to the spline tensor coefficients c i,j,k . One of the most common regularizations used in spline fitting is the Thin-plate spline model [START_REF] Craven | [END_REF]:

R( Ĩ) = Ω d i,j=1 ∂ 2 Ĩ(x) ∂x i ∂x j dx.
(1.18)

In this case

T i,j = Ω d k,l=1 ∂ 2 N i,p (x) ∂x k ∂x l ∂ 2 N j,p (x) ∂x k ∂x l dx (1.19)
where d is the dimension of the image. This regularization approach decreases oscillatory effects through the second order derivatives . We note that the data fitting term can be slightly modified in practice as it is sometimes desirable to put more weight on some specific voxels. To do so, a weighting matrix is usually introduced in the L 2 norm. When considering the Thin-Plate spline model, the regularization weight λ reduces oscillatory effects of the signal. If λ is set to 0, the L 2 projection allows to obtain a quasi-interpolatory function. We note that the numerical integration over Ω is performed voxel-wise as each B-spline is centered on the voxel but with the mimial number of Gauss integration points in order to accurately integrate the B-splines. In addition, as the multi-variate B-splines are obtained by tensor product, only integrals in each one-dimensional direction are performed which makes the assembly of the numerical operators more efficient. The integrands are computed as follows:

Ω N j (x)I(x)dx = n i ,n j ,n k i,j,k I i,j,k v i N x i (x)dx v j N y j (y)dy v k N z k (z)dz (1.20)
where v i , v j and v k are the one dimensional elements corresponding to the voxels and N x i , N y j , N z k are the univariate B-spline basis functions. The regularization term of Eq. (1.19) is assembled similarly as in Eq. (1.20). Fig. 1.6 shows the obtained result on the noisy sampled sigmoid function. We observe that when λ → 0, the fitting function nearly interpolates the data and is quite similar to a piecewise interpolatory function. In this case the initial signal dynamic is not preserved. When λ is very large, the solution tends to the least squares line. Setting λ = 0.3 allows to preserve the sigmoid dynamics which proves the interest of adding the Thin-Plate regularization term. The regularized L 2 projection would undoubtedly allow to obtain smooth cell boundaries when applied to two and three-dimensional images. It is however a costly method as it lieds to the resolution of a linear system which size is equal to the image size but the main drawback is that it does not preserves the edges as it nearly interpolates the noisy signal. In addition it is oscillatory even when noise is not considered. Thin-Plate regularization allows to alleviate the oscillatory effect but the slope of the obtained function is very low.

Filtering of the input signal

The authors in [Verhoosel et al. 2015] have suggested a relatively cheaper method by approximating the L 2 projection matrix M by its diagonal through lumping. Instead of solving the linear system, Mc = b (see Eq. 1.16 when λ = 0), one can consider the lumping of M by summing each row of the matrix and solving the diagonal system L(M)c = b where L is the lumping operator: (1.21)

(L(M)) i,i = m j=1 M i,j = m j=1 Ω N i (x)N j (x)dx = Ω N i (x) m j=1 N j (x)dx = Ω N i (x)
The authors in [Verhoosel et al. 2015] compute the coefficients c as it follows :

c i = Ω N i,p (x)I(x)dx Ω N i,p (x)dx . (1.22)
Each coefficient c i is the weighted average of the image over a region of a maximum of (p + 1) d voxels due to the compact support property of B-spline functions. The authors highlighted the fact that this approximates a Gaussian filter as it is known that the cubic B-spline function (illustrated in Fig. 1.1) approximates well a Gaussian kernel function. Suprisingly, the coeffcients c i set by the authors in [Verhoosel et al. 2015] are none other than a special case of a bilateral filter where the chosen kernel is a B-spline function. Bilateral filtering is an image restoration technique that aims at smoothing images while preserving edges. The general idea is to consider a weighted average of the input grey-level data. The general form of this filter is given by the following formula [Tomasi & Manduchi 1998]:

Ī(x) = Ω I(ξ)f (ξ -x)g(I(ξ) -I(x))dξ Ω f (ξ -x)g(I(ξ) -I(x))dξ (1.23)
where f is an intensity kernel used for smoothing the grey-level variations and g is a spatial kernel used for smoothing the distance variations. f and g are usually chosen as Gaussian kernels, i.e.:

f (x) = e -||x|| 2 2 /(2σ 2 d ) , g(x) = e -x 2 /(2σ 2 i ) (1.24)
where σ i and σ d are the bilateral filtering parameters (Gaussians standard variations). The intensity parameter σ i (in grey-level unit) is chosen so that the smoothing is performed for only the pixels that have nearly equal values (and this in order to preserve the edges). The spatial parameter σ d (in pixel or physical unit) specifies the radius of pixels that are averaged. In order to preserve the initial image intensity, a normalization by the weighting functions is performed. Fig. 1.7 illustrates how the bilateral filter acts on a two-dimensional noisy step signal.

Figure 1.7: Schematic explanation of the bilateral filter. Image taken from [Tomasi & Manduchi 1998]. The functions c and s in the caption stand for the kernel functions f and g.

For the sigmoid signal that we consider as an example, the edge is defined in a pixel regions therefore σ d must be greater than 1. The choice of σ i is general set with an estimation of the standard variation of the supposed white noise that acts on the image similarly to any fintering technique. Figs. 1.8-1.9 show the effect of the standard variation σ i on the signal using either the bilateral or a simple Gaussian filter. The more σ i is large the more the filter acts as Gaussian one. When choosing the correct parameter by setting σ i = 0.3, the bilateral filter acts independantly on each part of the step signal which is the desired behavior of an edge-preserving filter. Fig. 1.10 shows that the splines based on the bilateral filtered grey-levels are sharper than the ones with a Gaussian (with the optimal value). Therefore it clearly represents an adequat filter for correcly preserving the edges. Returning to Eq. (1.22), the coefficients c i are therefore the special case where the distance kernel is equal to 1 therefore, this function does not differentiate between the boundaries and homogeneous regions. Rather than using Eq. (1.22) to compute the coefficients c i,j,k , we just compute c = G σ * I where G σ is a Gaussian kernel. In fact this convolution can be efficiently computed using the Fourier transform with very robust algorithms (with available libraries) which is not the case for classic B-spline linear combinations.

We have seen in this section, that the formula (1.10) represents a general and practical framework for representing continuously grey-level images. It allows to obtain different smoothness levels depending on the chosen degree and the method to set the spline coefficients c i,j,k .

Level-set extraction from images

Now that the continuous representation Ĩ is constructed, the contours of the cellular geometry are extracted by simply taking an iso-value as follows:

x ∈ Ω if Ĩ(x) ≥ γ x ∈ Ω if Ĩ(x) < γ , (1.25)
where γ is the grey-level threshold value. The two numerical methods that will be investigated in chapter 2 will rely on this level-set representation. The finite cell method directly uses formula (1.25) as it is a level-set based method whereas the finite element method needs a meshing step. For both visualization and meshing, we consider extracting the contour segments and faces from Ĩ using the marching squares and the marching cubes algorithms. The marching method [START_REF] Lorensen | [END_REF] is the most famous and simple algorithm for extracting the iso-surfaces from a scalar distance field. This method consists in moving a window in the scalar field domain (see ). For each step, the scalar field is evaluated at the nodes of the window. The obtained result is checked in a look up table (see Fig. 1.12) that determines the position of the interfaces which is detemined by linear interpolation. Figure 1.12: The Marching squares and Marching cubes look up tables.

In this work, the general methodology for extracting the contours of cellular architectures is summarized as it follows: first, image denoising is performed depending on the level of noise. The bilateral filter or any other type of filters that do not alter the grey-level values in the neighberhood of the boundaries can be used. Secondly, if the image is sufficiently resolved, one can consider directly the spline coefficients as c = I. If the image is poorly resolved, we consider Gaussian smoothing in order to get a regular and smooth boundary of the cellular architecture. In this case, the spline coefficients are set to the image filtered with a Gaussian kernel (c = G σ * I). The regular image field Ĩ is afterward computed using formula (1.10). Finally the marching cubes is applied to the new field Ĩ in order to extract the boundary segments or faces. We make use of the library scikit-image [van der Walt et al. 2014] to perform this task in Python. The library ITK [Johnson et al. 2013] (Insight Segmentation and Registration Toolkit) can also be used for the same purpose. We illustrate in the following how the spline coefficients of the images influence the regularity of the obtained contours. We consider a virtual image corrupted with a Gaussian noise (of standard variation equal to 10 grey-levels) representing an oscilating curve with a very thin edge (see Fig. 1.13 left). As indicated previously, we compute the image spline coefficients by convolution with a Gaussian kernel. Fig. 1.13 also shows the spline coefficients for different standard variations.

Figure 1.13: From left to right: image 1: initial geometry, image 2: noisy and poorly resolved image of the geometry, images 3-4-5: spline coefficients equal to G σ * I with σ ∈ {0.6, 1, 1.5}.

Fig. 1.14 shows the obtained boundaries for the different standard deviations. The threshold value considered for the marching squares algorithm is equal to (max( Ĩ) + min( Ĩ))/2. We observe that when the Gaussian standard deviation is low, the boundary is more oscillatorty but it is the nearest to the exact boundary. We also observe that the curve smooths too much the very sharp boundaries which can be considered as a limit of this method. The boundary regularization effect of the spline coefficients is also highlighted in Fig. 1.15. We observe that extracting the contours directly on the original image leads to a very noisy and irregular boundary that can itself lead to significant errors related the mechanical solutions that is why we resort to Gaussian filtering as a strategy to reconstruct the smooth cell boundaries. In practice, acquiring images of cellular materials using micro-tomography can lead to poorly 1.2. Binary thresholding based on smooth representations of images resolved images. This can be due to either the experimental setup wich sometimes constrains the distance of the sample from the X-ray source or simply the energy of this source. As an example, Fig. 1.16 shows a slice of the reconstructed volume obtained from a scan of a Rohacell-51 polymetacrylimid foam. In this special case , we observe that the material's boundary is represented with very few voxels.

The cell walls have only a thickness of 2 to 3 voxels. This situation challenges both the gaussian and the bilateral filter when the noise level is very high. We remark in Fig. 1.16 that the noise corrupts two struts that are identifiable with the naked eye (see the two red cirles in the right of Fig. 1.16). This situation reveals that if we fix a low threshold, areas with high noise level are also taken into account in the geometry. In the other case, if a high threshold value is chosen the two cell struts (circled in Fig. 1.16) are removed and a non fully connected geometry is obtained (broken struts). This situation reveals the limits of our segmentation method. In this chapter, we have shown that it is possible to extract smooth and regular boundaries of cellular architectures using binary thresholding when applied to a B-spline representation of images. Provided a sufficient level of resolution, we are able to characterize the cellular domain using a distance field. This spline representation of voxel data has two main applications in this work. The first is related to the construction of image-based mechanical models and the second related to sub-voxel evaluation in DIC and DVC.

Chapter 2

From an image-based level-set geometry to mechanical analysis

In this chapter, we classify image-based mechanical modeling into two categories. The first class is based on the boundary conforming finite element method and the second one is based on fictitious domain approaches. We present a numerical implementation and investigation of these methods and their application as image-based models.

Three aspects play a crucial role in image-based modeling. The geometric aspect which is related to the reconstruction of an observed geometry. The analysis aspect related to the generation of analysis suitable meshes and finally the setting of boundary conditions which is the most difficult aspect especially if one need to fix boundary conditions that fit to an experimental setup in order to perform a dialogue between measurements and simulation.

Finite element image-based mechanical modeling

Concerning the FEM, one of the most common and direct methods is to threshold the image and convert the voxel connectivity into an hexahedral finite element mesh [Hollister et al. 1994, Ulrich et al. 1998, Chevalier et al. 2007]. This is illustrated in Fig. 2.1b where the voxel connectivty is constructed from a magnetic resonance image (MRI) of a trabecular bone. This approach leads to a very irregular geometry which itself can lead to significant errors such as high strain concentrations. It should be noted that one could hope avoiding the segmentation step by using the grey-level value to adjust the voxel behavior [Homminga et al. 2001] but this can lead to huge FE problems and excessive computational burdens especially when whe have very porous materials. When setting over-kill voxel meshes, the development of specific high performance computing tools are to be considered [Liu et al. 2019]. Finally, the most accurate FE approach is the use of tetrahedral finite elements that conform to the complex geometric boundary. To do so, we make use of the marching cubes algorithm to extract the surface and then perform FE meshing [START_REF] Frey | [END_REF], Ulrich et al. 1998] (see Fig. 2.1c). Before proceeding to FE computations, one needs to obtain analysis suitable FE meshes. This step can be time consuming and needs special treatments such as mesh regularization and cleaning, that is why as indicated in chapter 1, the smoothing of the grey-level field is perfomed before applying the marching cubes algorithm. In addition, in order to perform mesh generation, the surface given to the volumetric mesh engine needs to be watertight. In this work, a cuboid region is extracted and closed by assigning the boundaries a grey-level field value smaller than the threshold value γ. Even when denoising images, isolated groups of voxels in the image domain can be segmented as illustrated in Fig. 2.2. These small connected components are removed as they do not contribute to the total stiffness. This procedure can be either performed a priori on the image by binary thresholding the grey-level field and using connected component labeling [Gonzalez et al. 2009] or a posteriori by element search in the generated finite element mesh. Even when performing Gaussian smoothing, the Figure 2.2: Isolated components in a segmented surface mesh of a trabecular bone specimen. The red circles reveal isolated connected components in the mesh connectivity.

presence of poorly shaped and degenerated triangles is still observed. It can lead to instabilities in the mechanical solutions. To illustrate this, we consider a finite element simulation in elasticity on the trabecular bone volume of Fig. 2.1. The bottom surface of the volume is clamped and a compression (using Dirichlet boundary conditions) is performed in z direction. The irregular solution at the degenerated element is revealed in Fig. 2.3. This example shows that a manual intervention should be performed after the simulation in order to detect and "clean" this type of problematic regions. In this work, the different steps of mesh generation are defined as follows: first a smooth grey-level field is constructred using Gaussian smoothing and spline approximation as it was defined in chapter 1. Then the marching cubes is applied to this smooth grey-level field on an over-sampled grid (see Fig. 2.4b).

During the segmentation step, the thresholding value is chosen so that the segmented volume's density corresponds to the experimental density or porosity of the imaged sample (ρ = V /V t ) where V is the cell volume and V t is the total volume (void included)). We obtain afterwards a polyhedric surface mesh which is given to the mesh generator of the CGAL library [The CGAL Project 2021]. It is based mainly on Delaunay tetrahedralization. This library offers the opportunity to fix the element size by varying the mean facet size as shown in Fig. 2.4. Once the mesh is generated, we proceed to a mesh cleaning step that consists in removing the isolated connected components and/or the elements of very small volumes. We should note however the presence of other techniques for enhancing the surface reconstruction step before proceeding to the meshing step [Lizier et al. 2008, Dietrich et al. 2008] but these aspects were not further investigated. One of the main motivations for using fictitious domain approaches is to avoid this meshing step, (which is error prone and time-consuming) and to dissociate the approximation space from the complex geometry, see [START_REF] Parvizian | Jamshid Parvizian, Alexander Düster and Ernst Rank. Finite cell method[END_REF], Embar et al. 2010, Schillinger & Rank 2011, Burman et al. 2015, Lehrenfeld 2016, Fries & Omerović 2016, Legrain & Moes 2018, Kerfriden et al. 2020, Wei et al. 2021,Antolin & Hirschler 2021] to name few. The next section presents these alternative solutions and focus on one of these methods which is the finite cell method, see [START_REF] Parvizian | Jamshid Parvizian, Alexander Düster and Ernst Rank. Finite cell method[END_REF]] for the origin.

Fictitious domain approaches 2.2.1 Brief state of the art

Fictitious domain approaches are a class of methods for solving partial differential equations. They consist in immersing any complex geometry in a more simple domain (usually a cuboid) and represent the solution of the physical problem using the degrees of freedom related to the embedding mesh. When combined with higher order and regular functions, these embedded domain (also referred to as immersed boundary, or unfitted, or even cut) finite element techniques seem to meet their full potential. In fact, they are superior in terms of accuracy than the standard low order boundary fitted strategies (such as the one of previous section) for the same amout of degrees of freedom but at the price of increasing the matrices bandwith. The main advantage of these methods is that they avoid the meshing step. This facilitates the analysis of geometrically complex and timeevolving problems. The immersed boundary method, for example, is originally a computational fluid dynamics method that was initiated in [Peskin 1977] in order to model the blood flow in the heart. We can say roughly that all the fictitious domain methods have in common the same type of geometry which is given as a level-set function and all aim at integrating accurately discontinuous functions (whether related to geometric or material discontinuity). The use of the extended finite element method (XFEM) as a fictitious domain method for treating complex micro-structures was suggested in [Sukumar et al. 2001,Moës et al. 2003,Legrain et al. 2011,Lian et al. 2013] and was applied to porous and two phases materials. It consists in generating first a coarse mesh defined over the whole image domain and refining the mesh in the interface regions. Using the partition of unity, the displacement field is expressed as the sum of two terms; a first term composed of the classical finite element linear combination and a second term composed of the linear combination of enriched degrees of freedom. First, the interface level-set function is interpolated at the mesh nodes, then the degrees of freedom whose support is cut by the material discontinuity are enriched with the use of a special enrichment function that depends of the nodal level-set values [Moës et al. 2003]. The cutFem method [Burman et al. 2015, Hansbo et al. 2017] relies similarly on a background finite element mesh in which the computational domain is immersed. The elements of interest are those which are totally included in the domain and those which cut the interface. The integration of the cut elements is performed by linearizing the level-set function for the approximation of the cutting interface. Afterwards, the cut tetrahedral finite elements are subdivided into integration tetrahedrons on which exact Gauss integration is performed. The particularity of the cutFem method is that it uses a stabilization approach that consists in penalizing the jump of the normal derivatives of the displacement field accross cut edges (or faces). The cutFem has been very recently successfully applied as an imagebased model for solving contact problems between bones [Claus et al. 2021,Mikaeili et al. 2021]. Finally the finite cell method [START_REF] Parvizian | Jamshid Parvizian, Alexander Düster and Ernst Rank. Finite cell method[END_REF]] (FCM), consists usually in embedding the domain in a high order structured grid, integrating accurately the cut elements and penalizing the solution in the fictitious domain by introducing a fictitious stiffness (in linear elasticity). It has the same spirit as the cutFem method. In fact, the tetrahedral version of the finite cell method [START_REF] Varduhn | [END_REF],Duczek et al. 2016] shares many similarities with the cutFem method.

In this thesis, a focus is performed on the FCM because its isogeometric version based on B-splines [Schillinger et al. 2012b,Schillinger & Ruess 2015,Kamensky et al. 2015,de Prenter et al. 2017,Elfverson et al. 2019] opens multiple attractive applications thanks to the structured aspect of the embedding mesh especially for image-correlation as the high order B-splines allow a better description of strains. In the context of image-based modeling, the finite cell method has been successfully applied to the analysis of cellular materials using micro-tomographic images. In [Ruess et al. 2012] and [Schillinger & Ruess 2015], mechanical simulations of a human femure and a metal foam were conducted. Combined to hierarchical splines, the finite cell method was also used to model the micro-mechanical response of a trabecular bone in [Verhoosel et al. 2015]. Very recently, this method has also ben applied to the study of additively manufactured and architected structures [Korshunova et al. 2020, Korshunova et al. 2021].

A major part of our work when reviewing this method has consisted in re-implementing and revisiting the numerical tools related to the FCM so some aspects may be obviously redundant with respect to the existing litterature. However we present a synthesis of the FCM with multiple contributions. After revisiting numerically the FCM in section 2.2.2, in which we verify the convergence properties of the method and also suggest an efficient implementation for the assembly of differential operators, we will try in section 2.2.3 to adjust the modeling parameters (element size, polynomial degree, quadrature rule) so that the discretization error is comparable to that produced through the imaging (i.e., pixelation) of the real specimen. The proposed study in section 2.2.3.2 thus breaks with the usual practice in the field where the numerical convergence is most of time assessed with geometries deemed to be exact (such as in references [Düster et al. 2012, Verhoosel et al. 2015] for instance). This enables us to end up with a "fairly-priced" image-based model in the sense that going to a more refined model will not enable a better accuracy since the total error is dominated by the intrinsic geometry error coming from the image generation process.

The particular case of the FCM

In this work, only elasticity is treated but the method is general and can be applied to any type of mechanical problem.

Stress penalization in linear elasticity

Let Ω p be a domain characterized by a level-set function φ. We denote Γ d the boundary on which a displacement field u 0 is imposed (Dirichlet boundary) and Γ n the boundary on which a loading σ.n = f s is applied. We denote f v the volume forces that act on the domain Ω p . As illustrated in Fig. 2.5, the finite cell method consists in embedding Ω p in a simple rectangular domain denoted Ω and penalizing the stress and forces in the fictitious domain with a spatial penalization function α. The complementary of Ω p is called the fictitious domain and is noted Ω f . Let us consider the principle of virtual work (2.1):

         Ωp σ(u), ε(v) F dx + Ωp f v vdx = Γn f s vds u, v ∈ U ad u = u 0 on Γ d σ.n = f s on Γ n (2.1)
where U ad is the space of kinematically admissible solutions [Ern & Guermond 2002] and A, B F is the Frobenius product defined by A, B F = tr(A T B). tr is the trace of a matrix. The penalization of the stress tensor and volume forces [START_REF] Parvizian | Jamshid Parvizian, Alexander Düster and Ernst Rank. Finite cell method[END_REF]] is given by:

σ α = ασ and f α = αf v (2.2) with ∀x ∈ Ω = Ω p ∪ Ω f , α(x) = α p = 1 if x ∈ Ω p α f = 10 -n << 1 if x ∈ Ω f . (2.3)
The variational problem is subsequently solved using a classic Galerkin approach with high order basis functions defined over a structured mesh as indicated in the right image of Fig. 2.5. With the use of IGA [Hughes et al. 2005, Cottrell et al. 2009] (see appendix A for the details of the method), a set of B-splines can be easily defined over the structured grid and order elevation can be performed. Fig. 2.6 illustrates as an example a structured cubic B-spline grid that embeds the domain Ω p of Fig. 2.5. Similarly to the finite element method, a linear stiffness system is obtained. The elasticity problem is therefore reduced to a linear system Ku = f where u is the vector of displacements at all the control points of the embedding B-spline grid (defined on Ω), f is the force vector and K is the stiffness matrix. The interpolation of the displacement field in the fictitious domain is discarded and the solution has a meaning only in the physical domain Ω p . The challenge and the complexity of the finite cell method is related to the accurate integration of basis functions that have a support that intersects the physical domain, the setting of the boundary conditions when the boundaries are also immersed and finally the treatement of the ill-conditioning due to very small cut elements. If we set α f = 0, the basis functions that have a support lying completely in the fictitious domain must be removed otherwise the system would be singular. If α f = 0, the basis function makes the linear system invertible but still leads to a very ill-conditioned matrix.

Figure 2.6: A structured cubic B-spline embedding mesh.

Integration of cut cells

The integration of the penalized term coming from Eqs. (2.2)-2.3 with the standard Gauss quadrature is no longer valid in a cut element due to the arbitrary shape of the cutting interface. Fig. 2.7 illustrates a penalized integrand on a cut element. In the case of Fig. 2.7a, if the inteface is not taken into account, the integration can be obtained exactly by setting the adequat number of integration points corresponding the regularity of the integrand defined over a rectangle however; in the case of Fig. 2.7b, this problem becomes a general problem of integrating a function on an arbitrary domain. Multiple approaches have been proposed in the litterature in order to tackle this problem. For example, moment fitting approaches [START_REF] Müller | [END_REF], Abedian et al. 2013, Joulaian et al. 2016, Legrain 2021] have emerged in order to integrate functions over arbitrary domains by setting arbitrary integration points and finding the optimal weights of integration. They rely mainly on using polynomial approximations of the integrand. A plenty of other works emerged also in order to create geometrically faithful quadrature rules [Legrain 2013,Nagy & Benson 2015,Kudela et al. 2016,Fries & Omerović 2016,Stavrev et al. 2016, Fries et al. 2017]. In this work, and for computational simplicity, we rather use a spacetree decomposition scheme that allows to approximate the physical domain. Each element of the structured mesh (called a cell in the FCM terminology) is divided into four sub-cells if it cuts the boundary (see Fig. 2.8). The sub-cells that do not cut the geometric boundary are integrated with a full Gauss quadrature (with α p on Ω p and α f on Ω f ). This decomposition is repeated until a predefined maximum level is reached (denote l max ). The last cut cells are finally not taken into account into integration (or wrongly integrated with a Gauss rule). This leads to a very large number of space-tree levels in order to achieve the theoretical convergence (this will be highlighted in the numerical example in section 2.2.2.6). Fig. 2.9 shows the repartition of integration points on each sub-cell. It clearly highlights the fact that Gauss integration has no meaning on the last cut cells. Figure 2.9: Quadtree scheme: integration on a complex domain without accurate integration at the last cut-cells.

Remark 2.2.1. Let us note that, in the context of this thesis, when manipulating B-spline basis functions, a full Gauss-Legendre tensor product quadrature is used on each homogeneous (or uncut) element. The formula of Gauss-Legendre with n points (n ≥ 1) is exact for polynomial functions of degree 2n -1. A minimal number of integration points to integrate a p polynomial function is set to

n gp = p + 1 2 (2.4)
whrere ⌈x⌉ is the nearest integer larger than x (ceiling function). We consider n gp integration points in each direction for a multivariate polynomial function. For a problem of dimension d, we integrate an element-wise polynomial function with n d gp tensor product Gauss points. The multivariate Gauss weights are defined as the tensor product of the univariate weights.

In order to improve the accuracy of the integration at the last cut cell and to decrease the number of sub-cells, we choose to adopt a triangulation (or tetrahedralization) of the last cut cells using a strategy similar to the ones suggested in [Marco et al. 2015, Verhoosel et al. 2015, Kudela et al. 2016]. The oldest paper found and that uses the integration over cut-cells with triangulation is [Oevermann et al. 2009] and concerns the finite volume method. We proceed as follows for the tetrahedralization of cut cells: first, the level-set function φ is evaluated at the vertices of the cut cell (see Fig. 2.11a). Secondly, the intersection between the interface and the cell edges is estimated by setting the isozero of the level-set by linear interpolation exactly as it is done in the marching cubes algorithm (see Fig. 2.11b). Finally, the set composed of the interface vertices in addition to the vertices of negative level-set values are given to a tetrahedralization algorithm. We make use of the Delaunay algorithm provided in the CGAL library [The CGAL Project 2021]. It takes as an input a set of vertices and performs the triangulation or the tetrahedralization. The advantage of this method is its ease of implementation but the main drawback is that a set of tetrahedrons is obtained with no information on the boundary of the cutting interface. With this method, the setting of Neumann boundary conditions is not possible as its computation depends on the computation of interface surface integrals. To perform this need, we rather suggest performing the tetrahedralization manually by simply calling the Marching cubes sub-routines as they list all the possible cutting configurations (see again Fig. 1.12). We have to raise the fact that all this triangulation strategies add a certain computational complexity from an implementation point of vue and seem to discredit a little bit the argument that consists in describing the finite cell method as a "meshing-free" method. We must emphasize however that this is only related to the geometric aspect which is not the case of the finite element method.

Boundary conditions

When the embedding domain boundary does not conform to the physical boundary, Dirichlet boundary conditions can no longer be applied. In this work, the simplest case is treated, i.e. the boundary of the embedding mesh is aligned with the Dirichlet exterior boundaries. This corresponds to the situations illustrated in Figs. 2.12a-2.12b. In the general case, when the boundary conditions are imposed on boundaries that are included in the embedding mesh (for example, for the situation in Fig. 2.12c) , a weak enforcement of boundary conditions is used. When using the Nitsche method, the terms of the principale of virtual work are extended by a weighted (or penalized) residual term for displacement and unknown reacting forces (see [Schillinger & Ruess 2015] page 12 for the details of the equations). An additional stabilization term is also added in order to guarantee the positive definiteness of the stiffness matrix [Ruess et al. 2013, Schillinger & Ruess 2015, Antolin et al. 2019].

Finally Neumann boundary conditions (such as those displayed in red in Fig. 2.12b) are imposed using the integral term in the principale of virtual work of Eq. (2.1). This term is computed for the degrees of freedom whose support intersect the Neumann boundary.

Efficient assembly method for the B-spline finite cell method

When no parametric transformation is considered, the integration of B-spline basis functions can be performed in each one-dimensional direction thanks to the tensor-product definition . Let N be one multi-variate basis function (corresponding to one B-spline control point). The partial derivative of (c) Dirichlet boundary conditions on the immersed boundary (This case is not treated in the computational code).

Figure 2.12: Some configurations of boundary conditions when using a structured embedding mesh.

N can be efficiently computed by the following formula:

Ωe ∂ d i=1 β i N (x) ∂x β 1 1 ...∂x β d d dx = Ωe d l=1 d β l N l dx β l l (x i )dx = d j=1 Ωe j d β i N i dx β i i (x j )dx j , (2.5)
where Ω e j is the univariate element corresponding to the element Ω e and (N l ) l are the univariate basis functions corresponding to N . This general principale of efficient assembly for IGA is also known as sum factorization [START_REF] Antolin | [END_REF]. By exploiting the tensor product of B-splines, we only need to compute the integral of univariate basis functions. This makes the assembly of discrete differential operators very efficient. We can transpose this idea in order to assemble the stiffness matrix or any discretized symmetric form. When considering the symmetric form resulting from the weak formulation of a boundary value problem, a discretization using the Galerkin method leads to symmetric operators equal to the tensor product of shape function derivatives operators (with respect to the physical coordinates). Using the same idea of Eq. (2.5), one can efficiently assemble the block of the integrated differential operators. Each univariate element of the discretization is supported by a finite set of basis functions that can be assembled in one vector function denoted φ φ φ l (x) and defined by φ φ φ l (x) i = φ l i (x) where (φ i l ) i is the set of univariate basis functions of the dimension l. In order to define the multi-variate basis functions vector denoted φ φ φ, we introduce the flatten kronecker product denoted by the symbol ⊗ v and defined as follows:

(a⊗ v b) j+n(i-1) = a i b j ∀i ∈ 1, m ∀j ∈ 1, n , a ∈ R m , b ∈ R n . (2.6)
This operation is a flattening of the tensor product of two vectors defined by

(a ⊗ b) ij = a i b j (2.7)
and its general formula defined for matrices:

A ⊗ B =    A 1,1 B . . . A 1,n B . . . . . . A m,1 B . . . A m,n B    ∈ R mp×nq , for A ∈ R m×n and B ∈ R p×q .
(2.8)

For example, if tri-variate B-splines are considered, the vector of basis functions φ φ φ is defined by φ φ φ = φ φ φ 3 ⊗ v φ φ φ 2 ⊗ v φ φ φ 1 with φ φ φ 1 , φ φ φ 2 , φ φ φ 3 the vectors of univariate B-spline basis functions for the directions x, y and z.

For any symmetric form, a contribution to a block of an elementary differential operator can be expressed as follows:

Ωe ∂ d i=1 α i φ φ φ(x) ∂x α 1 1 ...∂x α d d ⊗ ∂ d i=1 β i φ φ φ(x) ∂x β 1 1 ...∂x β d d dx.
(2.9)

In order to simplify the computation of (2.9), we can remark that the flatten tensor product verifies the following property:

(a⊗ v b) ⊗ (c⊗ v d) = (a⊗c) ⊗ (b⊗d), for a, c ∈ R n and b, d ∈ R m .
(2.10) Using Eq. (2.10), the block integrated differential operator (2.9) can be simplified as follows:

Ωe ∂ d i=1 α i φ φ φ(x) ∂x α 1 1 ...∂x α d d ⊗ ∂ d i=1 β i φ φ φ(x) ∂x β 1 1 ...∂x β d d dx = Ωe l=1 v l=d d α l φ φ φ l (x l ) dx α l l ⊗ l=1 v l=d d β l φ φ φ l (x l ) dx β l l dx = Ωe 1 l=d d α l φ φ φ l (x l ) dx α l l ⊗ d β l φ φ φ l (x l ) dx β l l dx = 1 l=d Ωe l d α l φ φ φ l (x l ) dx α l l ⊗ d β l φ φ φ l (x l ) dx β l l dx l .
(2.11)

The classical elemental assembly consists in computing the term (2.9) by summing the n d gp block matrices (which are themselves computed by a tensor product) whereas with the suggested simplification, only d tensor products of the integrals of one-dimensional differential operators is needed. We can apply this method for the assembly of homogeneous finite cells (i.e those which are not cut) in the octree scheme which makes the assembly much faster. As an example, suppose we want to assemble a Laplacian over the domain shown in Fig. 2 (2.12)

The first step consists in computing once for all the one-dimensional integrals for all the quad-tree levels. Suppose that the quadtree decomposition algorithm is located at an homogeneous element (see the green element e shown in Fig. 2 We note that this approach can be applied to any tensor product basis. It allows to improve the octree scheme because all the homogeneous cells are precomputed for all octree levels as illustrated in Fig. 2.14. The assembly of the differential operators, including the octree decomposition is performed using a shared memory parallelism (with OpenMp [Chapman et al. 2007]) in a C++ environement and the interfacing with the Python environment is performed with the SWIG tool [Beazley et al. 1996] 1 .

Conditioning concerns and numerical solving

When setting α f = 0, the basis functions that have a support entirely located in the fictitious domain are removed. But this is still not sufficient when solving the problem's linear system using iterative solvers. The use of the conjugate gradient method fails in this case as the obtained matrices are (2p+1-2/d) where η is the smallest non-zero volume fraction of the cutelements, p is the order of the basis functions and d the dimension of the problem. This shows that the conditioning of the matrix becomes too large if very small cut-cells exist and/or the order of basis functions is high. The authors have therefore constructred an efficient schwarz preconditioner to tackle this problem [de Prenter et al. 2017]. Similar partial observations were also obtained in [Elfverson et al. 2018], where the authors suggested to remove basis functions based on a volume indicator. The same idea was defined earlier in [Verhoosel et al. 2015] where the authors remove basis functions based on the volume fraction of their support. In the case of our implementation of the finite cell method, we proceed to the fictitious stiffness strategy that consists in setting α f = 0 and only remove the basis functions that have a support that is completely in the fictitious domain (i.e. Ω p ∩ Supp(N i,p ) = ∅).

In this case we solve the obtained linear systems using either LU decomposition (using the SuperLU solver provided by the Python scipy library). If the boundary conditions are correctly set, i.e. the linear system to solve is still symmetric positive definite (see Eq. (A.73) in appendix A), we rather use a Cholesky decomposition (using the CHOLMOD solver provided by the Python scikit-sparse library). It was also indicated in [START_REF] Heinze | [END_REF], that a block-diagonal Jacobi preconditioned conjugate gradient method was applied in the case of using penalized stiffness (α f = 0) but this has not been tested in this work. For the basis function removal aspect, as the geometry can be complex, we sugggest here some general indicators that allow to quantify the influence of each basis function. We denote I the index set of the basis functions to remove.

• Volume fraction indicator: This indicator, suggested in [Verhoosel et al. 2015], consists in computing the volume fraction of the support of a basis function. It is defined as follows

I = i ∈ 1, m / s(i) = V(Supp(N p i ) ∩ Ω p ) V(Supp(N p i )) ≤ ε , (2.17)
where V(x) is the volume of the domain x and m is the number of basis functions of the B-spline mesh (we recall that it is equal to the number of control points).

• Basis function variations: The indicator of Eq. (2.17) is efficient only if the elements cut the domain without being "crossed" by the volume as shown in Fig. 2.15. As we rather choose relatively coarse B-spline elements, we generalize the indicator of Eq. (2.17) by taking into account the variation of each basis function over the cut domain:

I =          i ∈ 1, m / s(i) = Supp(N p i )∩Ωp N p i (x, y)dxdy Supp(N p i ) N p i (x, y)dxdy ≤ ε,          .
(2.18)

Figure 2.15: Example of a basis function that is removed with the indicator of Eq. (2.17) but not removed with the indicator of Eq. (2.18). The red square in the support of the basis function and the red node is the corresponding control point.

We note that this criteria is equal to the criteria (2.17) when N p i = 1.

• Stiffness summation: One calso consider a column summation of the stiffness matrix

I =    i ∈ 1, m / s(i) = m j=1 |K ij | ≤ ε,    . (2.19)
As we set ε = 0, all this criterias are equivalent but we make use of the criteria of Eq. (2.18). Once the index set I determined, the final system is reduced by removing the degrees of freedom corresponding to the nodes in I.

Numerical examples

In order to verify the different convergence properties of the finite cell mehod, we consider a first numerical example for an exact problem. For this purpose, we start by considering the geometry presented in Fig. 2.16. This test case is equivalent to solving an infinite open-hole plate in tension.

We fix the material properties as follows: E = 10 5 P a and ν = 0.3, where E is the Young's modulus and ν is the Poisson ratio. When considering a far-field traction denoted T x that we fix to 10P a, the

-1 0 1 2 3 4 5 6 -1 0 1 2 3 4 5 6 σ. n = σ ex . n σ. n = σ ex . n u y = 0 u x = 0 R = 1 L = 4
(a) Mechanical problem definition. stress field can be expressed in an analytical form :

     σ ex xx (r, θ) = cos(θ) 2 σ ex rr (r, θ) -2 cos(θ) sin(θ)σ ex rθ (r, θ) + sin(θ) 2 σ ex θθ (r, θ) σ ex xy (r, θ) = sin(θ) cos(θ)σ ex rr (r, θ) + cos(θ) 2 -sin(θ) 2 σ ex rθ (r, θ) -sin(θ) cos(θ)σ ex θθ (r, θ) σ ex yy (r, θ) = sin(θ) 2 σ ex rr + 2 cos(θ) sin(θ)σ ex rθ (r, θ) + cos(θ) 2 σ ex θθ (r, θ) (2.20) where                  σ ex rr (r, θ) = T x 2 1 - R 2 r 2 + T x 2 1 + 3R 4 r 4 - 4R 2 r 2 cos(2θ) σ ex θθ (r, θ) = T x 2 1 + R 2 r 2 - T x 2 1 + 3R 4 r 4 cos(2θ) σ ex rθ (r, θ) = - T x 2 1 - 3R 4 r 4 + 2R 2 r 2 sin(2θ)
.

(2.21) Formula (2.21) was extracted from [Cottrell et al. 2009] page 116. We solve this problem using the finite cell method by setting a structured B-spline grid which boundaries conform to the boundaries of the plates (see Fig. 2.16b). This makes the setting of the boundary conditions easier and no weak imposition of Dirichlet boundary conditions is considered. Based on the exact solution defined over the exact domain Ω ex , we consider the error in energy norm defined as:

E h = Ω ex σ ex -σ h , ε ex -ε h F dx Ω ex σ ex , ε ex F dx , (2.22)
where ex and h subscripts stand for, respectively, the exact and the approximated solutions. The numerical parameters related to the finite cell method are the element size, the degree of approximation, the penalization parameter α f and the number of quad-tree levels l max . Let us perform some numerical experimentations based mainly on element refinement.

• Influence of the stress penalization factor: We consider varying the element size for quadratic and cubic B-splines and look for the evolution of the approximation error for different penalization parameters. Fig. 2.17 shows the obtained results. We observe that for the quadratic approximation, the convergence is the same for values of α f smaller or equal to 10 -5 (see Fig.

2.17a

). For the cubic approximation, the convergence is the same for values of α f smaller or equal to 10 -7 . The penalization must be set very small when choosing a high order B-spline basis.

As the support of B-spline basis functions increases with the order, fictitious basis functions influence more the physical domain. Figure 2.17: Influence of the penalization parameter under h refinement.

• Influence of the quadtree integration depth: We investigate the effect of the integration of the last cut-cells by varying the number of quad-tree levels with and without the triangulation of the last sub-cells. Fig. 2.18 shows that, in order to achieve full convergence, a large number of subcells is needed when no triangulation is performed (l max = 9) which is very costly and only l max = 4 when the triangulation is performed. Fig. 2.18 also shows that if the element is sufficiently large, the error related to the mesh is dominant over the integration error and the number of sub-cells can be reduced.

• Element refinement and order elevation: Finally, we consider the convergence curves when setting α f = 10 -10 . Under h refinement, the theoretical convergence order with a degree p is equal to

O(h p ) = O(m) -p/2
, where m is the number of degrees of freedom and h is the size of the element [Dauge et al. 2015]. This is verified by the curves obatined in Fig. 2.19 where each curve has an approximate slope of -p/2. Also, for the same number of degrees of freedom, high order B-splines allow to attain very low errors compared to triangular finite elements. Fig. 2.19 shows that a very poor element discretization using quadratic B-splines allows to attain the same error obtained using a very fine finite element mesh with less than one order of degrees of freedoms. Fig. 2.20 illustrate this large difference in terms of number of elements. In this

10 3 10 4 Number of degrees of freedom 10 -5 10 -4 10 -3 10 -2 H 1 semi-norm error lmax = 0 lmax = 1 lmax = 2 lmax = 3 lmax = 4 lmax = 5 lmax = 6 lmax = 7 lmax = 9
(a) No accurate integration at the last cut-cells. work, we only consider uniform refinement, but the method can be even more ameliorated by considering local refinement strategies. Hierarchical B-splines, which is a method for nesting different B-spline spaces [Giannelli et al. 2012,D'Angella et al. 2018,D'Angella et al. 2018,Garau & Vázquez 2018] had been successfully applied in the context of the finite cell method [Schillinger & Ruess 2015, Verhoosel et al. 2015]. Figure 2.20: Two discretizations that give the same error value.

Remark 2.2.2 (Three-dimensionional extension). When using high order B-splines in high dimensions, the number of non-zero values of the sparse operators increases which may represent a practical limit. For a problem defined on a space of dimension d, the number of non-zero values is proportional to (p + 1) 2d whereas it is proportional to 4 × 4 = 16 for linear tetrahedral elements, that is why resorting to very high order B-splines seems to be costly for very large problems.

We finally consider a three-dimensional example that consists in analyzing the elastic behavior of a curved beam. The level-set function that describes a torus is given by the equation

φ(x, y, z) = (x -x 0 ) 2 + (z -z 0 ) 2 -R 2 + (y -y 0 ) 2 -r 2
, where R is the distance from the torus center (x 0 , y 0 , z 0 ) = (0, 0, 0) to the center of the tube and r is the radius of the tube. We set R = 4.5, r = 0.5. The bottom face of the geometry (plane z = 0) is clamped and the opposite face (plane x = 0) is submitted to a shear displacement in z direction (u z = u 0 = -1). Fig. 2.21 illustrates the geometry and the boundary conditions. We set the Young's modulus value E to 2 × 10 4 and the Poisson ratio ν to 0.3. We solve this problem using 6 elements for each direction and with 2 levels for the octree integration scheme. Fig. 2.22 shows the finite cell mesh in addition to the octree integration domain. We compare qualitatively the solution using the finite cell method with the solution obtained from a NURBS torus composed of 1400 quadratic elements. We refer the reader to section A.2.6 in appendix A for the details concerning the generation of the exact geometry. Fig. 2.23 shows the Von Mises stress field (defined by Eq. (2.23)) obtained for both solutions. We observe that one cubic element encompassing the beam thickness seams to be sufficient enough to accurately represent the bending of the beam.

σ vm = 3 2 σ - 1 3 tr(σ)I F = 1 2 [(σ xx -σ yy ) 2 + (σ yy -σ zz ) 2 + (σ zz -σ xx ) 2 ] + 3(σ 2 xy + σ 2 yz + σ 2 xz ) (2.23)

Application as an image-based model

Even though meshing a distance field can be quite automatic with existing mesh generators, manual intervention is still a problem when geometric irregularities exist in the finite element mesh. The finite cell method seems to be more robust for arbitrary complex geometries because the geometric artefacts (i.e. the critical finite elements) are somehow separated from the basis of approximation.

In practice, the finite cell method takes as an input a level-set function coming from any type of geometric representation (Stl file for example or an image or even a finite element-mesh). We can say roughly that all the complexity of any fictitious domain code lies in the leve-set evaluation which is located in the integration rootine. In practice, geometry requests to arbitrary complex geometries can be characterized by efficient data-structures such as Axis Aligned Bounding Boxes (AABB trees).

The CGAL library [The CGAL Project 2021] allows to perform this need efficiently. The idea is to give to the fictitious domain method a geometry tree object and call the distance computation rootine when needed. It is also possible to incorporate a pure voxel model as the one depicted in Fig. 2.1b in a finite cell mesh. To do so, one need to position the B-spline grid so that the elements of the mesh must contain an integer number of voxels, i.e. an element boundary cannot intersect a voxel. In this case, each cell is reduced to a binary matrix of zeros and ones. The integration of the cell is therefore reduced to the decomposition of a binary matrix into the union of rectangular homogeneous blocks which is a problem of pattern recognition [START_REF] Kulikov | Decomposition of binary images-A survey and comparison[END_REF], Subercaze et al. 2016]. The binary segmentation allows to approach the true geometry only if the image is highly defined (see again Fig.

2.1b that illustrates the irregular voxel geometry resulting from a poorly resolved image). In this work, we rather use the smooth level-set as described in chapter 1. This is the starting point of the finite cell method applied as an image-based model.

Parameter adjustements for fairly-priced image-based mechanical modeling

The content of this subsection is reproduced from [START_REF] Rouwane | [END_REF].

After presenting the finite cell method, its basics and its numerical implementation, we focus now on the application on arbitrary cellular geometries obtained from grey-level images. We will perform first a numerical investigation of the convergence properties when considering inexact geometric models, in particular reconstructed geometries from grey-level images. We rely on the same geometry as the one shown in Fig. 2.16. Let us assume that we capture a representative image I of the sample, i.e., in our case, that we record an undistorted fronto parallel image or reconstruct a tomographic image from X-ray projections. Such an image usually results from a relatively complex process, particularly in the case of tomography. However, the latter is generally composed of a set of regularly arranged pixels or voxels carrying information that may at first be related to a physical property. The contrasts in the image plane are here assumed to result only from the underlying architecture. We consequently assume that the observed texture is binary and that it is simply quantified by the voxel. The idea is to build, from the set of grey-level data set grabbed in a region of interest, a mechanical model that accounts at best for the assumed local behavior of the real sample. Fig. 2.24 summarizes the different steps for the construction of a mechanical image-based model. In the context of this section, where only synthetic tests are performed, the images are generated by a simple algorithm that consists in assigning, for a given pixel region, a grey-value that represents the volume fraction and then quantizes it. More precisely, this grey-value is determined by computing the distance of the pixels (or voxels) to the reference geometry. We must emphasize that this rasterization procedure does not represent the sophisticated mechanism of real image acquisition and that other more realistic models exist, those including filling factors for instance [START_REF] Orteu | [END_REF]]. Other external parameters such as noise can also induce non negligible errors. In this first study regarding the image-based model construction, the effect of noise is not examined for the sake of clarity even though it was partially discussed in chapter 1. It will, however, be taken into account later in the second part treating DIC. In the following, we first present the process of constructing a geometric model from grey-level data (see step 2 in Fig. 2.24). Based on this geometric model, we then construct an automated analysis model that approximates the elastic behavior of the image being observed (see step 3 in Fig. 2.24). Finally, a numerical convergence study is conducted in order to properly fix the parameters of the constructed mechanical model (element length, order of approximation and quadrature rule) such that the error obtained at step 3 of Fig. 2.24 is in the same order of magnitude as the one coming from the intrinsic geometry discrepancy produced through steps 1 and 2 of Fig. 2.24. In that sense, we will characterize this image-based model (i.e., including the correct parameter settings) as a fairly priced model.

Analysis of the geometry error and its propagation in the mechanical solution

We now apply the B-spline variant of the finite cell method to synthetic images using a given levelset function obtained from the spline representation of images. We start by conducting a geometry convergence analysis in order to fix the sufficient number of quad-tree levels that should be chosen practically and then proceed to a numerical investigation of the finite cell method applied to the level-set coming from an image. This investigation will consist in setting the other discretization parameters (mesh size h and polynomial degree p). We recall that our objective is to end up with a fairly-priced mechanical model in the sense that going to a more refined model will not offer better accuracy since the total error will become dominated by the intrinsic geometry error coming from the image generation process.

Numerical investigation of the level-set geometry error

To start with, we consider a very simple, one-dimensional example. More precisely, we undertake to compute the length of the support of a step function represented by a one-dimensional image. Given an analytical step function V defined by:

V (x) = v max x ∈ [0, c[ v min x ∈]c, l]
with v max > v min and c < l, (2.24)

we assume that we have acquired a one-dimensional image defined by a set of pixels (g i ) i∈ 0,n and that represents a sampled version of V (see Fig. Starting from a one-dimensional mesh, we apply the one-dimensional version of the quad-tree scheme of Fig. 2.10 in order to estimate the value of the true interface c. We denote this quad-tree approximation of the interface c a . In addition, we denote c the abscissa of the threshold value of the smoothing function (see Fig. 2.25 right). We recall that the smooth representation of the image is defined by formula (1.10). As (1.10) cannot be inverted directly, c is found using a Newton-Raphson algorithm in a small neighborhood near the interface. The geometric approximation of the interface induces that c = c = c a (see Fig. 2.26). We introduce the three errors listed below:

• Intrinsic geometry error:

E = c - c c , (2.25)
it represents the difference between the true interface and the exact value of the "false" (or modeling) interface.

• Total geometry error:

Ē = c a -c c , (2.26)
it represents the difference between the estimation of the false interface (using the onedimensional quad-tree scheme) and the true interface.

• Domain integration error:

Ẽ = c a - c c , (2.27)
it represents the quad-tree integration error with respect to the false geometry. More precisely, this error comes from the linearization of the interface in the last sub-cell. In other words, it measures how we accurately integrate the false geometry. The position of the interface c within one pixel influences the accuracy of the geometry approximation. Fig. 2.27 shows that the acquisition process induces an error when the pixel does not intersect the interface symmetrically. The error is minimal when the interface is located exactly at the center or at the boundaries of a pixel. The maximum error is about 5 × 10 -2 pixels (i.e. a 2% relative error). Another important user-defined parameter for estimating the interface is the threshold value. A pragmatic choice is adopted here: we set this value (see Eq. (1.25)) to the mean grey value:

γ = v min + v max 2 .
(2.28)

The threshold γ is dependent of the chosen level-set method (i.e. the manner how the physical representation of the image is constructed) and the image resolution. The determination of the threshold criteria is still a practical open question. Experimental calibration strategies could help adjusting the threshold criteria with respect to physical properties of the imaged specimen. The interface c a is calculated by summing the lengths of the integration sub-cells. We recall that we perform a linearization of the level-set at the finest level (one-dimensional version of the triangulation shown in Fig. 2.10). Denoting x min and x max the extremities of the last sub-cell, the linear approximation of the interface denoted c a can actually be computed as:

c a = γ -b a , with a = Ĩ(x max ) -Ĩ(x min ) x max -x min and b = Ĩ(x min )x max -Ĩ(x max )x min x max -x min (2.29)
where Ĩ is the continuous representation of V displayed in Fig. 2.25-right.

In order to fix the sufficient sub-cell level, we conduct a convergence study using the different errors E, Ē and Ẽ defined in (2.25), (2.26) and (2.27), respectively. We vary the quad-tree algorithm's maximum level. Due to the influence of the location of the boundary interface as shown in Fig. 2.27, each error considered is computed for the worst acquisition configuration, i.e., the configuration in which the geometry error is maximum (see Fig. 2.27 again). Fig. 2.28 shows the evolution of the two errors Ē and Ẽ for two image resolutions. The corresponding intrinsic geometry error E is also put on this figure. We observe that the error Ē converges to the intrinsic geometry error E and that the image resolution improves the geometry approximation (the intrinsic geometry error is reduced from left to right). The integration error Ẽ goes down to zero and intersects the geometry error E for a sub-cell size equal to the pixel size. This induces that there is no interest for refining the domain of integration beyond the size of a pixel since the intrinsic geometry error is already achieved. In other words, going to higher quad-tree integration levels would necessary lead to over-calculations. with respect to the size of the smallest sub-cell in pixel size units (one-dimensional test case).

We then reproduce the previous numerical experiment for a two-dimensional case by considering the geometry depicted in Fig. (2.16). The image is generated with a grey-level filling algorithm based on volume fraction similar to the one used in the one-dimensional example related to Eq. (2.24). The pixels that are respectively inside and outside the circle's boundary are assigned with a minimum grey-level value v min and a maximum grey-level value v max . The pixels that intersect the circle's boundary are assigned with a mean grey-level value. Using the level-set method, we approximate the circle boundary by setting again a threshold value defined by (2.28). In order to quantify the geometry error, we detect a large number of radii spaning θ min = 0 to θ max = π/2 of the smoothing function by applying a Newton-Raphson algorithm on the radial axis. The exact geometry for Fig. 2.16 is defined for the exact radius R = 1 and a length of the plate equal to L = 4.

We investigate the convergence of the quad-tree algorithm by considering the approximation of the area of the plate. Similarly to the one-dimensional study, we define the following errors, i.e. the intrinsic geometry error E, the total geometry error Ē and the domain integration error Ẽ such that:

E = | Ã -A | A , Ē = |A a -A| A Ẽ = |A a -Ã| Ã .
(2.30)

A a is the approximation of the area bounded by the level-set. It is equal to the sum of areas of the integration sub-cells and triangles. Ã is the area of the level-set geometry defined by:

à = L 2 - θ=θmax θ=θ min r=R(θ) r=0 rdrdθ, (2.31)
where R(θ) is the radius of the level-set approximated circle, which depends on the angle θ. Ã can be numerically computed by:

à = L 2 - θ max -θ min 2n θ n θ i=1 R 2 i , (2.32)
where each R i is the radius abscissa, i.e. such that the smoothing function is equal to the threshold iso-value, and n θ is the number of angles considered (here n θ = 10 6 ). A is the exact area of the plate defined by:

A = L 2 - (θ max -θ min )R 2 2 .
(2.33)

In Fig. 2.29, we plot the evolution of the errors E, Ē and Ẽ with respect to the size of the smallest integration sub-cell, similarly as in Fig. 2.28. We perform the numerical test for two different image resolutions. We obtain the same behavior as in the one-dimensional case. As a consequence, a sufficient quad-tree level can be chosen so that the smallest integration element has approximately a size equal to the pixel size. We generalize this observation, without any proof, to other image dimensions. For a dimension d, the maximum space-tree levels for an element e is given by the following rule:

V(e) 2 dlmax = ∆x ⇒ l max (e) = 1 d log 2 V(e) ∆x (2.34)
where V(e) and ∆x are respectively the volume of the element and the voxel. In the rest of this work, the number of integration sub-cells will be fixed by the constant defined in (2.34). In order to study local curvatures influence, the same analysis was undertaken on more complex geometries. For instance, polar curves defined by R(θ) = R+0.1 sin(ωθ) were considered (see Fig. 1.13 for the geometry). For the tests performed with ω = 8 and ω = 16, a similar error evolution as the one shown in Fig. 2.29 was obtained. However, as it was highlighted in chapter 1, the level-set function based on Gaussian smoothing fails to represent sharp or very highly curved boundaries.

Numerical convergence of the image-based model: adjustement of the degree and element size

After identifying the sufficient number of quad-tree levels, we study in this section the numerical convergence of the mechanical solution of the image-based level-set geometry. For this purpose, we consider the same geometry and boundary conditions as the one considered in Fig. 2.16. In the following, we investigate the convergence of the image-based immersed boundary model defined previously under mesh refinement and confront it to three other reference FE image-based models listed below:

• Q 4 voxel-based model: Threshold of the initial image and conversion of the pixel connectivity into Q 4 finite elements after a binary segmentation [Ulrich et al. 1998] (equivalent to Fig. 2.1b).

• Finite element mesh of the Marching squares boundary: Contour detection using the linear marching squares algorithm and meshing of the closed curve boundary (equivalent to Fig.

2.1c).

• Mechanical properties based on grey-levels: Converts each image pixel to a Q 4 finite element which has an elastic behavior based on the pixel's grey-level value. More precisely, the model's Young's modulus (denoted here Y ) can be defined as a function of grey-levels (here set linear):

Y (v) = v -v min v max -v min E max + v max -v v max -v min E min (2.35)
where v is the element's grey-level value, E max = E and E min = αE with α a penalization parameter equal to α f = 10 -8 . This method has always a number of elements equal to the number of pixels and therefore induces a large number of degrees of freedoms. However, it has a particular interest when dealing with heterogeneous materials that have different material properties [Liu et al. 2019].

Based on the exact solution defined over the exact domain, we consider the error in energy norm defined by (2.22). This test is exactly the same as the one conducted in Fig. 2.19 but this time with an erroneous level-set constructed from the image. Fig. 2.30 shows the results obtained in terms of convergence under mesh refinement for two image resolutions. More precisely, "deg 1 Fcm", "deg 2 Fcm", "deg 3 Fcm" and "deg 4 Fcm" denote the solution of the proposed image-based immersed boundary model when considering the quad-tree level of Eq. (2.34) and for polynomial degree p = 1, p = 2, p = 3 and p = 4, respectively. We observe that the theoretical convergence rate is obtained. The convergence is again of order O(h p ) ≈ O(m) -p/2 , where m is the number of degrees of freedom and h is the equivalent size of the element. However, the high order solutions converge to an asymptotic value. This result breaks with what is more usually encountered in the field where the convergence of numerical methods is most of time assessed from resulting geometries deemed to be exact (e.g., see [Düster et al. 2012, Verhoosel et al. 2015]). Here, since we take into account the additional error due to the pixelation of the image (see step 1 in Fig. 2.24), we are able to observe an asymptotic constant behavior which means that there is no need to go to too fine models for optimal accuracy. Indeed, after a certain refinement level, the total error becomes dominated by the intrinsic geometry error coming from the initial image generation process. As it was mentionned for exact geometries, Fig. 2.30 shows that the image-based finite cell model is significantly more accurate (for a given number of degrees of freedom) than the three voxel-based methods represented by the triangular dots and the cross dot. In fact, the cubic variant of the proposed model allows to gain approximately one order on the precision for the same amount of degrees of freedom. This accounts for the use of the finite cell method as an image-based model in order to obtain accurate solutions. More precisely, we will consider p ∈ {2, 3} for our model since it seems that we do not increase the accuracy with higher orders. This appears in contrast with what is usually advised in the current literature, which is, once again, due to the fact that we take into account the pixelation error in the process here. Regarding the size of the element, the idea would be to take the one corresponding to the beginning of the stagnation of the error. It depends on the test case. However, we will see in section 2.3 that for open-cell materials with one sub-cellular scale, we can consider a size of element such that it encompasses the thickness of a fiber to correctly reproduce local bending.

Towards the analysis of arbitrary complex geometries

In order to come closer to real cellular micro-architectures, we consider two other numerical examples in two and three dimensions, for which we compare the different solutions visually.

A cellular hexagonal structure

The content of this subsection is reproduced from [START_REF] Rouwane | [END_REF]. Fig. 2.31 shows an image of a repeated hexagonal pattern. One of the applications of image-based modeling would be for example the extraction of windows of different sizes from complex cellular architectures in order to determine an elementary representative volume. We consider in this example extracting from the image of Fig. 2.31 a hexagonal window. We investigate the mechanical problem defined in Fig. 2.32a where the geometry is subjected to the displacement load u 0 = -1. The center

ε vm = ||ε|| F = tr(ε T ε) = ε 2 xx + ε 2 yy + 2ε 2 xy .
(2.36)

We compare the strain norm for the different methods in Fig. 2.30. In addition, we put on the right of each figure the reference solution coming from a boundary fitted mesh that is more refined than the one depicted in Fig. 2.32a. Obviously, the strain field obtained using the proposed model has the highest smoothness compared to the other alternatives. We observe that the marching squares and the voxel based models allow to obtain similar strain fields but with more roughness. Finally, Tab. 2.1 indicates the number of degrees-of-freedom considered for each model. We emphasize that, in contrast to voxel-based approaches that are constrained by the image resolution, the immersed approach has the advantage of uncoupling the discretization of the solution space from that of the geometry, and therefore allows to keep the same B-spline approximation space for different image resolutions.

Cubic finite cell method (Fig.

An open-cell foam

As a second example, we consider a compression test on a small region extracted from a reconstructed X-ray µ-CT scan image of a polyurethane foam. This specimen will be later studied in details in the second part of this thesis where DVC measurements are performed (see chapter 5). We only illustrate here the application of the developed image-based model in three dimensions. We consider a linear elastic simulation in which we compress this foam sample. To do so, we fix the bottom surface (by fixing the nodes at minimal y height) and apply a Dirichlet displacement u y = u 0 = -1 voxels.

The considered finite element mesh has an approximate face size equal to 0.5 voxels (illustrated in Fig. 2.37b). We set some arbitrary mechanical properties to this sample by considering the Young's modulus E = 7 × 10 10 and ν = 0.3. A cubic B-spline discretization with an element size of 5 voxels 3 is considered and illustrated in Fig. 2.37c. These simulations have been run on a 64-bit laptop, equiped with Intel Core i7-3630QM processor (2.4GHz frequency, 8 CPUs), 15 GB of RAM. We compare in Table 2.2 the different parameters and results of the simulation. Figs. 2.38-2.39-2.40-2.41 the displacement and Von Mises strain fields. Qualitatively, the FCM seems to be more regular and smooth for representing strains.

First, the same observation can be made concerning the number of degrees of freedom of the finite cell method. Even though the FCM has less degrees of freedom, the resulting matrix has much more non-zero values. We remark that the direct resolution of a very large finite element system with 478 776 degrees of freedom is faster than the resolution of a small linear system obtained with cubic B-spline basis functions. The LU decomposition seems to be slower in this case. This experimental numerical test shows the limits of our implementation of the method and adapted numerical solving of the final linear system must be considered. Practically, as we are limited by the memory storage, we are unable to consider large regions of interests which is one of the goals of the second part treating DVC. In order to improve the numerical solving, multiple contributions exist in the litterature. We

Conclusion

In the first part of this work, a general framework for continuously representing images as spline functions is considered for two and three-dimensional images. The obtained smooth representation allows to construct level-set functions that can be either used for mesh-based finite elements or for immersed methods. The constructued smooth representations allow to preserve the initial image dynamic in addition to smooth image gradients which is desirable for image correlation.

We attempt to shed a new quantitative light on the construction of proper unfitted image-based models. More precisely, a special care is taken to analyze the numerical approximation error while taking into account the intrinsic geometry error resulting from the image generation process (sampling + quantization). This study thus breaks with the usual convergence analysis where the modeling error and the convergence of numerical methods are most of the time assessed with geometries deemed to be exact, which is not the case with image-based models in general. By means of various two-dimensional numerical experiments, we show that, in addition to being dependent on the modeling parameters (element size, polynomial degrees, quadrature rule), the numerical approximation error is related to the pixel size which controls the intrinsic geometry error. More specifically, the geometry error propagates into the mechanical model error. In this context, we propose rules that allow for automatic determination of the discretization parameters. For the quadrature rule in particular, the number of integration sub-cells is adjusted according to the definition of the images in order to achieve asymptotic precision with a minimal computational effort. Then the polynomial degree and element size can be properly chosen to reach the best possible accuracy while avoiding excessive and unnecessary calculations. In particular, it is shown that we do not need to go to higher orders than p = 3, while it is often stated in the literature (see, e.g., [Düster et al. 2012,Schillinger & Ruess 2015]). In that sense, the resulting image-based model is characterized as a fairly-priced model. The model constructed in this way is far more accurate with up to an order of magnitude less number of unknowns compared to all other classical (boundary fitted) approaches considered (Marching squares, voxel-based with segmentation or with mechanical properties related to grey-levels).

The term fairly priced does not mean however more computationaly efficient. In fact we have shown, using the three-dimensional test on a foam specimen, that the linear systems resulting from the finite cell method are very costly to solve with a direct solver because the sparsity of the obtained matrices decreases. As no further numerical developments were made concerning the linear system resolution aspect, we make use of the finite cell method for relatively small problems (in the limit of the provided memory storage) and for two-dimensional applications. When considering large regions of interest, we rather use the finite element method. This choice is purely practical as the main goal of this work is to perform measurement on real cellular architectures.

Digital image correlation (DIC) or image registration in the field of applied mathematics and computer vision is a general computational method for finding unknown kinematic transformations between different image configurations. It is considered as non-contacting method in the field of experimental mechanics [Schreier et al. 2009] and is now largely used to measure displacement fields on patterned surfaces (Stereo DIC [Fouque et al. 2021a, Fouque et al. 2021b, Fouque 2022]), reconstructed tomographic volumes (DVC) [Bay 2008, Gates et al. 2011, Leclerc et al. 2012[START_REF] Buljac | Digital volume correlation: review of progress and challenges[END_REF] or even in fluids containing particles. For fluid measurements, one of the most famous methods is particle image velocimetry (PIV) [Elsinga et al. 2006]. In DVC, the most common approach in the experimental mechanics community is the subset based approach [START_REF] Bay | [END_REF]. It undertakes to solve a collection of independent image registration problems in sub-volumes (the so-called subsets) distributed in the region of interest (ROI). In each subset, the displacement is usually assumed to be a polynomial function. Although widely used to study the targeted materials because of its simplicity and its computational efficiency, this approach may not appear relevant if the objective is to measure very complex and localized kinematics. As a result, a finite-element approach to DVC (often referred to as global ) was later introduced by [Roux et al. 2008] to conveniently bridge the experimental and simulation outputs. In this variant, the images are registered by using a FE displacement field defined in the whole volume. This approach provides access to globally continuous displacement fields, which can be naturally compared to the simulation results. A fine mesh would then be necessary to capture local kinematics when dealing with complex micro-structures, which would make the DVC problem ill conditioned since the kinematic basis becomes too rich with respect to the amount of (relevant) data provided by the imaging tool. However, this formulation of DVC can be easily regularized using Tikhonov like methods [Tikhonov 1963, Tikhonov & Arsenin 1977, Horn & Schunck 1981, Poggio et al. 1985] to compensate for a texture deficit under a given scale [Leclerc et al. 2011]. The idea is to add to the least-square correlation functional a regularization term introducing a new a priori on the displacement field [Poggio et al. 1985]. This can consist in a "smoothness constraint" [Horn & Schunck 1981] or even a mechanical constraint [Bajcsy & Kovačič 1989, Ferrant et al. 1999, Kybic & Unser 2003, Ruhnau et al. 2007, Réthoré et al. 2009, Leclerc et al. 2011]. Recent contributions have for instance overcome the texture issue in trabecular bone samples [van Dijk et al. 2019] or wood cell materials [START_REF] Patera | [END_REF] using regularization schemes based mostly on second-order Tikhonov differential operators.

The method established in this work extends the concept of the equilibrium gap method [Réthoré et al. 2009] to make it a true technology to measure sub-cellular, possibly inelastic, kinematic fields in speckle-free cellular materials. From a global point of view, we seek to exploit the information coming from the movement of cell boundaries (i.e., the only significant gradients available) and rely on a mechanical model to estimate displacements in the textureless micro-structure, e.g., in walls and/or struts.

From the general point of view of DIC, Fig. 2.42 illustrates the main ideas of our contribution. Given two images of the material captured during an in situ test, one in the reference state and the other once the sample has been loaded, we seek to identify the displacement field at a sub-cellular scale by DIC. To compensate for the poor texture, we suggest building an automated image-based model to regularize the ill posed DIC problem.

Mechanical test

Reference con✁guration

Deformed con✁guration Image This second part of the thesis will be organized as follows:

• Chapter 3:

In this chapter, we present a synthesis and overview of Digital Image Correlation in its two and three-dimensional versions. The fundamental optimization aspects are synthetized and a litterature of the different regularization schemes is presented. A fundamental aspect in DIC is the initialization of the optimization algorithm. This step becomes very challenging in the situation of real experimental DVC problems especially when large loading increments exist. We present the different numerical approaches that were developed in this work to tackle this important step. As our regularization model is maily based on the equilibrium gap method, we review the different numerical aspects of this approach in addition to its filtering properties following the works of [Réthoré et al. 2009, Leclerc et al. 2011].

• Chapter 4:

In order to assess the performance of the developped methodology, this chapter treats a twodimensional artificial cellular architecture (see Fig. 2.42 right). This geometry mimics a slice of a cellular foam as the one of Fig. 2 (in the general introduction). As the mechanical response of this architecture can be complex and local, the validation of the DIC method must be performed using general mechanical displacement fields, thus including non-rigid transformations. For this reason, the suggested DIC validation method consists in generating synthetic images from FE simulations and comparing the measured displacement fields to the FE reference displacement field.

The main contribution in this chapter consists in using automatic image-based mechanical modeling to feed the equilibrium gap regularization and in carrying out a L-curve study [Hansen 2000]. It allows to select the appropriate regularization weights and to interpret them as physical lengths. The mechanical regularization introduces a separation of scales above which DIC is dominant and below which it is assisted with image-based modeling. Regarding the validation, the novelty consists in assessing the method on synthetic linear elastic well as experimental non-linear elasto-plastic, possibly with local buckling, cases. Overall, we end up with a measurement method at the scale of the architecture (using the highest possible spatial measurement resolution) and basing it only on the texture of the sample. As it is based on the use of a regularization model representative of the micro-architecture of the material, we called this method Architecture-Driven Digital Image Correlation Technique (ADDICT).

Although this method could apply to non-linear regularization models [Réthoré et al. 2013], the model used here for weak regularization is elastic. The efficiency of the method to estimate local strain fields of samples undergoing possibly non-linear mechanical behaviours is analyzed considering three regimes (elasticity, elasto-plasticity and geometric non-linearity) for the generation of the synthetic images.

The Tikhonov-like terms used for the regularization of the DIC problem introduce two parameters that are trade-offs between data fidelity and regularity. This is where the L-curve study comes into play. Additionally, the influence of the regularization parameters on the true measurement error is performed. Finally, an experimental validation is performed by comparing the results of the proposed method on low resolution speckle-free images with those of a classic DIC on speckled high resolution images.

• Chapter 5:

This chapter treats the main purpose of this thesis: it constitutes the application of the above formalism and implementation to treat real three-dimensional image data. We generalize all the aspects investigated before to three-dimensions and we apply the elastic equilibrium gap regularization to DVC measurements on an open-cell polyurethane foam. We demonstrate the ability of the method to measure complex and localized phenomena, which is not possible using classical DVC approaches.

Chapter 3

Digital Image Correlation: optimization, discretization, regularization and initialization

In this chapter, the fundamentals of image correlation are presented along with the mechanical equilibrium gap regularization approach.

Given two grey-level images I r and I d representing the reference and the deformed configuration of a material sample, volume-based DIC can formally be expressed as the following:

Find u : Ω ⊂ R d → R d , such that I r (x) = I d (x + u(x)), (3.1)
where Ω is the region of interest defined in the initial volume. x defines the point coordinates in the reference region of interest Ω. d is typically equal to 2 for fronto-parallel DIC and 3 for DVC. It was underligned in [Leclerc et al. 2015] that the formulation (3.1) is a Lagrangian description of the optical flow while the formulation which consists in defining x in the deformed configuration is an Eulerian one. The Lagrangian description consists in following each point living in the reference configuration. By definition, the problem of registering two images has more unkowns than data at the voxel scale because an image is a scalar piece-wise constant function and the unknown displacement field is a vector field. The existence and uniqueness of the displacement is not guaranteed which makes the problem ill-posed. A general way for solving this problem is to use the minimization of an alignement metric as follows: arg min

u∈[L 2 (Ω)] d S(u, I r , I d ). (3.2)
Under the assumption of the conservation of grey-levels (see optical flow problem [Horn & Schunck 1981]), S can be chosen as the squared L 2 norm of the residual of grey-levels which is equal to I r -I d • (id + u), where • refers to the composition operator between two functions and id to the identity function. In this case, S is defined by [Lucas et al. 1981]:

S(u, I r , I d ) = 1 2 ||I r -I d • (id + u)|| 2 L 2 (Ω) = 1 2 Ω (I r (x) -I d (x + u(x))) 2 dx. (3.3) Problem (3.
3) is rather solved for Ĩr and Ĩd the two continuous representations of I r and I d (see section 3.4 for the definition). The L 2 distance is also known as the sum of squared differences metric (SSD) and it seems to be the most common method when the grey-level intensity is preserved or at worst varies homogeneously and linearly. Regardless of the choice of the window of alignement, multiple 3.1. Deforming images using arbitrary displacement fields other metrics and methods exist in the litterature and depend of the type of applications. Multiple examples of similarity measures were listed for example in [Modersitzki 2009[START_REF] Ozeré | Mathematical modelling of problems related to image registration[END_REF]]. We list here some of them:

• Linear correlation coefficient: consists in maximizing the correlation factor (or its equivalent named cross-correlation) [Lucas et al. 1981, Bajcsy & Kovačič 1989, Périé et al. 2002].

• Mutual information: Mutual information is a similarity measure that comes from information theory. It allows to quantify how much the information of the image I r is present in the image I d and vice-versa. It can be used in image correlation when the images are acquired with different imaging modalities (i.e. when grey-level conservation is not verified). Mutual information relies mainly on probabilistic concepts such as the joint histogram and the entropy of images [START_REF] Viola | Alignment by maximization of mutual information[END_REF], Rueckert et al. 1999].

• Gradient fields: This method consists in aligning the gradients of I r and I d . One can choose the norm of the cross product of the gradients of the two images [Haber & Modersitzki 2006a] in order to quantify the alignement. This approach can be also be used for level-set registration as the image gradient is orthogonal to the level-sets.

• Joint image segmentation and registration: Combines energy functionals in which both the alignement of grey-levels and the determination of the unknown contours appear [Yezzi et al. 2001].

We decide to choose the sum of squared distance functional because it is the most robust one when considering images acquired with the same modality. It also offers an elegent way for discretization using the same discretization tools for solving mechanical problems (this will be detailed in section 3.3).

Deforming images using arbitrary displacement fields

We present in this section how we can deform an image using a reference displacement field. This problem can be viewed as the direct problem of DIC. Better understanding this direct approach can help understanding the inverse problem. In addition, deforming images with reference displacement fields can be very helpful for the design of DIC algorithms as it allows to perform uncertainty analysis and virtual numerical experimentations.

Given an image I r and a continuous displacement field u : R 2 -→ R 2 , let us construct an image I d that verifies the grey-level transportation condition (3.1). The idea is to assign a grey-level value to each pixel (or voxel) in the grid of the image of the deformed configuration. For more clarity we denote respectively Ω r and Ω d the domains of the reference and deformed configurations. Given x d ∈ Ω d , we look for the grey-level that must be assigned to this point. This problem reads:

find x r ∈ Ω r such that x d = φ(x r ) = x r + u(x r ).
(3.4)

Deforming an image consists in inverting φ for a given displacement u. We note that the inverse transformation of φ is not reduced to the substratction of u(x d ) from a point x d in the deformed configuration. A hypothesis that makes the DIC problem solvable and that allows us to deform the image is that φ is a diffeomorphic map which means that φ is differentiable and admits an inverse map φ -1 that is also differentiable. This means that the determinant of the Jacobian matrix of φ is non zero (i.e. the initial topology is preserved). One can invert the mapping φ using a Newton algorithm:

x r k+1 = x r k + J -1 φ (x r k )(x d -φ(x k r )) with J φ (x) = I + J u (x) (3.5)
with an initial start x r 0 = x du(x d ) for example. Once the coordinate x r is identified, a value of I r (x r ) is assigned to I d (x d ) with an interpolation method (typically a piecewise polynomial spline interpolation function or a spline L 2 projection). When the transformation φ is affine (i.e. of form φ(x) = Ax + b which includes shear, scaling, translation and rotation), the image of the deformed configuration can be obtained by inverse Fourier Transform without resorting to interpolation. As

I d (x d ) = I r (φ -1 (x d
)), we have by change of variables the following result (coming from the Fourier dilation and translation theorems [Olson et al. 2017]):

F(I d )(ξ) = F(I r • φ -1 )(ξ) = R d I r (φ -1 (x))e -i2πξ T x dx = |det(A)|e -i2πξ T b F(I r )(A T ξ) (by change of variable)
.

(3.6) Formula (3.6) shows that in order to deform I r , we only need to perform a linear transformation in the frequency domain and multiply the fourier transform by the phase and the factor |det(A)|. The image I d is finally obtained using the inverse Fourier transform of (3.6). Depending on the type of transformations used to deforme the image I r , one can either choose to use the Fourier transform or the Newton method. In the context of this work, general finite element displacement fields are interpolated from the mesh to the image domain therefore we rather choose to use the Newton method.

In addition, as we will deform images of geometric objects (typically finite element meshes) of cellular architectures, we will also use a simple rasterization algorithm that creates super-resolved binary images of the reference mesh and the deformed mesh then downsamples them to desired resolution by image binning. We have to note that this approach reaches rapidly its limits in three-dimensions as full volumetric images become very memory consuming.

Resolution of the correlation problem

There are multiple manners of viewing the differentiation of the optimization problem (3.2)-(3.3). One can discretize directly the functional as a non-linear least squares problem and then use the classical formulas of discrete optimization schemes [Passieux & Bouclier 2019]. Instead, we keep the continuous differential forms so that we transit easily to the finite element formalism. The functional S in (3.3) is twice differentiable at u provided the regularity of I r and I d . By this definition, there exist a linear form S ′ (u) and a bi-linear form S ′′ (u) that verify the expansion (3.7). As the displacement is a vector field, we note

L 2 (Ω) = [L 2 (Ω)] d
, where L 2 is the space of square-integrable functions. L 2 (Ω) is equiped with the dot product defined as f, g

L 2 (Ω) = Ω f (x) T g(x)dx. S(u + v) = S(u) + S ′ (u)(v) + 1 2 S ′′ (u)(v, v) + •(||v|| 2 L 2 (Ω) ) = S(u) + ∇S(u), v L 2 (Ω) + 1 2 H S (u)v, v L 2 (Ω) + •(||v|| 2 L 2 (Ω) ) , ∀u, v ∈ L 2 (Ω). (3.7)
By applying the expansion (3.7) to the functional (3.3) then developping a first order Taylor expansion to the image I d and eliminating the terms that are least quadratic, we show that the first differential of S is defined as follows:

S ′ (u)(v) = Ω -(I r (x) -I d (x + u(x))) v(x) T (∇I d )(x + u(x))dx, ∀u, v ∈ L 2 (Ω). (3.8)
Using the dot product in L 2 (Ω), we obtain the following form for the gradient of S:

∇S(u) = -(I r -I d • (id + u))(∇I d ) • (id + u), ∀u ∈ L 2 (Ω).
(3.9)

Solving the DIC problem consists in finding the optimal displacement u that vanishes ∇S(u). Therefore the equation

f (x, u(x)) = -(I r (x) -I d (x + u(x)))(∇I d )(x + u(x)) = 0, ∀x ∈ Ω (3.10)
provides necessary conditions for u being a minimizer of (3.2). The authors in [Fischer & Modersitzki 2004] have defined it as an Euler-Lagrange equation in which f (x, u) is simply the driving registration force. Once this gradient determined, one can directly discretize the problem by defining a discrete displacement defined at voxels or windows and use directly the Newton method for solving a nonlinear equation of form F(u) = 0 where u is a vector that collects the degrees of freedom. We rather present in this section a continuous approach for solving the optimality condition by considering the fixed point algorithm on the Euler-Lagrange equation (3.10). It considers a sequence (u k ) k that converges to the solution of (3.10):

u k+1 = u k + v k with u 0 ∈ L 2 (Ω) an initial displacement guess. (3.11)
v k is the unique solution of the Linearized Euler-Lagrange equation:

J f (u k )v k + f (u k ) = 0 (3.12)
which final variational form is given by performing a Galerking projection:

J f (u k )v k + f (u k ), δv L 2 (Ω) = 0, ∀δv ∈ L 2 (Ω). (3.13)
In optimization, this scheme can be viewed as the direction of descent v k of the Newton algorithm and this is defined as the unique minimizer of the quadratic approximation of S [Passieux & Bouclier 2019]:

v k = arg min v∈L(Ω) D(v) = S(u k ) + S ′ (u k )(v) + 1 2 S ′′ (u k )(v, v) = S(u k ) + ∇S(u k ), v L 2 (Ω) + 1 2 H S (u k )v, v L 2 (Ω) .
(3.14)

The optimality condition of this problem yields

D ′ (v k )(δv) = 0, ∀δv ∈ L 2 (Ω) (3.15) with D ′ (v k )(δv) = H S (u k )v k + ∇S(u k ), δv L 2 (Ω) (3.16)
In fact, as J f = J(∇S) = H S , we obtain the same form as (3.13). By applying a second Taylor expansion to ∇S, we can show that the Hessian of S is defined by :

H S (u) = (∇I d ) • (id + u)(∇I d ) • (id + u) T -(I r -I d • (id + u))(H I d ) • (id + u) = (∇I d ) • (id + u) ⊗ (∇I d ) • (id + u) -(I r -I d • (id + u))(H I d ) • (id + u), (3.17)
where H I d is the Hessian matrix of I d defined by

(H I d ) i,j = ∂ 2 I d ∂x i ∂x j . (3.18)
The complete fixed point iteration is finally given by:

δv, H S (u k )v k + ∇S(u k ) L 2 (Ω) = 0, ∀δv ∈ L 2 (Ω). (3.19)
From a computational point of vue, two major approximations are considered in this work:

• Hessian approximation: only the linearized Hessian is considered by neglegting the image Hessian term. This algorithm is known as the Gauss-Newton method [Passieux & Bouclier 2019]: To summarise, the optimization algorithm can be written as follows:

H S (u) = (∇I d ) • (id + u) ⊗ (∇I d ) • (id + u) (3.20) • Image gradient approximation: (∇I d ) • (id + u) is
1. Determine an initial guess u 0 .

2. For each iteration k + 1:

Solve δv, H S (u k )v k + ∇S(u k ) L 2 (Ω) = 0, ∀δv ∈ L 2 (Ω) (3.23)
which is equivalent to solving the variational form The form a is symmetric bilinear and the form l is linear. The existence and unicity of the solution of problem (3.24) can be proved using Lax-Milgram theory (see for example [Fedele et al. 2013b]) commonly used for the analysis of PDEs. We note however that this problem is different as it is strongly dependant of the data used (namely the images I r and I d ). In general, the functional S is not convex and no guarantee of a unique minimizer exists. Therefore, solving the problem voxel wise is in general impossible using the previous scheme as simply the tensort product ∇I r ⊗ ∇I r is equal to the zero matrix. That is why, a window of correlation is generaly used in regularized or non-regularized methods. We will introduce in the next section the discretization of the correlation problem.

a(v k , δv) = l(u k , δv), ∀δv ∈ L 2 (Ω) (3.24) with a(v k , δv) = Ω v k (x) T ∇I r (x) ⊗ ∇I r (x)δv(x)dx (3.25) and l(u k , δv) = Ω (I r (x) -I d (x + u k (x))) ∇I r (x) T δv(x)dx. ( 3

Discretization

The discretization methods of the DIC problem are very diverse and each community have been using different implementations and manners for viewing this problem (whether in the experimental mechanics community or the bio-medical imaging community). We list here some discretization techniques. As the litterature is very rich and diverse, we have tried to do our best to find the oldest references:

• Finite differences:

The idea is to consider the uknown displacements as the voxel displacements (as it was originally defined in the optical flow problem [Horn & Schunck 1981]) and consider finite-differences schemes on rectangular or cuboid grids. As the problem is in general ill-posed at the voxel scale, these methods are in general assisted by regularization [Horn & Schunck 1981].

• Sub-set methods:

This class of methods consist in solving the DIC problem through voxel summation over windows called sub-sets [Sutton et al. 1983, Chu et al. 1985]. Each independent subset transformation is supposed affine [Roma et al. 2002]. With this technique, the resolution can be very efficient as the correlation problem is solved independently on each subset. The numerical efficiency of this method can be even more ameliorated using the Fourier transform when the cross-correlation is maximized [Périé et al. 2002]. The only costly step in these methods concerns the initial displacement guess estimation that is usually performed through voxel search. If very large increments are considered, starting from a zero displacement vector can lead sometimes to very long runtimes especially in DVC.

• Interpolation methods:

The general idea of interpolation methods is to solve the problem of Box. 3.1 by expressing the displacement field in a basis of functions exactly as it would be done for discretizing a PDE with a Galerkin approach. Roughly speaking, the displacement field is expressed by a discretization:

u(x) = n i=1 N i (x)u i (3.27)
where (N i ) i∈ 1,n is the total set of basis functions and u i is the correspondant degree of freedom.

There is a very large spectrum of approaches for choosing the discretization basis. We list for example the spectral approaches that consists in choosing the Fourier basis [Roux et al. 2002,Wagne et al. 2002]. In this work, we consider the two discretizations that were presented in the first part, that is to say, the finite element discretization [Ferrant et al. 1999,Sun et al. 2005,Besnard et al. 2006] using triangular and tetrahedral finite elements and the B-spline discretization of arbitrary degrees [Rueckert et al. 1999, Réthoré et al. 2010[START_REF] Kleinendorst | On the use of adaptive refinement in isogeometric digital image correlation[END_REF], Dufour et al. 2016, Patera et al. 2018] to name few. Both methods are equivalent and will be used for assessing the DIC algorithms. Typically, when considering a structured linear B-spline grid, it is strictly equivalent to structured bilinear (Q 4 ) or trilinear finite element meshes (C 8 ). The advantage of B-splines is that their higher regularity offers a better description of strains especially if coarse elements are used.

We go back to formula (3.27). In practice, as the shape functions have a compact support, the discretization is usually defined element wise using a local basis function operator N ∈ R d×dne where n e is the number of basis functions per element (see Eq. (A.39) in appendix A for the details of a classical finite element discretization). A Galerkin approach is applied to Eq. (3.24) by performing summation over elements and adding block contributions of elementary matrices in the same fashion as any finite element assembly algorithm. In this case, the unknown displacement vector u k is updated using the descent scheme of Box. 3.1:

u k+1 = u k + d k (3.28) with Hd k = b k , H ∈ R dn×dn , b k ∈ R dn (3.29)
where d k ∈ R dn is the displacement increment. n is the total number of nodes or control points. H is obtained from the assembly of elemental block matrices H e ∈ R dne×dne defined by: Given assumptions (3.20) and (3.21), the operator H is constant at each iteration. As it is symmetric positive definite, it is factorized using the Cholesky method in a pre-processing step. Only the righthand side is re-assembled during the descent iterations. In the computational code, the assembly of the right hand side is assembled in parallel (using shared memory parallelism) in order to speed up the iterations. The matrices H and b can be viewed directly as the Hessian and the gradient of the discretization of the functional S.

H e = Ωe N (
The integration of the terms in (3.30)-(3.31) is performed using a rectangle integration with a subdivision volume proportional to the element resolution, i.e. a number of integration points in one element greater or equal to the number of its voxels. Therefore the integrals are defined as folllows:

Ω f (x)dx = n i=1 ωf (x i ) (3.32)

Sub-voxel evaluation

where the set (x i ) i∈ 1,n is the uniform subdivision of the element, f is the integrand and ω is the constant integration weight. In three dimensions, a uniformly subdivided grid of points (x i , y j , z k ) i,j,k with i ∈ 1, n x , j ∈ 1, n y , k ∈ 1, n z is considered where

n x = α e x ∆x
, n y = α e y ∆y , n z = α e z ∆z (3.33) and e x , e y , e z are the univariate lengths of the B-spline element. α is chosen at least equal to 1 (in this case it is equivalent to a voxel summation). As the regularity of I r is increased through the continuous representation (usually cubic), the greater the value α is chosen the higher accuracy of summation is obtained. Fig. 3.1a illustrates the uniform subdivision of a structured B-spline element. In this case, the integration weight is set to :

ω = V n x × n y × n z , (3.34)
where V is the constant element's volume. When a tetrahedral element is considered, the uniformly subdivided integration grid is considered by setting

n x = n y = n z = α V ∆x∆y∆z 1/3 , (3.35)
where V is the average volume element of the finite element mesh (if the meshing is uniform). We then only keep the points that verify z i ≤ 1x iy i , see Fig. 3.1b. In this case, the integration weight is set to:

ω = 1 6n x . (3.36)
Finally, the iterative scheme (3.28) is stopped using a stagnation criteria defined by

||u k+1 -u k || 2 ||u k+1 || 2 < ε.
(3.37)

Sub-voxel evaluation

Sub-voxel evaluation is an important aspect in DIC as it can influence the measurement accuracy [START_REF] Schreier | [END_REF]. In this work, we make use of the cubic B-spline representation that was constructed in chapter 1 section 1.2.1:

Ĩ(x) = m i=1 c i N i,p (x). (3.38)
We therefore correlate Ĩr and Ĩd the continuous representations of I r and I d . The B-spline coefficients c k are obtained by the filtering of the images I r and I d with a Gaussian filter of standard deviation σ ≤ 1. This improves the sub-voxel accuracy as it allows to remove the image noise. It is also a type of regularization acting on the data. In the case of our images of cellular materials, the Gaussian filtering allows to regularize the edge effects and provides smooth image gradients. In this work, we make use of this technique without any proof or experimental evidence. The superiority of the method was only observed when confronted to other interpolation schemes such as piecewise polynomial interpolatory functions. Figure 3.1: Integration points for the structured cuboid and tetrahedral elements. A structured grid of integration points is defined for the B-spline element. For the reference tetrahedron, only the points that are under the plane 1xyz are considered.

Image contrast enhancement

In practice, the acquired images I r and I d do not have the same grey-level histograms. Grey-level variations can occur locally or globally due to multiple factors such as noise or local changes in illumination. In this respect, the images can be globally corrected by rather correlating I r and I d defined by

I r = a 1 I r + b 1 and I d = a 2 I d + b 2 . (3.39)
In the context of this work, a global contrast correction (by centering and standardization) of the images is performed in a pre-processing step. We note that this approach is only valid for global variations of contrast. When local variations occur, one need to use other local contrast correction techniques. We also note that for global brightess correction between different loading increments, the following strategy was suggested in [Genet et al. 2018]: after correlating I r and I t d , i.e. finding the optimal displacement u, the authors find the optimal parameters (a, b) that minimize S(u, I r , aI d + b) (it results into solving a 2 × 2 linear system). The obtained scaling is afterward used to correlate I r and I t+1 d . The subscript t denotes here the increment step t.

Regularization

The problem of correlating two images is ill posed due to multiple reasons. First, from a data view point, a voxel has only one information which is its grey level value but the unknown displacement for each voxel has three components. It results that there are more unknown variables than input variables at the voxel scale. Therefore the correlation problem cannot be solved at the voxel scale.

The discretization in Eq. (3.27) introduces a spatial regularization that can be characterized as a strong regularization in the sense that it is directly related to the size of the approximation subspace. Roughly speaking, to be able to solve the minimization problem, the subset or finite element size must be chosen so that the amount of grey-level data available in a subset or finite element is richer than the corresponding elementary kinematic basis. In the conventionally used subset-DIC framework, the usual rule in this respect is to set a subset size that contains at least 3 texture dots [Schreier et al. 2009,Bornert et al. 2009,Jones & Iadicola 2018]. When considering cellular materials, this would lead to a subset size as depicted in Fig. 3.2 (see also section 4.2.2 where further details regarding this image are provided). At best the subset size would be equal to the cell size. Obviously, the resulting approximation space appears too coarse in view of estimating the kinematic fields at the sub-cellular scale. A finite element mesh as fine as the one of Fig. 3.2 would be necessary instead but, in this case, the strong regularization would not be sufficient anymore, thereby leading to a singular matrix H in Eq. (3.29). As a first solution, one could use, for instance, the Levenberg-Marquardt algorithm with an update strategy for stabilizing the Hessian matrix [Szeliski & Lavallée 1996]. However, this lacks of physical meaning and thus will not provide relevant derivative fields (such as strains and stresses) within the strut. Figure 3.2: Size of subset (red rectangle) to properly regularize the DIC problem coming from images of speckle-free cellular type materials. The resulting approximation space appears too coarse in view of estimating the kinematic fields at the sub-cellular scale. A finite element mesh as fine as the one depicted in this figure would be necessary instead, thus leading to a severely ill-posed inverse problem.

In addition to the treatment based on the Galerkin approximation (i.e., given by the approximation space related to the basis functions gathered in N) adding an energy term R in (3.2) allows to convexify the alignement distance when the grey-level distribution is not sufficient at the measurement scale, i.e. when the image gradient is either small or has a dominant direction, which is the case of cellular materials. This yields to the optimization form given in (3.40). Such a term also allows to softly introduce a prior knowledge on the physics behind the deformation of the sample.

The general form of the regularized correlation problem is given by the following optimization problem: arg min

u∈L 2 (Ω) S(u, I r , I d ) + λ 2 R(u). (3.40)
Multiple approaches can be adopted in order to choose R. In a discrete way, many forms of regularization can be written using the Tikhonov approach:

arg min

u∈R n S(u) + λ 2 u T Tu (3.41)
where T is a linear operator. The finite element discretization offers a great flexibility as it allows to directly construct differential operators coming from physical models. The final form of the descent step for the regularized optimization problem is given by:

(H + λT) d k = b k -λTu k . (3.42)
In the following, a small overview of the different regularization models is presented.

Regularization models

A very rich survey concerning bio-medical registration is given in [START_REF] Sotiras | Aristeidis Sotiras. Discrete image registration: a hybrid paradigm[END_REF]]. The authors in [Fischer & Modersitzki 2004] also list a large panoply of regularization energies in a continuous framework. In these two references, registration is generally presented by an equilibrium equation in which the volume forces f drive the registration (see Eq. (3.10)) and the equilibrium is completed by a physical operator A. We have the following equivalence:

f (x, u(x)) + λA(u)(x) = 0, ∀x ∈ Ω (Euler-Lagrange) ⇔ arg min u∈L 2 (Ω) S(u) + λ 2 R(u) (Energy) , (3.43)
where f and S are defined respectively by Eq. (3.10) and Eq. (3.3). We list some methods for choosing the regularization model:

• Diffusion model [Horn & Schunck 1981] -Euler-Lagrane:

A(u) = ∆u (3.44)
-Energy:

R(u) = Ω ||J u (x)|| 2 F dx = d i=1 Ω ||∇u i (x)|| 2 2 dx. (3.45)
In this case, the correspondant Tikhonov operator T is equal to L the Laplacian operator defined at the mesh nodes. It is a second order regularization that penalizes large variations of the displacement.

• Curvature model [START_REF] Fischer | Curvature based image registration[END_REF] -Euler-Lagrane:

A(u) = ∆ 2 u = ∆(∆u) (3.46)
3.6. Regularization

76

-Energy:

R(u) = Ω d i=1 (∆u i (x)) 2 dx. (3.47)
This fourth order regularization needs at least a quadratic discretization mesh with a C 1 regularity. However, it can be well approximated by setting T = L T L when considering a linear finite element mesh. This approximattion was first used in the DIC community in [Leclerc et al. 2011]. Both authors in [START_REF] Fischer | Curvature based image registration[END_REF]] and [Leclerc et al. 2011] have highlighted the fact that affine transformations belong to the kernel of this regularization and therefore are not penalized which does not allter the measurement of coarse (or global) transformations.

• Quadratic variation model [Rueckert et al. 1999]:

-Energy:

R(u) = Ω d i=1 ||H u i (x)|| 2 F . (3.48)
This energy also known as the Thin-plate energy represents the bending energy of a plate (when d = 2) . This regularization was considered for example in [START_REF] Patera | [END_REF], van Dijk et al. 2019] for DVC. This energy takes into account the cross partial derivatives whereas the curvature model takes into account only the square of direct second order partial derivatives.

• Elastic model [Bajcsy & Kovačič 1989, Ferrant et al. 1999]:

-Euler-Lagrane

A(u) = ∇.σ(u) = ∇. E 2(1 + ν) ε + E (1 + ν)(1 -2ν) tr(ε)I = E 2(1 + ν) ∆u + E 2(1 + ν)(1 -2ν) ∇(∇.u), as ε(u) = 1 2 (∇u + ∇u T ) (3.49)
where E and ν are respectively the Young's modulus and Poisson coefficient of the material. ∇. stands for the divergence and σ and ε represent respectively the stress and strain tensors.

-Energy:

R(u) = 1 2 Ω σ(u), ε(u) F dx - Γn u T f s ds, with σ(u)n = f s on Γ n (3.50)
where Γ n are the boundaries that are not free (Neumann boundaries).

When discretized, the algebraic form of the energy is equal to R(u) = 1 2 u T Kuu T f , where K is the stiffness matrix and f is the external force vector.

• Elastic equilibrium gap: Instead of considering the elastic energy, one can further increase the order of the regularization in an analogue way of transitioning from the diffusion model to the curvature model. The idea is to consider the squared norm of the equilibrium gap. The equilibrium gap seems to be written only in a discrete form. In this case, the discrete energy is defined by the squared norm of Kuf :

R(u) = ||Ku -f || 2 2 .
(3.51)

This approach was initiated in [Claire et al. 2004] for damage identification and was later introduced as a regularization in [Réthoré et al. 2009] where it was used for the identification of cracks in a silicon carbide specimen. It was later detailed in [Leclerc et al. 2011] where a physical interpretation of the regularization parameter λ was established. The advantage of the equilibrium gap regularization is that it is a fourth order filter so it is a less strong regularization compared to the pure elastic model of Eq. (3.49). Until this date, the equilibrium gap regularization have never been used for the measurement in cellular materials at the micro scale. This motivates its application in this work. It can be therefore applied in order to measure general transformations that are not necesserily elastic. This will be further investigated in chapter 4.

Practically, when acquiring images of deformed samples in solid mechanics, f can be reasonably assumed to be zero everywhere expect on the Dirichlet and non-zero Neumann boundaries. Indeed, we do not know well the distribution of the reaction forces or loading on such boundaries (we can only experimentally access to one component of the resultant force vector). That is why the equilibrium is prescribed only at the nodes of the bulk and at free boundaries. To do so, we introduce a binary selection operator D K that selects theses nodes. The equilibrium gap is written:

R K (u) = ||D K Ku|| 2 2 .
(3.52)

As the Dirichlet nodes are removed from the regularization, they are guided only by the correlation. This can lead sometimes to a very irregular solution when the element size is very small. That is why a special treatment of the Dirichlet/non-zero Neumann boundaries (i.e. those where f is not well known) must be considered. This is detailed in the next section 3.6.2.

• Total variation

We should also highlight the fact that all the previous listed regularization models can be effective for measuring smooth and regular displacement fields. In the case of irregular displacement fields (i.e. with steep gradients), strategies based on total variation regularization were suggested in the litterature [Frohn-Schauf et al. 2008, Chumchob et al. 2011, Vishnevskiy et al. 2016] to treat this case. The total variation of the displacement field is of the same form as the diffusion model but without the square of the norm:

R(u) = T V 1 β (u) = d i=1 Ω |∇u i | β dx, with |∇u i | β = d j=1 ∂u i ∂x j 2 + β (3.53)
where β is a very small parameter introduced for avoiding singularities when ||∇u i || 2 = 0. The equilibrium corresponding to the total variation energy is defined as follows:

A(u) i = ∇. ∇u i |∇u i | β . (3.54)
All these regularization schemes can be viewed as a certain filter that acts on the measured displacement so that large strain variations are attenuated. In this work, we make the choice to investigate and use the equilibirum gap approach (3.51) for regularizing DIC and DVC as it is justified by the mechanical behavior of the observed material. In addition, the highest the order of the regularization scheme, the less a priori is put, which is prefered if we want to use the measured displacements for identification purposes.

Boundary conditions for the equilibrium gap regularization

As the equilibrium gap regularization is used only on the degrees of freedom located at the bulk and the free edges (or surfaces), a second regularization term is introduced for the Dirichlet boundaries.

To do so, the complementary of the nodes selected by D K are also penalized using the model (3.47).

In , the red nodes that belong to non-free surfaces or edges are selected by the operator D L . In total, two regularization terms are considered: R K (u) and R L (u). In the litterature, two boundary stabilization strategies were considered:

• Boundary restriction of the volume Laplacian operator [Leclerc et al. 2011]:

R L (u) = ||D L Lu|| 2 2 (3.55)
where L is the volume Laplacian operator.

• Boundary restriction of the surface Laplacian operator [Mendoza et al. 2019]: Consists in penalizing each boundary (i.e. a union of the tetrahedral faces in three dimensions or the tiangles edges in two dimensions). For each Dirichlet boundary Γ d i an edge regularization term is considered as follows:

R L i (u) = u T (D L i K) T L i (D L i K)u (3.56)
where L i is a surface Laplacian defined on the boundary Γ d i . It is different from the volume Laplacian operator as it is defined using two-dimensional shape functions that map the two-dimensional reference element to a three-dimensional surface triangle. This regularization approach can be seen as a stabilization of the unkown external forces.

In this work, the considered optimization consists in employing the regularization (3.52) along with the stabilization (3.55). We make the hypothesis of an homogeneous and isotropic material at the sub-cellular scale with a constant Young's modulus E and Poisson ration ν. As the stiffness matrix is proportional to E, this parameter is set to 1. The final optimization scheme for the mechanically regularized image correlation is given by (3.57): arg min where λ K and λ L are the bulk and boundary regularization weights.

u∈R ndof S(u) + λ K 2 ||D K K(E = 1, ν)u|| 2 2 + λ L 2 ||D L Lu|| 2 2 , (3.57)
In our case, where the region observed represents a cellular material, the stiffness operator K is obtained by an automatic image-based mechanical model constructed in part I. The image-based model describes the elastic behavior of the cellular micro-architecture using either the FE method of the FCM method.

It has to be underlined that the (homogeneous and isotropic) elastic behavior at the sub-cellular scale is not prescribed in a strong way in Eq. (3.57). It is only used as a low pass filter to alleviate oscillatory effects [Leclerc et al. 2011, Leclerc et al. 2012]. From a global point of view, we exploit the information coming from the movement of cell boundaries (with S(u) in Eq. (3.57)) and weakly prescribe a locally elastic behavior to softly regularize DIC in the textureless microstructure, which makes sense in continuum mechanics, even for measuring inelastic fields as will be demonstrated in chapter 4. In some sense, such a procedure enables mitigating the tradeoff between the FE interpolation error (sometimes referred to as model error in DIC) and so-called ultimate random error (that is related to the ill-posedness of the inverse problem) [Bornert et al. 2009, Passieux et al. 2015]. Overall, when using this regularization, three a priori input parameters (λ K , λ L , ν) influence the DIC measurement quality. In theory, a correct estimation of ν must be provided which remains a problem for this class of methods. However it can be updated [Réthoré 2010]. The problem thus focus on the fine tuning of (λ K , λ L ), which will be addressed in chapter 4.

Filtering properties of the equilibrium gap regularization

As the different optimization residuals are not normalized in Eq. (3.57), typical values of λ L and λ K range from 10 1 to 10 9 and their sensitivity to the measured field is not constant across this interval. Besides, the link between λ L and λ K and physical lengths is not obvious. As a remedy, a mechanical interpretation of these regularization schemes has been introduced in [Leclerc et al. 2011, Leclerc et al. 2012]. To start with, a normalization of the functional in Eq. (3.57) can be considered using a reference shear wave displacement v of a wavelength T . It consists in dividing each term in the optimization in Eq. (3.57) by its evaluation at the displacement v. By denoting v the DOF vector associated to the finite element discretization of v, the descent direction using this normalization is therefore given by the following linear system:

H + λ K v T Hv ||D K Kv|| 2 2 K T D K K + λ L v T Hv ||D L Lv|| 2 2 L T D L L d (k) = b k -λ K v T Hv ||D K Kv|| 2 2 K T D K K + λ L v T Hv ||D L Lv|| 2 2 L T D L L u (k) .
(3.58)

Let us note at this stage that the left-hand side operator still remains constant and only the righthand side is updated during the optimization iterations. Using spectral analysis, it can be shown that the linear operators L and K used for regularization can be interpreted as low-pass filters (see, again, [Leclerc et al. 2011,Leclerc et al. 2012]). More precisely, regularizing using the squared L 2 norm of the second-order differential operators K and L can be seen as a fourth-order low-pass filter acting on the measured displacements on both the bulk and boundary regions. As a result, the regularization weights λ K and λ L can be related to cut-off characteristic lengths denoted l K and l L , respectively, which verify:

λ K = l K T 4 , λ L = l L T 4 . (3.59)
As λ K and λ L are dimensionless, the characteristic lengths l K and l L have the same unit as the period T of the shear wave which is in pixels (or voxels). For a proper study and a mechanical interpretation of the implemented methodology, the regularization weights will be tuned in chapter 4 by changing the values of the cut-off wave-lengths l K and l L (see section 4.3 in particular). The value of parameter T has no real influence on the results: it is just requested to take it large enough so that the wave v can be accurately described by the considered finite element or B-spline mesh.

Initialization

The Gauss-Newton algorithm is a gradient-based optimization algorithm. This induces that we can only converge to a local minimum near the initial guess, that is why a good initialization is essential in order to obtain accurate measurements. It is common to use efficient and fast strategies for finding the initial displacement field u 0 . The idea is to rapidly reach a displacement field that satisfies small values of the residual of grey-levels. Multiple approaches for the initialization exist in the litterature. For example, the maximization of the cross-correlation in the Fourier space using the Fast Fourier Transform can first be used [Périé et al. 2002, Forsberg et al. 2010] before to proceed to accurate measurements. Other matching techniques based on the optimal transport theory were also proposed in order to tackle alignement for large scale problems [Haker et al. 2004, Feydy et al. 2017, MacNeil et al. 2019]. In this work, we resort to multi-level initialization. It consists in constructing a pyramidal refinement strategy for both the image resolution and the measurement mesh. It is a general method that can be applied to a large variety of inverse problems [Ascher & Haber 2001]. For the mesh aspect, element size is gradually decreased until we attain the last desired measurement resolution [Bajcsy & Kovačič 1989,Rueckert et al. 1999,Haber & Modersitzki 2006b,Réthoré et al. 2007,Fedele et al. 2013a, Gomes Perini et al. 2014]. We also note that it is also possible to incoroporate a variation of the regularization parameter during the multi-level refinement so that the problem is adequatly convexified at each level. In this work, this approach was used but is considered more as an heuristic as no proof of the efficiency of the method was established.

Multi-level initialization

The multi-level refinement algorithm is synthetized in Algorithm 1. It is also illustrated through Fig. 3.5 for one multi-level step when using a finite element discretization.

Algorithm 1 Multi-level refinement algorithm Multi-level algorithm: → Setting the sequence of images at different resolutions (I r,k , I d,k ) k∈ 1,l → Setting the sequence of meshes (M k ) k∈ 1,l → Setting the sequence of regularization parameters (λ k ) k∈ 1,l k = 1 u i = 0 (Zero initial displacement field vector)

u k = Correlate I r,k and I d,k with (M k , u i , λ k ) while k < l do
Displacement exchange between the mesh M k and the mesh M k+1 :

u i = p M k+1 (u k ) u k+1 = Correlate I r,k+1 and I d,k+1 with (M k+1 , u i , λ k+1 ) k = k + 1 Optimal solution = u (k-1)
In order to perform the displacement exchange between different meshes, we list below the different possible configurations depending on the used mesh (structured B-spline mesh or finite element mesh). Given two meshes M 1 and M 1 and a displacement field u 1 at the nodes of M 1 , the objective is to project u 1 on M 2 . The different possibilities are listed in the diagram shown in Fig. 3.6. For a structured B-spline mesh, the location of a point is immediate. For any point belonging to the mesh domain, only a location in the univariate spaces is needed to find the three-dimensional element containing this point. It is simply performed using a one dimensional search in a ordered set. If the one-dimensional nodes are uniformly subdivided with a constant element size h then the index of the one-dimensional element containing the point x is ⌊(xx 0 )/h⌋ where x 0 the origin. Once the point is located, we just need to evaluate the displacement with the B-spline linear combination (3.27).

Isogeometric knot insertion (Case 2)

When considering a B-spline mesh, uniform element refinement can be perfomed efficiently thanks to algebraic operators called "knot insertion" operators. Given two B-spline meshes M 1 and M 2 of equal order such that M 2 ⊂ M 1 , there exists a linear operator (called knot insertion operator) and denoted can be performed using AABB trees (for Axis Aligned Bounding Boxes). As shown in Fig. 3.7, efficient trees that group the mesh componenents into boxes is constructured in order to accelerate the search operations. This allows us to detect the nearest face to the point.

• Location of the tetrahedral element:

After determining the nearest face, the location test is performed on the tetrahedrons that share this same face (they are at most two). To do so, one can consider two methods:

-Method 1: Computation of the barycentric coordinates by resolution of the linear system:

           x = m i=1 λ i t i m i=1 λ i = 1 (3.65)
If λ i > 0, ∀i ∈ 1, m then the point x belongs to the tetrahedron bounded by the nodes (t i ) i∈ 1,m . m is the number of nodes per convex set (4 in our case).

-Method 2: A faster method which does not consist in solving a linear system can be considered. We only evaluate the signed distance of the point to each of the tetrahedron faces. First, the orientation of the faces must be first determined so that all the face normal vectors point towards the same direction. This is given by an orientation matrix denoted O ∈ R 3×4 that depends of the used mesh. Each column of O represents the indices of the nodes of the terahedron faces. The point x belongs to the tetrahedron if the distances to all the faces have the same sign. This is written as follows:

φ(x) = (x -t O 2,i ) T n i (t O 2,i ) < 0, ∀i ∈ {1, .., 4} (3.66)
where n i (t O 2,i ) is the normal vector at node t O 2,i (which is the second node of the face i).

It is defined by n

i (t 2 ) = (t O 2,i -t O 1,i ) × (t O 3,i -t O 2,i ).
• Determination of the iso-parametric coordinates:

Once the tetrahedron containing the point is determined, the isoparametric transformation (x = i N i (ξ)t i ) is inverted in order to find the isoparametric coordinate of the point x. The finite element interpolation formula can afterwards be applied to evaluate the desired displacement field (u(x) = i N i (ξ)u 1 (t i )).

Reduced-basis initialization

Another interesting approach for initialization in image correlation could be the use of model reduction techniques such as the reduced basis method. The reduced basis method is a general technique that can be applied to any type of inverse or direct problem that can lead to the resolution of a linear system, in particular, in our work, the correlation linear system (3.29). To do so, we define a reduced displacement field space B (for reduced basis). Basis reduction of the correlation problem consists in projecting the displacement residual r k = Hd kb k on B:

r k ∈ B ⊥ ⇔ ∀v ∈ B, v T r k = 0 (Orthogonal space definition) (3.67)
where H and b k are defined with the final mesh refinement. By denoting α α α the coordinates of the displacement in the reduced space (also called the reduced modes vector), the projection of the residual displacement leads to the following reduced system:

R T HRα α α k = R T b k (3.68)
where α α α 0 is either estimated from another initialization method, set to 0 or itself projected from an initial solution u 0 using the projection R T Rα α α 0 = R T u 0 . The idea here is to define a reduced measurement model with minimal effort directly on the correlation's linear system (3.29) using fewer degrees of freedom. We list in the following some possibilities for choosing the reduced modes:

• Scaling of a pre-load displacement solution:

When a linearity hypothesis can be made between different states of loads, one can use a previous load measurement solution denoted v and look for a displacement correction that scales to it. In this case, the scaling matrix is reduced to

R =   v x 0 0 0 v y 0 0 0 v z   with α α α ∈ R 3 . (3.69)
v x , v y and v z ∈ R n are the vector displacement blocks extracted from v (corresponding to the x, y, z directions respectively). 0 is a zero vector of dimension n (n is the number of mesh nodes).

• Affine transformation (Scaling, schearing, rotation and translation):

R =  
x y z 0 0 0 0 0 0 1 0 0 0 0 0 x y z 0 0 0 0 1 0 0 0 0 0 0 0 x y z 0 0 1   with α α α ∈ R 12 .

(3.70)

1 is a vector of dimension n containing only ones.

Analysis of the correlation 86

• Free-Form Deformation (FFD):

This method consists in controling the correlation problem (or in general any geometric problem) with few degrees of freedoms using a B-spline morphing box. It was applied to shape and displacement measurement for Stereo-DIC in [Chapelier et al. 2021]. In the case of DVC, it has an interest if we desire to deform the final finite element mesh with minimal computational effort.

• Subset initialization: We note that the use of the subset DIC/DVC method can be used in order to initialize the finite element based DIC algorithm. In this case, we get the structured grid of subsets and interpolate the displacement field linearly between each nodes of the subset. This seems to be sufficient to initialize accurately a global finite element approach.

• Finally, from a practical point of vue, one can combine these three previous strategies by selecting charateristic modes that can improve the measurement of global transformations. It has to be underligned again that the choice of the modes is not fixed by a very clear rule and usually manual tests are performed in order to find the characteristic modes of the global transformation before proceeding to the final correlation.

Analysis of the correlation

In order to assess the quality of the correlation at convergence, two criterias are considered in this work: a global scalar indicator related to the alignement functional S and the residual map:

• Global residual indicator:

This criteria can be considered as the correlation score.

η = 1 max(I r ) -min(I r ) S(u) V (Ω) . (3.71)
It is normalized with respect to the image dynamic and can be defined either on the whole region of interest or on specific elements.

• Analysis of the residual map:

The residual map is the grey-level difference I r -I d • (id + u) (for a measured displacement field u) displayed on the reference configuration domain. It is fundamental to use this map in addition to its histogram in order to investigate the quality of the correlation. It is also a tool for detecting regions where non-linear phenomena occur (i.e. where the correlation fails).

Chapter 4

Architecture-Driven Digital Image Correlation Technique (ADDICT) for the measurement of sub-cellular kinematic fields in speckle-free cellular materials

After presenting the image correlation framework and its numerical details, we use the equilibrium gap regularization and automated imagebased mechanical modeling in order to measure subcellular kinematics in two-dimensional artificial cellular geometries. We named the proposed methodology ADDICT for "architecture-driven digital image correlation technique". We assess the performance of this method by considering a set of synthetic and real images. The considered two-dimensional geometry mimicks a cellular foam. In order to validate the results of the correlation algorithm, we first compare the measured displacements with the reference finite element solution used to generate the synthetic images. We then compare the ADDICT approach and high resolution FE-DIC. The content of this chapter is partially reproduced from the results of [START_REF] Rouwane | [END_REF], Rouwane et al. 2022].

Problem setting: study of synthetic images

The performance of the proposed speckle-free ADDICT is assessed by analyzing a set of three synthetic test-cases. Namely, given a fine FE mesh fitting the architecture of the cellular material, wisely chosen constitutive properties, and boundary conditions, a displacement field u fem is computed from a standard FE analysis, as detailed in section 4.1.1. Then, synthetic images of the reference and of the deformed configurations are generated, as described in section 4.1.2. The interest of such virtual tests lies in the fact that the measured fields u meas can be compared with the ground truth u fem using appropriate measurement errors, see section 4.1.3. Fig. 4.1 summarizes the process of constructing and analyzing images for our virtual experiment. In addition to performing a virtual elastic test, we will also investigate the ability of our method to estimate local kinematic fields in non-linear regimes (in particular, plasticity and/or geometric non-linearities). We proceed as follows for the discussion of the results: in section 4.2.2, it is shown how challenging it is to estimate sub-cellular kinematic fields with classical subset DIC or macro elements from such images. The latter are then analyzed with the proposed method. Finally, for each of the three test cases, the influence of the regularization cut-off wave-length is analyzed in section 4.3 based on the so-called L-curves of the regularized optimization problems and their relation to the true measurement errors.

Construction of the three virtual tests

For the construction of the reference displacement field u fem , we considered the mechanical problem depicted in Fig. 4.2. The geometry was generated by considering first a Voronoi partition of random points. The cells strut were afterward created manually using segments and arc circles in the computer aided design interface CATIA. This software allowed us to simulate machining before manufacturing the sample (see the experimental study in section 4.4). The left boundary of the sample was fixed (u x = u y = 0) and a homogeneous displacement was prescribed at the right boundary (u y = 0 and u x = u 0 ). The top and bottom boundaries were assumed traction-free (σ.n = 0). The finite element mesh was chosen fine enough to correctly represent the local behavior of the cell struts: approximately six triangular finite elements in a cell strut were considered. show the global force-displacement mechanical response for the three test cases (i), (ii) and (iii), respectively. The red dots correspond to the mechanical states chosen to generate the images I d in the deformed configuration.

Generation of the synthetic images

The virtual DIC testing consists in generating a virtual image of the digital specimen of Fig. 4.2 in the load-free configuration I r , and another one after loading I d from the above computed displacements fields u fem . Let us stress that we did not generate the image in the deformed configuration by advecting a given reference image. Both images were generated from a digital object (given its geometry and color). In order to mimic the generation of grey-scale images from the geometry of the sample, a first high-resolution binary image is generated using a cartesian grid of pixels over the rectangle with vertices (x min , x max ) and (y min , y max ). Afterwards, a pixel grey-level value is assigned proportional to its surface fraction to meet the desired low resolution (about 4 pixels in the strut thickness). The same treatment is performed in order to generate the image of the sample in the reference and deformed configurations. This simple rendering method was sufficient in our two dimensional DIC analysis whereas other more complex physically sound rendering models could have also been considered, (see, for instance, [START_REF] Orteu | [END_REF], Passieux et al. 2015, Sur et al. 2018] in the context of Stereo-DIC).

Let us recall that the images are chosen for the loading states corresponding to the red bullets in Fig. 4.3. For the non-linear regimes (see,in particular,Figs. 4.3b and 4.3c), this ensures that the behaviour has clearly entered a non-linear regime. The corresponding images I r and I d are shown in Fig. 4.4 for each of the three mechanical problems.

Error quantification

As indicated in the overview of the synthetic experimental setup in Fig. 4.1, the computation of the measurement errors was performed by comparison with the reference FE displacement u fem used for generating the synthetic images. Since the reference FE mesh is consistent with the cell geometry, we choose to compute the error between the measured u meas x , u meas y and simulated u fem x , u fem y displacements at the n p Gauss points defined on all triangular elements of the simulation mesh. In Fig. 4.5, a zoomed window is provided to see the FE mesh and corresponding integration points located in the image domain. In order to quantify the measurement errors, we consider the measurement uncertainty denoted U . For instance, for the x-component of the displacement it is defined as follows:

U (u x ) = 1 n p -1 np i=1 (u fem x i -u meas x i ) 2 , (4.1)
where u xi stands for the evaluation at the i th Gauss point. The uncertainty U will be used for characterizing the measurement error for u x and u y with respect to ground truth.

Specifications of the ADDICT method

The proposed regularized method consists in first generating the automatic image-based mechanical model using either the finite element method or the finite cell method. This provides a measurement mesh that allows to build the differential operators needed for regularization. We then solve the regularized problem (3.57).

Finite element approach

The idea here is to reconstruct the finite element mesh from the reference image I r . To do so, we first extract the contours using the marching squares algorithm and then define a closed surface. The obtained surface is given to the finite element mesher based on Delaunay triangulation as it was detailed in chapter 1. Fig. 4.6 shows the mesh that was automatically generated from the image I r . Again, the threshold value chosen for extracting the contours was set to γ = (v min + v max )/2 where v min and v max are the minimal and maximal grey-level values of I r .

Figure 4.6: Automatic mesh generation of the cellular geometry using the marching squares algorithm and Delaunay triangulation provided by Gmsh.

Finite cell method approach

The different steps of the construction of the B-spline immersed model are recalled in Fig. 4.7. We describe them briefly these: first the continuous representation of the image is build using the spline representation as recalled in section 3.4. This allows to obtain the smooth level-set function that is equal to 1 if the evaluated point is in the region of interest and 0 in void areas (see Fig. 4.7a). In a second step, the region of interest is embedded in a structured smooth and higher-order B-spline grid for the discretization of the measured displacement field (see Fig. 4.7b). In order to build the differential operators on the cell geometry (i.e. the interior of the domain defined by the level-set), the integration and the assembly is performed by means of the quad-tree decomposition detailed in section 2.2.2.2.

(a) Level-set characterization of the material's boundary.

(b) Embedding of the region of interest in a smooth and higher-order cartesian B-spline grid.

(c) Quad-tree scheme with a closure linear tessellation for defining the domain of integration corresponding to the cell struts. The three fictitious domain parameters are adjusted following the conclusions of chapter 2: the maximum level of quad-tree decomposition is taken so that the minimal size of an integration element is about the same as the pixel size, and smooth cubic B-spline elements of size approximately equal to the cell strut thickness are employed. For illustration purposes, the considered cellular-like specimen is shown in Fig. 4.8 along with the chosen B-spline mesh that is composed of n x = 87 and n y = 64 elements in the x and y direction, respectively. The corresponding approximate element size is equal to 2.5 pixels. In the end, we make use of the B-spline grid and constructed fictitious domain integration rule to compute the sitfness matrix K, the correlation Hessian matrix H, the gradient at each iteration b k and the Laplacian operator L which intervene in Eq. (3.58). In addition, we interpolate the images by using the smooth spline representation defined in section (3.4) constructed at the first step of the image-based model to define the level-set function, which is interesting from a noise and gradient computation point of view [START_REF] Pan | Bias error reduction of digital image correlation using Gaussian pre-filtering[END_REF], Chan et al. 2017]. The remaining issue to address is that these operators are in general severely ill-conditioned, due to the fact that some basis functions can have their support that do not, or only slightly, intersect the physical domain. As a remedy, we remove the degrees of freedom corresponding the basis function N i such that:

s(i) ≤ ε (4.2)
where s(i) is the basis function indicator defined in Eq. (2.18). In this work, we fix ε = 10 -4 in order to obtain a good compromise between the conditioning of the left-hand side operator and the accuracy of the solution. We note that we can also keep all the basis functions and choose a high stifnness penalization parameter for the stiffness matrix which itself re-controlled with the regularization parameter λ K . In Fig. 4.9, we show the retained control points after applying Eq. ( 4.2) with the considered geometry and mesh. Overall, the strategy of Eq. (3.57) can be seen as an optimized version, using advanced image-based model techniques, of the mechanically regularized DIC scheme [Réthoré et al. 2009, Leclerc et al. 2011, Leclerc et al. 2012]. For the qualitative comparison, we will consider the compression test in its elastic regime (i.e. before that the buckling happens). The compared images I r and I d are illustrated in Fig. 4.10. Here, the imposed displacement is approximately equal to 10 pixels (corresponding to u 0 = -3 mm in the orginal mesh dimensions). In this first study, we will not insist on the physical length related to the regularization parameters but rather vary their values in a heurstic way.

Global B-spline DIC without regularization

In this numerical test, no regularization is considered (λ = 0). When considering cubic B-splines, the idea is to start by a mesh and gradually decrease the elements size using the multi-level refinement scheme detailed in section 3.7.1. Starting from a coarse mesh (see Fig. With the multi-level refinement approach, a global correlation residual η = 0.77% is attained (see section 3.71 for the definition of η). This residual is reached at a level of 4 subdivisions corresponding to the mesh shown in Fig. 4.12c and in which the element resolution is equal to (hx/∆x)(hy/∆y) ≈ 11.11 pixels, where h x , h y , ∆x, ∆y are respectively the element size and the pixel size in x and y directions. Beyond this mesh resolution, the algorithm diverges and becomes unstable. The displacement field solution has the same homogeneity as the finite element solution, however, the differentiation of the solution, gives a strain field solution that presents high strain concentrations with values that drastically surpass the reference finite element strain field (see Fig.

Regularization using the diffusion model

We consider the regularization scheme of Eq. (3.45). This time, we combine the multi-level refinement strategy and the tuning of the regularization parameter λ (see Eq. (3.43)-(3.45)). We start from a coarse mesh and vary λ from 10 5 to 10 1 . When using this regularization, the Gauss-Newton algorithm converges independently of the element size. When considering for example a mesh that has an element size equal to 4.74 pixels, we obtain at convergence the following residuals: η = 3.62% for λ = 10 5 , η = 0.93% for λ = 10 4 , η = 0.49% for λ = 10 3 , η = 0.38% for λ = 10 2 and η = 0.34% for λ = 10 1 . This numerically confirms the existing trade-off between the conservation of grey-levels and the regularity of the solution given by the cost function of Eq. (3.45), which will be discussed more in details in section 4.3. However, in all cases, the registered solution using this type of regularization presents high strain concentrations in the cell struts (see Fig. 4.14). It thus seems that the diffusion regularization enables to get a good solution in terms of residuals but the solution clearly does not live in a mechanically consistent space. The curvature regularization defined in (3.47) was also tested in its approximative version by considering T = L T L where L is the Laplacian operator and also by setting T = B where B is the bi-Laplacian operator assembled with the second order derivatives. No clear differences were observed and the same noisy strain field is obtained. Finally we consider the mechanical regularization scheme of Eq. (3.57) with all the regularization terms: λ K = 0 and λ L = 0. λ M and λ T are fixed in an iterative way (gradual decrease). Obviously, we could consider larger elements for the B-spline mesh and still make the algorithm converge. However, we would obtain less accurate results in this case since the B-spline mesh would be too coarse to properly describe the mechanical behavior at the sub-cellular scale of the sample. The chosen B-spline grid here has an element size equal to 4.74 pixels (the same as the one chosen in section 4.2.1.2). As explained previously, the mechanical regularization term defined by the stiffness operator K describes the elastic behavior of the observed architecture. A choice thus needs to be done for the associated Poisson coefficient ν. Here, we take the same parameter as the ones assumed for the behavior of the true micro-structure. However, we emphasize that since the stiffness matrix is proportional to E, the stiffness is tuned by the regularization parameter. In addition, we conducted different numerical tests with different Poisson ratios (ranging from 0.1 to 0.5) for the regularization term. Indistinguishable results were obtained which is consistent since, once again, the elastic behavior in the regularization is prescribed in a weak sense, only to alleviate noise effects in DIC. More precisely, the obtained solution fields are shown in Fig. 4.15. With λ K = 10 5 , the displacement field solution achieved a correlation residual η = 0.72%. We observe in Fig. 4.15, that the strain field solution is very close to the one used for generating the images I r and I d at the bulk of the sample. Thanks to the high accuracy of the mechanical model, regular and local bendings are observed at the cellular beams (see Fig. 4.16). This shows that, under elastic deformations, the suggested mechanical regularization surpasses the two previous DIC methods and allows to estimate local strain fields at the sub-cellular scale. These results appear very promising for the goal of treating real 3D image data of cellular materials. We report the uncertainty for each method in Tab. 4.2. These results clearly show that the mechanically regularized method is the most precise in terms of measurements of displacement and strain fields compared to the two standard DIC methods. In addition, when using standard methods (with or without Tikhonov regularization), it can be observed that the strain measurement uncertainty is always greater than the reference strain we want to measure, which proves quantitatively what could be observed qualitatively on Figs. 4.13 and 4.14.

U (u x ) (pixels) U (u y ) (pixels) U (ε vm ) Standard multi-level DIC 4.5 × 10 -1 1.9 × 10 -1 1.9 × 10 -1 Diffusion regularization 1.6 × 10 -1 1.1 × 10 -1 5 × 10 -2 Mechanical regularization 2 × 10 -2 3 × 10 -2 4.7 × 10 -3
Table 4.2: Uncertainty of the measurements.

Noise sensitivity analysis

In classical two-dimensional imaging, the noise sensitivity of the DIC algorithm is usually studied by considering a white noise, i.e. more precisely, by perturbing the images I r and I d with a random Gaussian variable of zero mean and with a standard deviation σ ∈ {0, 1, 2, 3} grey-levels. For each value of σ, 10 random variables are generated. We therefore consider the mean uncertainty over these 10 random tests. show, respectively, the evolution of the global displacement uncertainty defined by 1 2 (U (u x ) + U (u y )) and of the Von Mises strain uncertainty U (ε vm ) with respect to the noise standard deviation. . stands for the mean over the 10 noise tests and U is defined in Eq. (4.1). No significant error change is observed with respect to noise (same value along the Y axis). This is consistent with the sub-pixel evaluation based on the smooth spline defined in (3.38). In fact as we choose the spline coefficients as the filtered images, the continuous representation acts as a pre-filter for the images. We should note that this noise sensitivity study is very naive as in reality, the level of noise in tomography is much greater than the value fixed in the example. In addition, other more complex noises can be encountered in three-dimensional image acquisitions (we can quote, for example, ring artifacts and beam hardening [Boas & Fleischmann 2012]). That is why, acquired images are generally pre-processed with noise reduction algorithms.

A second analysis vs Subset based DIC

As mentioned in section 3.6 and illustrated in Fig. 3.2, the usual practice when considering macro elements is to set a subset size or finite/B-spline element size according to the characteristic length of the image texture. Based on the auto-correlation function of the image, we can first estimate the microstructure's characteristic lengths. The auto-correlation function can help achieve this goal and is defined as the cross-correlation (or shifting product) of a signal with itself. It is defined by the following formula:

∀τ ∈ R, A I (τ ) = R I(t + τ )I(t)dt. (4.3)
The Fourier transform of the auto-correlation of a signal is equal to its power spectrum [Olson et al. 2017], so one can compute the normalized auto-correlation efficiently using the following formula:

∀τ ∈ R A I (τ ) = 1 n p F -1 F I - Ī σ I 2 (τ ) ∈ [-1, 1] (4.4)
where τ is a shift in the image domain, F is the Fourier transform, n p is the total number of pixels, Īr and σ Ir are respectively the mean and standard deviation of the image I r . At τ = 0, a perfect match is obtained for the signal and its shift therefore, the auto-correlation is equal to its maximum value (perfect correlation equal to 1). When I is for example a rectangular window of support equal to 1 and centered at 0 (i.e. I = 1 [-1/2,1/2] ), its auto-correlation function is a triangular window function equal to:

A I (τ ) = (1 + τ )1 [-1,0] (τ ) + (1 -τ )1 [0,1] (τ ).
When considering a peridodic rectangular wave by summing non overlapping window functions, the obtained auto-correlation function is a periodic triangular wave. Fig. 4.18 illustrates a rectangular wave. The length of each rectangular window is equal to ∆ 1 and the period distance between each window is equal to ∆ 2 . For this perfect rectangular signal, the auto-correlation allows to determine both ∆ 1 and ∆ 2 wich are the support of the triangular windows and their period. This motivates the use the auto-correlation for images in order to determine the size and distance between particles of images of repeating patterns. In the case of two-dimensional images, and in order to take into account the geometry in all directions, the radially averaged normalized auto-correlation function is rather considered:

∀r ∈ [0, R] ÃI (r) = 1 n θ n θ i=1 A I (x 0 + r cos(θ i ), y 0 + r sin(θ i )) , (4.5)
where for each radius value a mean value over all angles θ i is computed with (θ i ) i∈ 1,n θ a sampling of [0, 2π]. The point (x 0 , y 0 ) is the image center and r is the radius shift. The reader interested by using this tool can find a plugin within the Fiji software [Schindelin et al. 2012]. As a second example, we consider a grid of disks uniformly distibuted on an image domain. The diameter of a disk is equal to ∆ 1 = 12 pixels, the horizontal distance between two disks is equal to ∆ 2 = 22 pixels and the diagonal distance is equal to ∆ 3 = 32 pixels. Let us return now to our practical DIC problem. For this second study, we consider testing the subset DIC softwate VIC-2D. We rather consider another deformation state, always in the elastic regime but for the tensile configuration (i) in Fig. 4.4b. By performing the analysis of the evolution of the radially averaged normalized auto-correlation, we can estimate an averaged texture size in the image and the periods existing thanks to the auto-correlation peaks. The secondary peak at around 30 pixels characterizes the mean cell size or the distance between two cell pores. Based on the usual practice in subset DIC [Schreier et al. 2009, Bornert et al. 2009, Jones & Iadicola 2018], it is stated that the subset should contain a minimum of three DIC texture features, which leads, in our case, to choose very large subset sizes incapable of reconstructing the local kinematic associated to strut bending (see also discussion related to Fig. 3.2).

The subset-method was applied with affine subset shape functions. In the case of using the image of Fig. 4.4a, the subset DIC tool VIC-2D suggests an automatic subset size so that the measurement unicertainty is minimized in the software's criterias. A subset size of 63 pixels is suggested in this case (approximately 3 pores per subset as shown by the orange square in Fig. 4.22), which is consistent with the usual practice. The step size was set to 1 pixel. The measurement points are marked by the red dots in Fig. 4.22. It should be noted that such a large subset size only allows measurement in an area relatively far from the edges. As we are interested by the measurement within the cell struts only, we show the post-processed results in the cell regions using a a posteriori binary segmentation. In Fig. 4.23, it can be seen that the displacement field estimated with the subset method is consistent with the reference field, at least at the macroscopic scale. But when analyzing the field measured by the subset approach in more detail, by looking in particular at the strain field in Fig. 4.25, we notice that the strain provided by the subset method is completely inconsistent and very far from the reference strain field. More precisely, the obtained strain fields are homogeneous at the scale of the cell-struts and the local bending observed in Fig. 4.25a is not identified. In addition, when analyzing the relative error with respect to the reference field (see Fig. 4.24), most measured regions using the subset method exceed an error of 20% which is not the case for the ADDICT method (error less than 10%). This shows that large subsets only allow to identify macroscopic (or homogenized) displacements and strain fields.

This problem is due to the difficult compromise in choosing the subset size. Indeed, this parameter alone is used to set both the regularization length and the measurement resolution. This motivates the use of a richer kinematic (small resolution) associated to an alternative regularization technique to better capture the sub-cellular displacement field gradients. are presented in Figs. 4.23c and 4.25c. It can be observed that the displacement field is much better resolved. It shows typical bending gradients which are quite similar to the reference fields. This is a clear illustration of the interest of the FE approach in DIC in its ability to use a mechanical model to improve DIC and to break the aforementioned trade-off. Until now, when using the mechanically regularized approach ADDICT, we have just tuned the mechanical regularization parameters λ K and λ L heuristically without using the normalized optimization problem. In the following section, we will study the two main parameters of our method: (a) the choice of the regularization lengths l L and l K (see Eq. (3.59)), and (b) the relevance of the model (here linear elastic) used for the regularization operator with respect to the nature of the non-linearity of the measured behaviours. We solve now the optimization problem in its normalized form given by Eq. (3.58).

Numerical investigation of the influence of the model and parameters used for the regularization

In this section, the influence of the regularization lengths l L and l K for different linear and nonlinear mechanical regimes is investigated using L-curves. The L-curve study of regularized leastsquares problems helps finding the optimal regularization parameter as the one related to the highest curvature point in a log-log plot of the regularization term versus the data fidelity term [Hansen 2000].

In other words, we usually end up with a L-shaped curve when performing this log-log plot and therefore the optimal regularization parameter corresponds to the one of the corner point of the "L". For our mechanically regularized scheme (see Eq. (3.57)), we thus consider on the horizontal axis the dimensionless data-fidelity term defined by S(u)/(max(I r )min(I r )), and on the vertical axis the variation of the mechanical equilibrium, i.e. ||D K Ku|| 2 2 . In order to investigate the filtering properties of the equilibrium gap based regularization, the plots are performed for different values of the characteristic lengths: l L and l K are respectively varied over 0, 40 pixels and 0, 200 pixels. The L-curve corresponding to the curvature regularization of Eq. (3.47) (i.e. ||Lu|| 2 2 with L is the Laplacian operator) is also given for comparison purpose regarding the employed regularization model. In a next step, to account for the relevance of the regularization parameters selected with the Lcurve approach, a measurement error study (w.r.t. ground truth) is carried out. Eventually, several deformed configurations of the material sample are provided with different values of regularization parameters to appreciate visually their influence on the results.

Linear elastic case

First, let us consider the L-curve when regularizing DIC with our approach of Eq. (3.57) in case (i), i.e. where the synthetic images were generated with a linear elastic model (corresponding to Figs. 4.3a and 4.4b). The obtained plot is shown in Fig. 4.26. The left (a) and right (b) sides of this figure exactly correspond to the same plot, only the colour of the markers changes. On the left (a), the colour depends on the value of the edge regularization length l L , and on the right (b) on the bulk elastic regularization length l K . The first thing that stands out is that the parameter l L has very little influence on the L-curve. It only has an effect when the bulk elastic regularization parameter l K is very large (see bottom zone in the figure), which corresponds to very strong regularization. In such a situation, it can be seen as an integrated type DIC method [Neggers et al. 2017] which gives good results provided that (a) the imposed mechanical behaviour in the bulk is the right one (which is the case on this test) and (b) the edge displacements are relevant. This is the reason why edge regularization has an effect in this zone. Fig. 4.26a shows that l L should be considered very small (1 to 5 pixels) in order to get an accurate measurement. In the remainder of this numerical investigation, we will choose l L = 1 so that the influence of the Laplacian-based regularization is restricted to its minimum.

Concerning the influence of the bulk regularization given by l K , while increasing this regularization weight, the equilibrium term keeps decreasing without a significant increase of the grey-level residual (the curve somehow plunges down). Therefore the usual "L" shape is not truly met here: in other words, the L-curve does not present a clear local convexity. The horizontal branch of the "L" shape disappears which means that putting more weight on the regularization term does not alter the greylevel residual. This was actually expected since the regularization kernel is perfect (here elastic) with respect to the prescribed mechanical response of the cellular micro-structure. This will not be the case when the regularization kernel does not exactly correspond to the mechanical response to be observed, see below.

Non-linear cases

The proposed ADDICT with elastic regularization is now applied to the images of test cases (ii) and (iii), i.e. with elasto-plastic constitutive relation without and with geometric non-linearities, as shown in respectively. In Fig. 4.27a, the corresponding L-curves are presented for the three input models (elastic, elasto-plastic and elasto-plastic with possible geometric non-linearities). Only the influence of l K is considered, l L being fixed to its optimal value following previous discussion (l L = 1).

Numerical investigation of the influence of the model and parameters used for the regularization 107

We can now observe three main regions in the L-curve (denoted R1, R2 and R3 in Fig. 4.27a). On the region R1 (i.e, l K < 25), the weight is put more on the grey-level conservation and the standard deviation is higher, the obtained solution is not accurate as will be shown in Fig. 4.28. Conversely, on the region R3 (i.e, l K > 30), the weight is put more on (elastic) regularity. In this case, the grey-level residual increases as the elastic regularity is no longer valid for describing the actual mechanics (here plasticity without or with geometric non-linearities). The choice of l K must be a compromise between regularity and grey-level conservation. The optimal value for the regularization length is at the point of maximum curvature [Hansen 2000], i.e. between 25 and 30 pixels, which defined region R2.

Through this study, it can also be emphasized that the L-curve is proving to be an excellent indicator of the relevance of a model in the context of validation [Neggers et al. 2017]. If the L-curve tends to plunge down as the regularization length increases, then the model is probably compatible with the observed mechanical field. are missing data for the Laplacian-based regularization over M3 because the convergence of the minimization algorithm was not achieved on these points with our residual tolerance and maximum number of iterations allowed.

Comparison with the curvature regularization

As mentioned above, the choice of the model used for regularization is one of the two important parameters of the approach. Here, the less physically sound curvature model of Eq. (3.47) was used to regularize the same set of images. Note that operator L is built by integrating only on the physical cell struts (i.e. avoiding the holes), which differs from the current practice in other fields where such regularization operators are used in both strut and void parts [Dall'Ara et al. 2017,Patera et al. 2018]. The corresponding L-curves are given in Fig. 4.27b. Looking closely at the L-curves of Fig. 4.27a with the different regularization operators, we can see that the L-curve is clearly more sensitive to the increase of the regularization length when using Laplacian-type regularization as compared to the 4.3. Numerical investigation of the influence of the model and parameters used for the regularization 108 elastic one.

Link between L-curve and error

In this section, the L-curves are compared to the true errors in order to numerically validate the optimality of the regularization length associated to the maximum curvature. In Fig. 4.28, the evolution of the measurement error is plotted as a function of the regularization lengths. We recall that, to compute the measurement error defined by Eq. (4.1), the displacement fields are computed at the Gauss-integration points that belong to both the reference finite element geometry and the constructed geometry using the level-set function. First, this figure provides numerical evidence that the optimal value of the regularization calculated from the maximum curvature point also corresponds to the minimum error. Second, this figure also provides numerical evidence that a weak elastic regularization, even when it is not representative of the actual mechanics of the observed specimen, is better than all the other less physical regularization techniques considered in this study, either in a strong way based on polynomials (subset) or in a weak way based on the Laplacian of the solution.

(U (u x ) 2 + U (u y )) 2 /2.
Overall, the interpretation that can be made of these results is that the term associated with the grey-level residuals (S(u) in Eq. (3.57)) captures the low frequency part of the solution, here associated with characteristic lengths higher than the cell length (≈ 30 pixels), i.e. the meso scale. In other words, it helps computing the part of the displacement field that aligns the mesh to the edges of the struts. The local part of the displacements, i.e. inside the struts or at the micro-scale, which do not modify the grey-level conservation term, are driven by the regularization. It therefore seems consistent that the optimal regularization length is close to the characteristic cell size.

Numerical investigation of the influence of the model and parameters used for the regularization 109

Deformed configurations with different values of regularization parameters

In order to visually appreciate the above interpretation, we will show several deformed configurations with different regularization weights. First, considering the elasto-plastic case (ii) (Figs. 4.3b-4.4c), we superpose the reference (red) and measured (green) cloud points for a very low regularization (see Fig. 4.29a) and for an optimal regularization (see Fig. 4.29b). Following previous discussion, the low regularization allows to satisfy more data fidelity (region R1) and the optimal regularization corresponds to the inflexion point obtained from the results of Fig. 4.27a (region R2)). When putting more weight on data-fidelity, Fig. 4.29a shows that non-physical displacements are observed within the cell-struts as the green points move differently than the reference points. Conversely, when considering the optimal regularization weight, the movement inside the cell struts is closer to their reference value, see Fig. 4.29b where the red and green point clouds are superimposed. This same effect can be observed on the deformed finite element mesh with the measured displacement when using the finite element version of ADDICT. When considering for example the Levenberg-Marquardt which consists simply in making the hessian H non-singular, the obtained displacement field is too irregular inside the cell struts. Fig. 4.30 shows that the nodes move in non-physical way inside the strut. This proves that any regularization at the sub-cell scale allows to align the struts but the mechanical regularization is the best among all as it only allows smooth mechanical deformations inside the cell-struts.

Secondly, in the case of the geometrically non-linear elasto-plastic model (iii) ), when putting a very large weight on the mechanical term (region R3), the correlation fails to correctly represent the geometric non-linearities (see Fig. 4.31a). In fact, we observe that the regularization model forces the cell struts to bend in an elastic way whereas they should exhibit a post-buckling behavior. When choosing the optimal weight l K (region R2), the buckling is correctly measured using case (geometrically non-linear with elasto-plasticity). Even if the value of the local strain is not totally 4.4. Application to a two-dimensional experiment 111 correct, it is much better than with the other regularization technique considered in this study, and it allows at least the determination of the location of high gradient areas even if the magnitude is not correct.

Application to a two-dimensional experiment

We now propose to demonstrate the potential of the mechanically regularized approach in an experimental situation where inelastic strains take place. To this end, we have chosen to perform a tensile test on a macroscopic two-dimensional cellular like specimen and to compare the two-dimensional kinematic measurements provided by ADDICT using low-definition speckle-free images of the main side with those obtained by a FE-DIC measurement based on high definition images of the opposite speckled side, considered as the reference (see Fig. 4.33). A classic FE-DIC approach is here preferred for the reference to obtain a dense continuous displacement everywhere in the struts. Figure 4.33: Specimen and preparation for DIC -The 50 mm large sample is milled from a 4 mm 2024-T3 aluminum sheet, then painted white between the regions where it will be fixed in the jaws.

One side is simply left as it is, while on the opposite side, an artificial texture is deposited by means of an airbrush.

We first choose a suitable geometry, material and production method to build our model material. The geometry adopted is identical to the one used in the previous section (see Fig. 4.2). The total width of the specimen is 50 mm, and the minimal struts thickness is approximately equal to 0.5 mm. The sample was machined from a 4 mm thick 2024-T3 aluminum sheet from the CAD file using a 5 axis CNC milling machine1 . This process was preferred to waterjet and laser cutting in order to obtain the desired geometry while minimizing the heat affected zone and avoiding the need to deburr the part. The minimum radii of the fillets were therefore limited in the CAD by the radius of the cutting tool.

Once machined, the sample was first cleaned with water and then dried. Afterward, it was painted with white matt acrylic paint. Fig. 4.34 shows the painted sample before the speckle deposit. Three approaches were considered for depositing the speckle. A first attempt consisted in using digital printing. This was performed by a printing company that prints images over plane surfaces. As a 300 dpi resolution was imposed by the printer, unsatisfactory results were obtained (see Fig. 4.35 for the considered Perlin noise based speckle). As digital printing is based on uniform ink stippling, a 300 dpi resolution only allows 7 ink dots per 0.5 mm (smallest cell strut thickness). The second attempt consisted in using manual spray-painting using a matt black spray (similar to the one displayed in Fig. 4.34). A finer speckle was obtained by spraying the paint at the farest distance (1.5 meters) but an unsufficient number of speckle particles was obtained in the strut thickness (in addition to a large dispersion of the spots distribution). The third promising option was to consider airbrush spraying. Airbrush spraying consists in passing fast compressed air through the paint. This allowed to obtain a very fine spray. Fig. 4.36 shows the components of the considered airbrush. An Instron 8561 100 kN electromechanical tensile machine equipped with a 10 kN load cell was used for this test. This machine was equipped with hydraulic jaws, which avoids accidental twisting of the sample during clamping. Particular care was taken to align the jaws beforehand. The test was carried out under displacement control at a constant displacement rate of 0.12 mm/min2 .

Using the third speckle deposition technique, a first tensile test was performed. Unfortunately, the speckle was unsatisfactory from another perspective. During the loading, the paint started to peel off as displayed in Fig. 4.38. This was explained by the quality of the paint but also by the fact that that the sample was not properly prepared. Vivien Murat, an engineer specialized in the treatment of surfaces at IRT Saint-Exupéry have highlighted the fact that it is preferable to wash the aluminum samples with acetone and dry them. The samples were kindly prepared by him3 . First a sandblasting treatment was peformed. This technique, consisting in projecting air mixed with sand on a surface, allowed us to clean the samples and to make them sufficiently rough to improve paint adhesion.

Finally, thin matt black spots were sprayed with a professional air-brush on the side where FE-DIC measurements were planned (see Fig. 4.33b). Fig. 4.39b shows the distribution of the textures obtained on the cell sample. The average diameter of the spots is estimated to be around 0.1 mm.

The experiment was monitored by multiple cameras triggered using an external TTL square signal. The frame rate was set at 0.2 fps. Fig. 4.37 shows the basic optical setup chosen for the present analysis. It consists of 2 systems that were very carefully positioned on either side of the sample and oriented (using laser devices) so that the optical axes were perpendicular to the imaged faces. In order to remove depth effects, a telecentric lens (Opto Engineering TC ZR 072-C) was used to image the speckle-free side of the sample. The dimensions of the obtained field of view are 70.4mm ×52.8 mm. The telecentric lens is equipped with a 5Mp CCD camera (Camera 1: Allied Vision Pike). On the opposite side, a 29Mp CCD camera (Camera 2: Allied Vision Prosilica GT6600) equipped with a macro lens (ZEISS PLANAR T 2.0/100 ZF MACRO) were selected to retrieve high resolution images of the textured surface. In this case, the intention was to correctly resolve the small texture created on the surface. The working distance of the macro lens was set to encompass almost the same region of interest (see Fig. 4.41). The resulting image has a resolution of about 78 pixels/mm. The zoom presented in Fig. 4.39b allows a better assessment of the type of texture which are later used for the FE-DIC. Note that the spots are on average more than 7 pixels, which is a little larger than the value recommended for DIC [Jones & Iadicola 2018]. The lighting during such an experiment is a problem in itself. It was indeed tricky to light correctly one side without dazzling the cameras placed on the opposite side. Fig. 4.37 illustrates how this problem was solved: 2 LED spotlights were used on each side. This same figure reveals an additional stereo DIC bench in the background. The latter allowed us to verify that there was no significant out-of-plane movement during sample clamping or during the test (the maximum out-of-plane displacement measured is at most a few tenths of a millimeter in the gauge region). This feature will consequently not be discussed further since it serves only the purpose of verification of the testing system.

The macroscopic force ( F ) -displacement ( Ū ) curve recorded from the testing machine during the experiment is plotted in Fig. 4.40. The dots indicate when the images were captured. For the DIC analysis which follows, we set the reference image I i r (i = 1 untextured face, i = 2 textured face) as the first images captured after the mechanical jaws were clamped (point ( Ū , F ) = (0, 0) of the curve in Fig. 4.40). Up to about 3 kN, the sample exhibit an elastic macroscopic response. Beyond that, the sample undergoes an irreversible strain, highlighted by the unloads. From now on, we will limit ourselves to present the DIC measurements only for the deformed state indicated by the red dot on Fig. 4.40 (point ( Ū , F )=(1.05 mm, 4.73 kN ). The total macroscopic strain is then estimated at 1.5%, while the corresponding residual macroscopic strain is about 0.8%. The corresponding images are then noted I i d . We now propose to measure the displacement fields by image correlation between the reference state (I r ) and the deformed state (I d ) images. The recorded images on the speckle-free side (I 1 r and I 1 d ) are processed by ADDICT. As we want to test our method in conditions similar to those described above (i.e. with only a few pixels in the strut thickness), the images are downsampled before being processed. Here, we proceed to three successive data binning leading to images of 256 × 306 pixels definition (see Fig. 4.39a). The resolution of the resulting images is then about 4.5 pixels/mm. We then automatically define the implicit geometry of the region of interest by building the automatic FCM image-based model (see Fig. 4.41a). The binary threshold value for the level-set segmentation is here simply set to (min(I 1 r ) + max(I 1 r ))/2. Since plastic strains are expected, the regularization parameter l K is taken approximately equal to the cell length (same reasoning as in Fig. 4.28 of section 4.3). When taking into account the resolution of the experimental images (which is slightly higher than that of the synthetic images of the previous section), the corresponding cut-off wave-length corresponds to l K = 50 pixels. This is confirmed by a new study based on the L-curve. Fig. 4.42 shows that the optimal regularization length lies indeed in the interval 25, 75 pixels. For their part, the high-resolution images (I 2 r and I 2 d ) of the textured side of the specimen are analyzed using the open-source FE-DIC library Pyxel [Passieux 2018]. The unstructured T3 measurement mesh is generated from the very same CAD data used for machining. The average element size is set to 0.2 mm to ensure theoretically that any element encompasses at least one spot. In this 2D configuration, the transformation between the mesh reference frame and the image reference frame (designated projector in this library) is described here with 4 parameters: one rotation around the optical axis, two in plane translations and one scaling. Those parameters are automatically identified by imposing that the projection of nodes on the edges must be aligned with the corresponding edges detected in the images (see Fig. 4.41b). Alternative techniques for positioning a mesh relative to an image exist using synthetic images [dell 'Isola et al. 2019, Baconnais et al. 2020] or levelsets [Passieux 2018]. In practice, we can check that only a few elements do not benefit from grey-scale gradients (see Fig. A quantitative analysis based on the hypothesis of 2D kinematics is now proposed. In the present situation, as it was precised in section 4.1.3, we can indeed directly project the displacement fields provided by ADDICT on the integration points of the FE-DIC technique (see Fig. 4.45). Fig. 4.46 presents the relative difference between the ADDICT and the FE-DIC measurements |u 1u 2 | Ū , where Ū stands for the imposed grips displacement. In no case do the observed differences exceed 3% of Ū . The local fluctuations for both components are explained by the uncertainty of the FE-DIC measurement. To complete these comparisons, we propose to look at the strains inside the struts (see Fig. 4.47). Not surprisingly, the regularized measurement leads to less noisy strains and less sharp gradients. Nevertheless, ADDICT allows us to correctly locate the most severely strained regions. In general, we note that the largest deviations are observed on the left and right edges of the ROI. This was expected and is due to the non-physical regularization required on these edges to force ADDICT to converge. The information provided in the immediate vicinity of these regions should therefore be taken with caution.

In addition to the relevance of the results provided, it should be noted that the use of ADDICT does not require any wizardly parameterization. Indeed, it should be remembered that the behaviour chosen for the regularization is elastic, and no optimization of the gray level threshold to adjust the position of the level-set has been performed (i.e. the description of the geometry has not be optimized -see Fig. Since the ADDICT measured a relevant displacement field, it becomes possible to validate a simulation by using only the low resolution speckle-free images. Consider, for example, the FE model introduced in section 4.1.1. The constitutive parameters adopted to describe the elasto-plastic behaviour of the struts are those presented in Table 4.1. Simple boundary conditions such as those presented in Fig. 4.2 are adopted. The imposed displacement u 0 is fixed at the value of the measured grips displacement u 0 = -Ū . The longitudinal and transverse displacement fields computed with Abaqus are respectively compared to the measurements in Fig. 4.43 and Fig. 4.44. The observed differences between the simulated and measured fields are much greater than the difference between the measurement fields. The simulated resultant F = 0.96 kN is also very different from the force measured at this stage (Fig. 4.40). This means that there is clearly room for an improvement of the simulation (ie. discretization, model, constitutive parameters). Considering that the mesh is sufficiently fine, and that the selected model is relevant, we could consider identifying the constitutive parameters. A classical Finite Element Model Updating (FEMU) approach, such as that proposed by [Xu 2018], but again based on measurements carried out with speckle-free images, could be adopted. Other identification strategies, entirely in line with the approach initiated here with ADDICT, could also be adopted [Réthoré 2010, Réthoré et al. 2013]. Although fascinating, this topic is beyond the scope of this work and would require a separate study.

Chapter 5

Volumetric kinematic measurements in cellular materials using image-based modeling

In this chapter, we demonstrate the applicability and potential of the mechanically regularized image correlation approach for the measurement of three-dimensional displacement fields in an open-cell foam material subjected to an in situ compression test.

The previous chapter proposed a proof of concept based on real and synthetic tests limited to twodimensional cases. Specific specimens whose geometry mimicked the architecture of cellular materials were used. The adopted methodology allowed us to understand how the correlation behaved in the presence of such architectures without speckles exploitable at the strut scale, both in its classical version and in its mechanically regularized version. A study based on the L-curve also allowed us to determine an optimal regularization length for the regularized DIC approach and that corresponds to a characteristic length of the cell architecture. It was indeed shown that the cut-off length l K could be set to the cell size to achieve relevant kinematical measurements at the architecture scale. Instead of setting the element size as equal to the cell size, it is the regularization length which is set to this value. Of course, this conclusion from a two-dimensional study could not be extended in three-dimensions without a clear and rigorous mathematical proof or at least a practical demonstration. However, we have directly exploited this "convenient" result to perform our first regularized DVC measurements in real cellular materials (foam) from micro tomographic images captured during in situ tests [Buffiere et al. 2010].

General context of the compression test under consideration

The experimental data set was kindly provided by Pascal Doumalin 1 of the Photomecanics and Experimental Mechanics (PEM) group of the P' Institut in Poitiers (France). The corresponding experiments were performed in the framework of a project (ANR SOFITT 2 "Saturated Open-pore Foams for Innovative Tribology in Turbomachinery"). The idea is to develop a new biomimetic inspired lubrication system called Ex-poro-hydrodynamic (XPHD). The standard lubification based on fluid 5.1. General context of the compression test under consideration 123 film (or hydrodynamic or thick film) consists in injecting a fluid film between solids in order to decrease the friction between two moving solid parts. As [Kunik et al. 2020] stated,"the essential idea of XPHD lubrication is to replace the antifriction material and the thin fluid film, traditionally used in classical sliding motion, with a porous layer imbibed with a fluid that provides a higher load carrying capacity. This type of lubrication represents a completely new technological solution that breaks with current practices and can replace petroleum lubricants, creating self-lubricating and therefore more ecologically sustainable and less expensive tribological systems". In this context, the researchers of the Mechanical Engineering and Complex Systems3 (PEM and Tribolub teams) work together in order to investigate the complex interaction between porous materials such as open-cell foams and fluids. The idea is to create stiffness and damping by dynamically modifying the geometry of an imbibed porous layer. "The knowledge acquired in this project can not only be used to find innovative technical solutions in turbomachinery but, in the longer term, can also be applied to studies dedicated to biological tissues or to human joints where the behavior of soaked porous materials (human cartilage) is also a key point" 2 . The selected porous material is a polyurethane foam. The in situ tests carried out in the X-ray microtomograph and presented hereafter aim at identifying the spatial distribution of the mechanical properties in relation with the local architectural properties (e.g. the local porosity). A first step towards the determination of local stiffness in porous media, consisted in conducting an identification of local mechanical properties from subset DVC measurements [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]. A finite element model updating (FEMU) technique [Petureau et al. 2019] was applied to this cellular material in order to identify a distribution of elastic constitutive parameters in a given ROI. The idea of the FEMU method is to search for the optimal set of parameters that, when injected into the simulation, the calculated displacement field is similar to the measured displacement field (similar in the least squares sense). For the FE model, the material is assumed to be element-wise homogeneous. Within each element, the behavior is assumed to be isotropic and elastic (two parameters to be identified per element: E, the Young's modulus and ν, the Poisson ratio). The nodes of the structured mesh correspond to the DVC measurement points. The DVC displacements measured on the edges are used as boundary conditions. In [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]], the DVC measurements were achieved using subsets whose size is larger than the cell size, for the reasons that were highlighted in section 4.2.2.

If this approach provides interesting data on the initial elastic behavior, it comes up against the problem of the reliability of the measurement as soon as the kinematics at the scale of the architecture becomes more complex (e.g. strut buckling). This question is the subject of work by the PEM team in Poitiers. On our side, we think that the ADDICT approach would allow to revisit the question of measurement and identification in the presence of such architecture and, why not, to allow the identification of more complex models.

In this last chapter, we present the very first results obtained from a DVC measurement regularized by an image-based model. The idea here is to transpose the same methodology that was presented in chapter 4 but this time on a data set obtained from a real in-situ mechanical test.

Experimental setup 124

Experimental setup

The experiment consisted in a compression test performed on an open-cell polyurethane foam. The foam had a porosity V (void)/V (total) equal to 0.96 [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]]. The mean diameter of the cells was 500 µm with a standard deviation of ±200 µm [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]]. This clearly shows that the material is highly heterogeneous. The tested sample had a cylindric shape with a diameter of 9 mm and a height of 10 mm. A compression machine was specifically designed and used in order to perform the in situ test (see Fig. 5.1). For this setup, an electric cylinder was used in order to apply the mechanical load. The capacity of the load cell used was 50 N . A reference scan was carried out, specimen mounted, before starting the test. Then the loading was interrupted several times in order to allow the sample to be scanned. Each scan lasted approximately 40 to 50 min. The four stages of loading considered corresponded to an overall strain of the order of 3.4%, 5.2%, 12.3% and 28.5%.

Acquisition and threshold value setting

In total, the starting point of our DVC analysis is the five reconstructed volumes of the foam (i.e. the image of the reference configuration I r and the four images of the deformed configurations). The voxel size of the reconstructed images is 6.67 µm. The reference image I r definition is 1523 × 1514 × 1654 voxels cube. Originally, the images are encoded in 16 bits and have a size of 7.6 GB. In order to reduce the memory footprint, the images that we will use later are converted to 8-bit grey-scale and weight only 3.8 GB. Fig. 5.2 shows the rendered volumes of these configurations. The visualization of the different loading states is performed using Paraview [Ahrens et al. 2005[START_REF] Ayachit | The paraview guide: a parallel visualization application[END_REF]]. It takes, in this case, as input, the different VTU (Unstructured VTK) files of the loading configurations. The segmentation, that was performed for rendering and visualization purposes, was done using the marching cubes algorithm (see chapter 1 for more details). The first step consists in computing the 5.3. Acquisition and threshold value setting 125 spline coefficients which are the coefficients c = G σ * I where I is the reconstructed volume and σ = 1. We then evaluate the continuous representation Ĩ at an over-sampled grid of points (in the image domain). The threshold value for both rendering and mesh generation was optimized using a bisection algorithm so that the porosity of the segmented geometry reaches the actual porosity. The evolution of the porosity of the rendered volume with respect to the threshold value is shown in Fig. 5.3a. In our case, when setting the threshold value corresponding to a porosity of 0.96 (that of the foam as provided by the material supplier), the obtained cell architecture is fragmented and clearly does not correspond to the real one which is a connected architecture (see Fig. 5.3b). In practice, without objective data on the real porosity of the sample, we choose the largest threshold value that allows to obtain a visually connected region. When setting γ = 56, we obtain the closed surface shown in Fig. 5.3c. After performing the volumetric meshing of this segmented surface, a specimen porosity of 0.93 is obtained which corresponds to a relative error equal to 3.1%. As the images have a size of 3.8 × 10 9 voxels in total, the treatment of the full volume induces multiple computing performance challenges such as memory management, the speed of input/output requests and post-processing. Multiple small computing details may sometimes lead to errors that sometimes may be difficult to detect. As an example, when allocating the memory of the full volume, the computational DVC code crashed and it was difficult to find the source of the problem. It was found afterward that as the full image memory size allocated had a value larger than the value encoded by an integer variable, an arbitrary integer value was assigned to the size. Therefore it was necessary to declare the size of the image as an unsigned integer in order to get sizes larger than 2.147 billions. In addition, as the assembly of the correlation right hand side b k (defined by Eq. (3.31)) is performed at each iteration of the optimization scheme, the integration which is a voxel summation becomes very costly. That is why the assembly was parallelized over the finite elements. In this work, and in order to reduce the memory cost, the measure was performed on the image stacks binned with a factor of 3. The treated images had therefore a size equal to 507 × 504 × 551 corresponding to a spatial resolution of 20 µm. In [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]], the authors have also considered this same resolution when using subset DVC.

The following numerical demonstrations were performed on a 64-bit station equiped with Intel Xeon(R) CPU E5-2637 v2 processor (3.5 GHz frequency, 16 CPUs), 125.8 GB of RAM. The parallelizable rootines were computed over the 16 CPUs using OpenMp directives. Concerning the linear system factorization, as the numerical solvers were used as a black box and no further investigation concerning this aspect was performed, the CHOLMOD solver exploited only 8 CPUs on this machine.

First DVC analyzes based on "macro" elements

The first step of the DVC analysis consists in setting a measurement mesh in the region of interest. It was indicated in [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]] that the authors have used subsets of a size 31 × 31 × 31 voxels with a uniform step size of 40 voxels. The subset has therefore an approximate size equal to the cell size. On our side, the approximate cell size was determined using the normalized radially averaged auto-correlation function (see again section 4.2.2 for the definition). The first maximum in Fig. 5.4a corresponding approximately to 28 voxels. This value provides an estimate of the mean cell size in very good agreement with the value of 500 µm indicated in [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]]. For this analysis, 8 × 8 × 12 subsets are evenly distributed inside a parallelepipedic core region of interest. This type of approach does not allow in any case to treat the whole volume (here cylindrical). This is one of the drawbacks of the subset based method which is not geometrically flexible. In this work, we have used the subset approach only to validate the global DVC results. The subset measure was performed using the TomoWarp2 software [Tudisco et al. 2017]. The correlation was successfully performed using a correlation window of 30 × 30 × 30 voxels and a node spacing of 30 voxels. All the subsets were correctly correlated for the two first loading states and the same approximative displacement magnitude was obtained with the finite element based DVC based on macro-elements. However, we have not been able to successfully run the TomoWarp2 code for the third state. This can be explained by the fact that we have not correctly set the initial guess (based on voxel search). The three first states were however correctly correlated in [START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]] (see page 134). The subset approach was not further investigated in the context of this work.

Let us now focus on the finite element DVC approach. To perform a measurement without any additional regularization, coarse meshes are required. In the following, and to best cover the volume of the sample, we will therefore rely on an unstructured mesh based on tetrahedral macro elements (macro in the sense of the underlying architecture, i.e. with a sufficient grey-level gradient distribution). Fig. 5.5 shows the three mesh discretizations that were considered. The corresponding mean element size (given by the third root of the mean element volume) is approximately equal to 33, 23 and 12 voxels. As only the first discretization allowed to perform a correlation without regularization. The two other discretizations were regularized using the diffusion model defined in Eq. (3.45). For each mesh, four levels of image binning were performed in order to gradually filter the high frequency components of the measured displacement. For the two first loading states, the initial displacement u 0 was set to zero and the Gauss-Newton algorithm converged to a solution. However for the third and fourth loading states, initializing the displacement field by zero is not sufficient to correctly converge the algorithm. After multiple attempts, consisting for example in setting initial guesses based on a scaling and a translation, the only successful strategy consisted in using a structured mesh with few elements in the height direction and running a DVC to estimate a first coarse displacement. This structured mesh, displayed in Fig. 5.6 covered the whole image domain (including a part of the bottom and top platens). The image signal given by the bottom and top platens helped the correlation algorithm as they are subject to translation during loading. To summarize, the first step of the correlation is performed with the mesh presented in Fig. 5.6. The obtained displacement is then evaluated at the nodes of the tetrahedral finite element meshes of Fig. 5.5 in order to get the initial guess u 0 . This displacement exchange corresponds to the first case explained in section 3.7.1.1. Finally, a tiny correlation (that converged in three iterations) was performed using the reduced basis method based on one scaling (corresponding to Eq. (3.69)). After the determination of the initial guess displacement u 0 , the multi-level image resolution strategy allowed us to obtain the results displayed in Fig. 5.7. The first mesh resolution of 33 voxels fails for the last loading state. Tab. 5.1 also shows the global residual indicators obtained for the different correlations. It clearly shows that the finer the discretization is chosen, the better the correlation is which is an obvious expectation in image correlation, provided that a well-suited weak regularisation is employed. For the fourth loading state, the correlation algorithm converges probably to a local minimum (because the correlation score is very high). A very large strain can be reported (28.5%). This clearly shows the limit of the correlation approach when choosing relatively large discretizations. For the first and second loading state, as we are still in the linear regime, the correlation based on large elements allows to get a "good" correlation score which is not the case for the third and fourth states. The measured displacement fields are shown in Fig. 5.8. This first study based on macro elements confirms the limits of a classical FE DVC approach similarly to the conclusions of chapter 4. This motivates the idea of finding a way for adapting the measurement discretization so that it correctly represents the underlying kinematics. When increasing the measurement resolution (i.e. decreasing the element size with weak regularization), we easily come to the conclusion to choose an element size corresponding to the voxel size which is equivalent to the optical flow or the voxel scale DVC suggested in [Leclerc et al. 2011]. Choosing a finite element mesh with a voxel size resolution will induce a huge amout of degrees of freedom. In the case of the chosen region of interest which volume is equal to 80 millions, this will induce approximatively 80 millions elements. This is clearly an inadequat solution as it is a very costly and not necessarly a very accuracte solution as shown in the numerical example of chapter 2. This type of discretization requires however the development of robust high performance computing strategies such as those suggested in [Liu et al. 2019]. This partially motivates the use of the ADDICT approach as it will allow to obtain in the same time a geometric characterization of the cellular architecture and a geometrically justified DVC displacement measurement.

Analysis using the ADDICT approach

After considering a usual FE DVC measurement based on "macro" elements, we will now treat the images using the ADDICT approach. As our implementation of the FCM is quite limited from the numerical solving point of view and cannot be used for treating large scale problems, we rather make use of the finite element approach.

For the mechanically regularized approach, we extract a region from the whole volume. The idea here is to use the ADDICT approach only on regions where highly localized phenomena occur. We choose to extract a slice from the volume of the reference configuration I r where the grey-level residuals were higher. The full cylindric slice has a volume of 17.865 million of voxels (a radius of 240 voxels and a height of 100 voxels). As it was indicated in chapter 1, when extracting the volume for meshing, a value smaller than the threshold value γ is assigned to the grey-level field at the boundaries so that we obtain a closed watertight surface. Figs. 5.9-5.10 show the closed extracted slice from the cylindric volume in addition to its finite element mesh. The mesh was generated using the CGAL library [The CGAL Project 2021]. The process took 30 minutes to be completed. 91 isolated connected components (noisy isolated voxels or groups of voxels) were removed. The considered mesh has 3 660 568 elements and 1 212 616 nodes (equivalent to 3 637 848 degrees of freedom). The average element size is around 0.34 voxels cube.

Once the finite element mesh is assembled, we apply the ADDICT approach by solving the mechanically regularized problem defined in Eqs. (3.57)-(3.58)-(3.59). It is worth noting that the mechanically regularized approach increases significantly the numerical cost of DIC as we add to the correlation hessian H the matrix that has the same pattern as K T K. The higher the regularization order, the denser the correlation linear system. Fig. 5.11 illustrates this remark by considering the finite-difference three dimensional Laplacian operator L of a scalar field. The number of extra diagonals increases when considering the Tikhonov scalar product of L T L instead of L only. This increases the resolution cost of the correlation system. This is one of the main reasons why we don't resort to the high order FCM method in this section, as we want to treat the largest possible region of interest in the limit of our available memory ressources.

From a general viewpoint, our implementation of ADDICT which consists in performing classical finite element assemblies over the whole domain, cannot be seen as the solution for treating very large problems. The DVC method based on macro elements and the ADDICT approach are two complementary approaches that should be used concurrently, for instance in a multiscale framework. In a first step, the first is used as the initialization of the second. The ADDICT approach was developed in the spirit of analysing regions where localized phenomena occur. That is why it would be interesting to find in the future, an automatic and efficient way to transit from the macro discretization to the micro discretization during the measurement process.

Let us apply now the ADDICT approach to the region of interest defined in Figs. 5.9-5.10. The idea is to apply the elastic equilibrium regularization scheme not on the whole image domain but only on the targeted slice. The ADDICT approach allows to remove all the voxels in the voids from the computation and as the considered foam is very porous (i.e. there is far more void than material), the correlation based on the ADDICT allows to treat only 3.6 million fine tetrahedral elements that represent the micro-architecture instead of the whole cylinder which volume (≈ πr 2 h) represents more than 17.8 million voxels.

After meshing the region of interest, the boundary condition regularization strategy must be defined. This is still the most tricky aspect for which no rigorous solution has been given in this manuscript. As we consider a compression test, and since the specimen is cylindric, the side surface of the slice should be free. In the experimental setup, the foam sample is put inside a rigid transparent polymer tube that plays the role of the load frame (see Fig. 5.1 right). Here, the internal diameter of the tube is the result of a compromise aimed on the one hand at minimizing the scan resolution (because of the cone beam technology, the X-ray source -sample distance must be minimized) and on the other hand at leaving a functional clearance (to get the sample move laterally). This solution has the merit of simplifying the setup of the reconstruction. If the lateral movements of the sample was not hindered, the side would be free. Fig. 5.12a illustrates the corresponding boundary conditions if the external boundary was free during loading. In reality, this hypothesis is not verified. Unfortunately, during loading, contact occurs between the loading cylinder and the external boundary of the foam sample. In order to avoid making the hypothesis of traction free in regions where it is not the case, we decide to penalize the nodes that belong to the external boundary in addition to the top and bottom nodes (as shown in Fig. 5.12c).

From the conclusion of chapter 4, the optimal regularization length l K in two-dimensions can be chosen equal to the cell length. As indicated at the very beginning of this chapter, we generalize this choice to the considered three-dimensional cellular geometry. By varying the regularization length l K , the minimal value that allowed the Gauss-Newton to converge is equal to 35 voxels which is quite near to the cell length. This re-confirms the characteristic length choice determined by the numerical experimentations performed in two dimensions. For the regularization length of the boundary given by l L , the choice performed in chapter 4 and which consisted in setting l L to its smallest value was not verified in this case. We proceeded by bisection with a large value of 50 voxels and gradually decreased it. l L was set equal to l K . The stiffness matrix K of the architecture was assembled with a Young's modulus E = 1 and an arbitrary Poisson coefficient ν = 0.28. In total, the mechanically regularized scheme started from the initial DVC solution obtained with the macro discretization shown in Fig. 5.5c (12 voxels element size and regularized with the diffusion model of Eq. (3.45)). This solution will be denoted as the macro solution in the following results and figures. The total time including loading the images, assembling the numerical operators, running the Gauss Newton iterations and post-processing is approximately equal to 12 minutes (the time for determining the initial guess is not included). For the initial guess measure, which corresponds again to the DVC measure with the finite element mesh of Fig. 5.5c, the total running time was equal to 8 minutes. In total, all the process consumed 55 GB of RAM and lasted 20 minutes. The first step for the analysis of the correlation consists in analyzing the grey-level residual distribution. We recall that it is defined by Eq. (5.1). r(x) = I r (x) -I d (x + u(x)), ∀x ∈ Ω r , with u being the measured displacement field.

(5.1)

To do so, we consider the histogram of the residual field in Fig. 5.13 along with the residual itself post-processed on the cellular architecture, see Fig. 5.14. First we observe that when the correlation converges (i.e. for the three first states), the residual distribution is a Gaussian centered at zero greylevels. We can see that the correlation accuracy decreases with the increase of the loading increments when macro-elements are used. This can be explained by the fact that the kinematics becomes more and more complex, but the discretization is the same. In parallel, we can also observe the grey-level residual field in Fig. 5.14 which confirms this statement. Fig. 5.13b shows that the uncertainty is divided at least by a factor of 2 for the first and second states and a factor of 4 for the third state. The mechanically regularized approach clearly improves the correlation accuracy for the third state. However, with multiple initialization attempts, the fourth loading state (which has a 28.5% strain level) could not be correlated correctly. This is detected with the second peak in the histograms and corresponds to the very dark circular region observed on the grey-level residual field. This failed correlation is due to the fact that the Gauss-Newton algorithm is very sensitive to the initialization and none of the proposed initialization techniques in chapter 3 succeeded. Another source of the problem can be related to the chosen cylindric mesh that contains the foam sample. If we look closely to the shape of the external boundary of the sample shown in Fig. 5.2, Fig. 5.9 or even Fig. 5.1, then we notice that the boundary is not completely straight. This induces that a large part of void was correlated in the elements which has perhaps increased the measurement uncertainty in this region. Let us focus only on the third loading state as it is a state with already a non-negligible strain level (12.3%). We display in Figs. 5.16-5.17-5.18 the obtained displacement fields in x, y, z directions (we recall that y is the compression direction). From a macro point of view, both solutions (micro and micro) are similar. However, if we look very precisely at some struts, displacement variations are measured at the subcellular scale with the micro approach which is not the case for the macro solution.

One of the main problems encountered with the micro measure is the irregular solution at some finite elements/nodes. The displacement seems to diverge to very large nonphysical values. This can be either explained by the mesh quality or by the fact that no grey-levels are integrated under the corresponding tetrahedrons. This problem illustrated in Fig. 5.15 is the main problem that alters the convergence of the optimization algorithm. Figure 5.17: Axial displacement u y .

In order to look for the differences between the two DVC approaches in more details, let us extract from the foam slice another small region of interest as indicated in Fig. 5.19. We compare the displacement fields in order to get a planar view of the cell region. The displacement fields on this It can be seen that the strain is homogeneous at the cell scale for the macro approach, while it becomes possible to clearly locate the spans that concentrate the largest deformations and bendings. The ADDICT method makes it possible to extract much more detailed and valuable information from the same raw data. Finally, in order to see even more clearly the differences between the macro and micro solutions, we proceed to a visualization approach similar to the one shown in chapter 4. As we can also segment the images of the deformed configurations, we can verify visually if the displacement field obtained by a DVC method allows to align the initial configuration on the deformed configuration. Fig. 5.24 illustrate this verification approach. We can see in the extracted region a large number of cell struts that undergo geometric non-linearity like buckling or large rotations. We also observe some kind of localized bands. The difference between the micro and macro DVC measurements is clear: the macro solution only allows rigid transformations at the cell scale and is completely unable to represent the complex local kinematic of each individual cell strut. This reconfirms again the same observations made in the two-dimensional study. Finally a zoom on some regions where buckling occurs is given in in order to appreciate the differences between the macro and micro solutions. These figures consist in superimposing, on the deformed image I d , the mesh constructed on the reference image I r , advected by the measured displacement field. It can be seen that the macro approach (global DVC but the result would be the same with the subset approach) does not allow to describe the particularly complex kinematics of each span. Here again, we can see that the mechanical regularization does not act as a strong regularization that constrains the kinematic measurement in the space of purely elastic solutions. Indeed, the measured displacement does not correspond at all to the elastic model. Only the high frequency component of the displacement (the displacements internal to the spans) are affected by the elasticity.

Conclusion

In the absence of texture at a scale smaller than the cell struts, the grey-scale conservation functional alone is unable to estimate local kinematic fields even roughly. Nor can it alone identify a strut (that localises in term of strains) more than others. However, this functional makes it possible to estimate the macroscopic component of the displacements provided that any sufficiently strong (subset or element size) or weak (low order Tikhonov like) regularization is used.

In the first two-dimensional study, we showed that it is possible to complement this macroscopic estimate obtained by the grey-level functional with an estimate at the microscale by relying weakly on an a priori assumption of the underlying physics. Although not limiting, the assumption used here was linear elasticity, even if the observed behaviour was non-linear. More interestingly, the regularization weighting parameter l K acts as a flexible way to separate the scales: the parts of the displacement of wavelength greater than l K are handled by the grey-scale metric (if sufficient gradients exist in the image) while the ones smaller than l K by the model.

We provided the numerical evidence that (a) the L-curve technique allows objective selection of the regularization parameters, (b) the optimal length coincides with the minimum of the true error and (c) the optimal length predicted with this technique is fully consistent with the lengths involved in the architecture of the material studied. The L-curve study can thus be avoided in practice to find a suitable parameter, since observations of the architecture of the material (with possible computation of the auto-correlation) may be sufficient as a first approach.

By studying numerous synthetic and real test cases, both linear and non-linear, and with the aim of producing, each time, a reliable reference to compare with, we have been able to show that this method provides reliable local information on the distribution of strains. Indeed, even if the reconstructed geometry does not perfectly match the actual specimen geometry, even if the behaviour of the regularization model is not perfect with respect to the true response of the architectured material at the sub-cellular scale (elastic vs. nonlinear), we have shown that the method allows to estimate complex local kinematic fields (displacements and strains) in a robust way in very poorly defined images and in the absence of texture. More than that, the method allows identification of the distribution of strains in the various struts and the zones within each strut where the strain localises, despite the poorly corresponding input data.

The three-dimensional extension of the ADDICT approach made it possible to measure particularly complex three-dimensional kinematics (such as buckling or bending) in the absence of pattern at the strut scale. The efforts on the implementation allowed to treat a real case representative (with several million degrees of freedom) of an in situ mechanical test. of a mechanically stressed cellular sample using DVC without texture at the sub-cell scale. Therefore our approach, consisting in exploiting the sample-specific model and the image gradient information contained in the domain boundaries by a grey-level dispersion type metric, appears as a promising strategy to investigate the mechanics of structured materials at their micro-structural scale. Such a tool should undoubtedly be valuable for studying the behaviour of a large number of cellular materials (metallic/polymeric foams, bones, wood, additively manufactured lattice structures, etc.). The work presented in chapter 5 will constitute the content of a paper and its writing is in progress.

Prospects

We list in the following some prospects and improvements:

• Numerical solving for the finite cell method When performing the convergence analysis under element and order refinement, an asymptotic behavior is observed which proves that too much refinement can lead to unnecessary computations in image-based modeling. One of the prospects of the observatory study performed in [START_REF] Rouwane | [END_REF]] would be to find rules (based on grey-levels, noise, porosity, resolution) that estimate a modeling error from which refinement criterias can be deduced automatically.

Concerning the numerical solving, the direct resolution of the obtained linear systems reaches rapidly its limits in three dimensions. Iterative solvers must be considered as it was suggested in [de Prenter et al. 2017, Jomo et al. 2017, Jomo et al. 2019]. They could contribute to treating large scale problems.

• High performance computing

The problem of high performance computing for regularized DVC is of double difficulty. First, from a memory point of vue, the equilibrium gap regularization increases the numerical cost of the DVC algorithm as the number of extra diagonals increases in the left hand side of the resolved linear system. The computational cost issue may become a real concern in three dimensions. Domain decomposition techniques or model reduction techniques particularly adapted to the tensor structure of B-splines could then be used advantageously [Passieux & Périé 2012, Gomes Perini et al. 2014, Bouclier & Passieux 2017]. In addition, adequate iterative solvers must be adapted for the special Gauss-Newton algorithm of DIC/DVC [START_REF] Tournier | [END_REF]].

• Adaptive isogeometric refinement for DVC Isogeometric hierarchical (local) refinement is a powerful technique that is very used in structural analysis [Giannelli et al. 2012,D'Angella et al. 2018,D'Angella et al. 2018,Garau & Vázquez 2018, D'Angella & Reali 2020]. When using B-splines, truncated hierarchical B-splines (THBS) were specially designed for getting analysis suitable spaces defined as the union of nested B-spline spaces which basis functions have all the desired properties for analysis (i.e. linear independance, partition of unity and minimal support). This refinement technique has been used in direct simulation and in the particular context of the finite cell method [Schillinger et al. 2012a, Verhoosel et al. 2015]. It has been also used in DIC in order to decrease measurement uncertainty in special targeted regions [START_REF] Kleinendorst | On the use of adaptive refinement in isogeometric digital image correlation[END_REF]]. The authors have constructed in this work a refinement indicator that allows to automatically adapt the mesh size in regions where the grey-level residual is high which allows to increase the measurement resolution without refining all the mesh. For the FCM version of the ADDICT approach, the idea would be to use a coarse mesh as usual when performing classical DVC (i.e. using macro elements). Using the refinement indicator defined in [START_REF] Kleinendorst | On the use of adaptive refinement in isogeometric digital image correlation[END_REF], the regions where the macro correlation fails to represent the kinematics can be detected. We then come enriching this region by refining the mesh and adding a suitable regularization scheme only where it is needed. This approach seems to be perfect for DVC as it allows to perform an easy and elegant numerical transition between micro and macro DVC. It can also be seen as a solution for approximating irregular displacement fields.

• Extension to multi-phase materials

Concerning the structural analysis aspect, in order to model multi-phase (> 2) materials, it can be interesting to investigate coupling methods [Nguyen et al. 2014, Mikaeili et al. 2022] and other advanced fictitious domain methods such as those based on the cutFem method [Kerfriden et al. 2020].

Concerning the measurement aspect, very recently, the equilibrium-gap regularization was used for DVC measurements in heterogeneous materials, for example in [Naylor et al. 2019] where a carbon fibre and epoxy matrix composite was studied and [Tsitova et al. 2021] where the damage in mortar was quantified.

As highlighted in the mechanically regularized optimization schem of Eq. (3.57), the Young's modulus E is set to 1 and the elastic effect is guided by the regularization parameter. We can get the E term out of the norm only because K is proportional to E (i.e. only one considered phase). In the case where a material with multiple mechanical properties is considered, for example with two material phases, the formulation (3.57) is not adapted because it introduces a competition between the two phases. Therefore, it is better in this case to consider as many meshes as phases and assemble independant stiffness matrices (especially if there are more than two phases). A first step would be to segment the different phases, determine the boundaries common between them. If two phases remain connected during the deformation, then a condition of equal displacements at the intefaces must be considered [Bouclier & Passieux 2017].

• Boundary regularization model

As higlighted in chapter 5, Fig. 5.15, multiple extreme irregularities were obtained at certain finite element nodes when considering a treahedral mesh that conforms to the architecture. This is either related to the quality of the constructed finite element mesh or to the limits of the regularization model at certain finite element nodes. In addition to the unclarity of the action of the boundary regularization defined by the restriction of the volume Laplacian operator, more rigorous investigations must be perfomed in order to find a solution to this problem (in particular investigate the boundary regularization approach suggested in [Mendoza et al. 2019]).

• Improvement of the regularization model Among other perspectives, a very interesting avenue concerns the regularization operator. It is indeed possible, with exactly the same formalism, to consider more advanced models (in particular non-linear ones) [Réthoré et al. 2013]. In particular, it would be interesting to update the constitutive parameters and the geometry stiffness of the regularization model, which is possible within the very same framework [Réthoré 2010, Réthoré et al. 2013].

• Surface correlation

To some extent, our strategy shares many similarities with the newly-introduced mechanical shape correlation [START_REF] Kleinendorst | Mechanical Shape Correlation: A novel integrated digital image correlation approach[END_REF] which is an integrated DIC approach based on the contours of deforming curvilinear objects. Similarities can also be made with the virtual image correlation which aims at measuring also complex curved boundaries of objects [START_REF] Réthoré | [END_REF], Jiang et al. 2019].

As no speckle is given under the cell-scale, one can actually view the correlation problem for cellular architectures as a problem of registering two surfaces (here defined by two level-sets or two surface meshes). It would be therefore interesting to use the finite cell embedding mesh (in a free form deformation (FFD) fashion [Chapelier et al. 2021]) as a measurement support to guide the initial configuration surface by using the boundary grey-level information in addition to the normal in order to align the boundary on the deformed configuration. This approach consisting in correlating two arbitrary surfaces (defined by a level-set of a set of polyhedrons) can be efficient from a numerical point of view as it allows to reduce the cost because the volume tetrahedrons are not taken into account as no information is available under the cell scale. A resulting question would be the dialogue between the obtained surface measurement and a FE computation that defines the behavior inside the cell struts.

• Initialization of the Gauss-Newton algorithm One of the most difficult aspects in image correlation is the accurate initialization of the gradientbased optimization scheme. Even though some numerical strategies were suggested in chapter 3, their use is quite manual and heuristic in the sense that multiple attempts are performed before running the final correlation. Another prospect of this work would be to investigate other robust initialization strategies so that minimal intervention is needed by the correlation user.

• Multiple acquisitions (resolutions)

In order to reduce the geometric errors related to noise and poor resolution (imposed by the imaging tool), it would be interesting to perform high resolution scans of the studied samples with high energy tomographs so that a reference geometry in which a high trust is placed can be set once for all. From this reference image, one can build faithful geometries for both the simulation and correlation aspects. For correlation, one can go back to conventional tomographs in order to perform in situ testings. A perspective to this work would be to find a robust framework (defined for example by three-dimensional projectors) that allows to make the dialogue between different acquisitions of the same object.

Appendices

Appendix A

Isogeometric Analysis: Basics and numerical implementation

This appendix is a brief summary of the isogeometric analysis (IGA) method applied to elasticity (which constitutes our home-made Python IGA library). It also details some computational aspects related to our implementation of IGA.

This appendix might be seen as another redundance but as we re-implemented the method, we have prefered writing this appendix for the sake of facilitating the understanding of the important concepts of IGA to the reader and show some important implementation details. A reference book for understanding B-splines, NURBS and the basic algorithms from a geometrical point of vue is the NURBS book [Piegl & Tiller 1996]. The basic algorithms for IGA: knot definition, B-spline functions (and their derivatives) evaluation, curve and surface definition, knot refinement, degree elevation are all defined in this book in the C programming language. For IGA, in addition to the book [Cottrell et al. 2009], a large number of references is available online. We also emphasize that a large number of researchers have already made the effort of implementing open-source IGA libraries such as PetIga [Collier et al. 2013], GeoPdes [Vázquez 2016], Nutils [van Zwieten et al. 2018] and tIGAr [Kamensky & Bazilevs 2019] to name a few. IGA was introduced in [Hughes et al. 2005] and later detailed in [Cottrell et al. 2009] as a new method for bringing closer the fields of computer aided design (CAD) and analysis. The idea of IGA is to use the same basis functions that generate the geometry for finite element analysis purposes. Many geometrical objects in CAD can be described using parametric curves, surfaces or volumes that are based on B-splines or Non uniform rational B-splines (NURBS). Fortunately, these set of basis functions have many mathematical properties that make them suitable for analysis. Furthermore, their high regularity provides high accuracy compared to classical finite elements. IGA has been applied to multiple fields. The introduction in [Borden et al. 2011] lists for example applications in turbulence, fluid-structure interaction, incompressibility, structural analysis, shells, phase-field analysis, large deformation with mesh distortion, shape optimization, topology optimization and immersed boundary methods which is its main application in this work.

A.1 B-splines: Definition and properties B-splines can be seen as the set of functions that generalize Bernstein polynomials. Bernstein polynomials are used for generating Bézier curves and emerge when applying a subdivision scheme based on piece-wise linear interpolation. Let us illustrate this through an example. Given four control points P 0 , P 1 , P 2 and P 3 in R 2 , we can generate a regular curve by successive linear interpolations (see

• p = 2                        N 0,2 (ξ) = 1 2 ξ 2 1 [0,1[ (ξ) + 1 2 ξ(2 -ξ) + 1 2 (3 -ξ)(ξ -1) 1 [1,2[ (ξ)+ 1 2 (3 -ξ) 2 1 [2,3[ (ξ) N 1,2 (ξ) = 1 2 (ξ -1) 2 1 [1,2[ (ξ) + 1 2 (ξ -1)(3 -ξ) + 1 2 (4 -ξ)(ξ -2) 1 [2,3[ (ξ)+ 1 2 (4 -ξ) 2 1 [3,4[ (ξ)
.

• p = 3 N i,p (ξ)N j,q (η)N k,r (ζ)P i,j,k (A.12)

with P i,j,k = (X i,j,k , Y i,j,k , Z i,j,k ) ∈ R 3 , where X, Y, Z ∈ R n×m×l are the control lattice coordinates.

All the notions that are defined for curves are directly applicable to surfaces and solids by tensor product generalization. NURBS surfaces and volumes are defined similarly to Eq. (A.10) by dividing respectively Eq. (A.11) and Eq. (A.12) by the linear combination of weighted B-splines. The tensor product nature of B-spline surfaces and solids implies that local refinement is impossible. That is why multiple techniques were developed in order to surpass this obstacle. We can list for example two interesting strategies: hierarchical B-splines [Forsey & Bartels 1988], truncated hierarchical Bsplines [Giannelli et al. 2012] and many other methods such as T-splines [START_REF] Sederberg | [END_REF].

A.2 Generation of some simple surfaces and volumes

We present in this section some very simple surfaces and volumes. In each figure of the geometries, the red points are the control points of the geometry, the green lines are the segments of the control grid and the black lines are the boundaries of the physical elements (i.e. the mapping of the parametric elements).

A.2.1 Rectangular plate

A rectangle is the most trivial geometry that can be generated using IGA. In order to get a plate defined by C p-1 elements and discretized with n x × n y elements, the idea is to start from one linear Q 4 element, perform degree elevation up to the desired degree p and then perform knot insertion (by This geometric parametrization leads to obtaining multiple coincedant nodes with the same coordinates. This creates a geometric fracture. As the treatment of identical degrees of freedom was not taken into account in our computational code, we avoid this irregularity by using another parametrization that was suggested in [START_REF] Kiendl | Isogeometric analysis and shape optimal design of shell structures[END_REF]] and that consists in creating a disk by degenerating a square. The degenerated disk of center (x 0 , y 0 ) and of radius r is defined as the following quadratic patch: ξ ξ ξ = (0, 0, 0, 1, 1, 1), η η η = (0, 0, 0, 1, 1, 1), (A.23) We should note that Nurbs are not practically used. Usually, circle arcs are well approximated using just B-splines (see for example [Apprich et al. 2017] where the authors approximate well a Nurbs circle using a non-rational cubic B-spline).

X =   -1 -2 -1 0 0 0 1 2 1   r + x 0 , Y =   -1 0 1 -2 0 2 -1 0 1   r + y 0 , W =   1 1/ √ 2 1 1/ √ 2 1 1/ √ 2 1 1/ √ 2
A.2.6 Quarter of a Torus: extrusion along the Nurbs circle arc.

Let us consider now a more complex Nurbs solid geometry. In general, B-spline solids are not immediate to generate using manual control (as in manual free-form design tools), therefore they are generaly defined as a geometric transformation of a B-spline surface (for example through extrusions along a certain curve). In order to generate a cylindric volume for example, the idea is to extrude either the singular disk or the degenerated disk.

A.3 Isogeometric analysis applied to linear elasticity

After giving a brief summary of the major aspects of B-spline modeling, we will see in this part how to use the same basis functions that define the geometry for the analysis of a partial differential equation. We will treat examples from linear elasticity applied to planar surfaces but the methodology remains the same for any boundary problem and can be applied to other fields. The main idea of IGA is to consider the B-spline basis functions as a basis for the space of approximation. Given a weak formulation of a linear partial differential equation of the form a(u, v) = l(v) where u is a solution that depends of the spatial variable x, the approximation of u at a physical point x(ξ), denoted by u h is given by the linear combination given in (A.25). where m is the number of basis functions which is equal to the number of control points and u(p i ) is the displacement of the control point. We note that the number of degrees of freedom is the number of basis functions multiplied by the dimension of the vector field. The starting point of the IGA program is a B-spline surface of degree (p, q), defined by two knot vectors ξ ξ ξ = (ξ 0 , ..., ξ r ) and η η η = (η 0 , ..., η l ). We note m the total number of basis functions which is equal to the number of the control points of the surface. By definition (A.4), m is equal to (rp) × (lq). The B-spline surface is defined by its control points and weights X, Y, Z and W that are matrices in R (r-p)×(l-q) . For the solid case, they are control tensors (of order three) but we only focus on two-dimensions as the generalization to three-dimensions is exactly the same. We define p x , p y , p z , p w ∈ R m the vectors of x, y ,z coordinates of the control points in addition to their weights. We define, as it is common in the finite element method, the vector of unknown displacements u = (u x 1 , ..., u xm , u y 1 , ..., u ym , u z 1 , ..., u zm ) T which represents the displacement of the control points vector (p x , p y , p z ). The following two equalities of (A.26) are the fundamental relations of IGA. They express the definition of the mapping between the parametric space (Ω p ) and the physical space Ω: The main difference with the finite element method is that a node does not belong to the physical geometry but rather controls it. In finite elements, the geometry is constructed using the basis of approximation of the solution. IGA, however, constructs the approximation space that fits to the geometry. Let us consider the problem (A.27) of plane stress linear elasticity without volume forces defined on a domain Ω: The operator D in (A.28) is called the stiffness operator.

                     ∇.σ = 0 in Ω ε = 1 + ν E σ - ν E tr (
The isogeometric analysis is founded on the Galerkin method which consists in transforming the continuous problem into a discrete problem using a weak formulation. Using Green formula, it can be shown that the weak formulation of problem (A.27) is given by (A.29) (see for example [Ern & Guermond 2002]). Solving the weak formulation consists in finding the displacement u that satisfies the dirichlet condition for all the v's that are also in the same admissible space [Ern & Guermond 2002]. The product σ(u), ε(v) F is the Frobenius inner product defined by A, B F = tr(A T B) for two matrices A and B and tr is a the trace of a matrix. The equation (A.29) represents the principal of virtual work. The linear form l(v) is interpreted as the work of the external forces in the virtual displacement v. The bi-linear form a(u, v) represents the work of the constraint at the equilibrium state in the deformation ε(v). The quantity j(u) = 1 2 a(u, u)l(u) represents the potential energy of the mechanical system. 1 2 a(u, u) is the energy of deformation and l(u) is the work of external forces.

Using the approximation in the B-spline space, (A.29) can be reduced to solving a linear system expressed by (A.30). .30) where u is the vector of displacement at the different control points of the geometry. The number of degrees of freedom denoted ndof is equal to the number of control points multiplied by the dimension of the problem. The stiffness matrix K and the force vector f are assembled using the contribution of each element. The next section details the computational aspects related the implementation of the IGA method for linear elasticity.

Ku = f , K ∈ R ndof ×ndof , f ∈ R ndof (A

A.4 Plane stress implementation

The general form of Hooke's law is defined by Eq. (A.31).

σ =          σ xx σ yy σ zz σ xy σ xz σ yz          = H          ε xx ε yy ε zz 2ε xy 2ε xz 2ε yz          (A.31)
where H is the Hooke stiffness matrix defined by:

H = E (1 + ν)(1 -2ν)          1 -ν ν ν 0 0 0 ν 1 -ν ν 0 0 0 ν ν 1 -ν 0 0 0 0 0 0 (1 -2ν)/2 0 0 0 0 0 0 (1 -2ν)/2 0 0 0 0 0 0 (1 -2ν)/2          (A.32) and          ε xx ε yy ε zz 2ε xy 2ε xz 2ε yz          =          ∂ x 0 0 0 ∂ y 0 0 0 ∂ z ∂ y ∂ x 0 ∂ z 0 ∂ x 0 ∂ z ∂ y            u x u y u z   .
(A.33) Formula (A.33) is a rewriting of ε = 1 2 (∇u + ∇u T ) in a vector form. In this section, we make the plane stress hypothesis (i.e. σ zz = σ xz = σ yz = 0). In this case, we have the following relation: .34) where H is the Hooke matrix given by (A.35) for an isotropic material: .35) The relation between strain and displacement at a physical point (x, y) is expressed by the equality (A.36) where u x and u y represent respectively the displacement in x and y directions. In order to discretize the forms in Eq. (A.29), we need to introduce some differential matrices that are very useful. Before that, let us note that as the approximation of the displacement at a point (x(ξ, η), y(ξ, η)) is given by the B-spline linear combination (A.26), it is better to define N w i the bi-variate rational function corresponding to the control point of index i. N w i is defined as follows: In this case, the displacement is expressed as follows: u(x(ξ, η), y(ξ, η)) = N(ξ, η)u.

H = E 1 -ν 2   1 ν 0 ν 1 0 0 0 1 2 (1 -ν)   . (A
N w i (ξ, η) =
We also need to define the strain differential operator B defined by: When considering the approximation in the B-spline space, it is important to keep in mind that the derivatives in (A.40) are in the physical space and the basis functions are defined in the parametric space. It is therefore necessary to compute the spatial derivative by inverting the Chain rule:

∂ ξ ∂ η = ∂ ξ x ∂ ξ y ∂ η x ∂ η y ∂ x ∂ y (A.41)
Inverting Eq. A.41 leads to: .42) After defining the different operators, we can now establish the stiffness system. Starting from the weak formulation (A.29), we are looking for a displacement u at a point (x, y) which is a linear combination of the B-spline basis functions. We start by writing the weak formulation:

∂ x ∂ y = 1 ∂ ξ x∂ η y -∂ ξ y∂ η x ∂ η y -∂ ξ y -∂ η x ∂ ξ x ∂ ξ ∂ η . ( A

A.4.1 Assembly

In elemental assembly, the parametric element should be be mapped to the reference element on which the position and the weights of the quadrature points are known. The mapping between the isoparametric element and a parametric element is given by:

ϕ : [-1, 1] × [-1, 1] -→ [ξ i , ξ i+1 ] × [η j , η j+1 ] (u, v) -→ ϕ(u, v) = ξ i + (ξ i+1 -ξ i ) u+1 2 η j + (η j+1 -η j ) v+1 2 . (A.50)
The integral on a parametric element is therefore written as follows: Let E be the matrix that indicates the indices of the basis functions that do have a non zero support on en element e. For a given index element e, E e,i is the global index in the vector u of the i-th basis where f is the evaluation of the source term at the integration points. Using this approach, the stiffness operator K can be assembled similarly by a vectorization operation. The discretization of the bi-linear form a can be defined as follows: Finally the discrete version for the bi-linear form is written:

a u x u y , v x v y = v T Ku (A.66)
with K is the stiffness matrix defined by:

K = H 0,0 × ∂ x N x T W g ∂ x N x + H 0,1 × ∂ x N x T W g ∂ y N y + H 0,2 × ∂ x N x T W g (∂ y N x + ∂ x N y ) + H 1,0 × ∂ y N y T W g ∂ x N x + H 1,1 × ∂ y N y T W g ∂ y N y + H 1,2 × ∂ y N T y W g (∂ y N x + ∂ x N y ) + H 2,0 × (∂ y N x + ∂ x N y ) T W g ∂ x N x + H 2,1 × (∂ y N x + ∂ x N y ) T W g ∂ y N y + H 2,2 × (∂ y N x + ∂ x N y ) T W g (∂ y N x + ∂ x N y ).
(A.67) This method of assembly has a practical advantage because it does not require to perform loops over elements but only requires using simple dot products. These operations are efficiently implemented in programming languages like Matlab or Python. Once the differential operators are created, this allows to integrate any function in an easy way. We should note that the matrices N and its derivatives can be constructed with any type of basis functions. The only drawback of this method is that we are saving in memory the differential sparse matrices and this becomes very memory consuming for high order (p > 1) three-dimensional problems that is why we only use this method in two dimensions and sometimes in three dimensions for linear finite elements. However, all this assumes that we have an efficient and fast method for constructing the differential operators N, mathbf partial x N and ∂ y N. In the case of finite elements, the definition of basis functions on the iso-parametric element is performed once for all and the construction of these sparse differential operators can also be performed in a vectorized way (see for example the integration routines in the Pyxel library [Passieux 2018]). In the case of B-splines they can be also be obtained simply by tensor product of univaritate operators. Roughly written N = N η ⊗ N ξ .

A.4.2 Imposing boundary conditions

In the contrary to finite elements where the nodes of the FE mesh define the geometry, imposing Dirichlet/Neumann boundary conditions is very complex in IGA. In fact, up to this date, we use C -1 regularity at the boundaries by considering open knot vectors (first and last knot vectors have a multiplicity equal to the degree + 1). In this case, the boundary control points of the surface are conforming to the physical surface and therefore we can directly apply the loadings over these control points. However, if we want to impose conditions on regions on the physical surface that are non conforming to the control mesh, then it is not possible to apply strong boundary conditions. A solution would be to perform inverse knot insertion which consists in performing knot insertions so that the obtained control point conforms to the physical point on which we would like to impose a boundary condition. In the context of this work, boundary conditions are directly imposed on conforming control points. In all examples, the loadings are applied on the image of one of the four boundaries of the parametric space (i.e. ξ = ξ min and η varies or ξ = ξ max and η varies or η = η min and ξ varies or η = η max and ξ varies). Concerning how boundary conditions are imposed, we use in this work two methods:

• Lagrange multipliers method:

The idea, is that the Dirichlet boundary conditions can be written by introducing a boundary condition matrix denoted here C. We can group the boundary displacement relations by writing:

Cu = u d .
(A.68)

Imposing the Dirichlet boundary conditions using the Lagrange multipliers method consists in solving the following system:

K C T C 0 u λ = f u d (A.69)
where the number of rows of C is the number of boundary conditions and the number of columns of C is the total ndof . 0 is a square matrix of size equal to the number of boundary conditions. Afterward, the values of λ are discarded. This comes from the solution of a constrained quadratic optimization problem which consists in minimizing the discretized energy under the Dirichlet boundary condition constraint:

arg min

u∈R ndof Cu=u d 1 2 u T Ku -u T f . (A.70)
It can be shown using Lagrange duality that the solution is equivalent to solving: arg min

u∈R ndof λ 1 2 u T Ku -u T f + λ T (Cu -u d ).
(A.71)

The system (A.69) represents the optimality condition of Eq. (A.71) with respect to (u, λ).

We have to note that type of boundary conditions makes the final system not positive. Thus the resolution with the Cholesky decomposition does not remain possible. To keep a symmetric positive definite system, we rather proceed to the strategy that consists in solving the equilibrium only on internal nodes. A reference strain energy of problem (A.22) is given by (A.74) and can be found in [Zienkiewicz et al. 2005].

E ex = 1 π (log(2) -0.6) (A.74)
We can see in Fig. A.24 that the theoretical order of convergence O(h p ) is obtained. The relative error curves have a slope equal to minus the order of approximation. For a fixed degree and for the same number of degrees of freedom, the isogeometric analysis solution is more accurate than the finite element solution. This comes to the end of this appendix. We have summarized in this appendix, the IGA method, its numerical implementation and showed two basic numerical examples that were used to validate the developed computational code. We will do our best to share this library, in addition to the developed DVC code (as open source) in the near future.

Figure 1

 1 Figure 1.2: Linear B-spline basis functions on the image domain [x min , x max ].

Figure 1

 1 Figure 1.3: Quadratic B-spline basis functions on the image domain [x min , x max ].

Figure 1

 1 Figure 1.4: Cubic B-spline basis functions on the image domain [x min , x max ].

Figure 1

 1 Figure 1.5: Ĩ when setting c = I.

  (a) Quadratic approximation.(b) Cubic approximation.

Figure 1

 1 Figure 1.6: Smoothing L 2 projection.

  dx (Partition of unity property of B-splines).

  Figure 1.8: Effect of gaussian filtering.

  Figure 1.10: Comparison between the bilateral and gaussian filter.

  Fig. 1.11 illustrates the different steps of the marchinq cubes algorithm in two dimensions. (a) Image field. (b) One moving window. (c) Zoom on the moving window. (d) Evaluation of the scalar field at the four nodes. (e) Determining of the interface positions by linear interpolation.(f) Linking the two linearized intefaces to obtain the edge contour.

Figure 1 .

 1 Figure 1.11: Illustration of the different steps of the marching squares algorithm.

Figure 1 .

 1 Figure 1.14: Iso-values of the continuous spline representation of the star image of Fig. 1.13.

  (a) Marching cubes without Gaussian smoothing. (b) Marching cubes applied with Gaussian smoothing (σ = 1).

Figure 1 .

 1 Figure 1.15: Marching cubes on a celullar architecture.

Figure 1 .

 1 Figure 1.16: Extraction of a region from a slice of a closed cell foam. Left image size: 934 × 629. Right image size: 126 × 140.

  (a) Grey-level image (voxel size is equal to 78µm). (b) Hexahedral finite element mesh. (c) Tetrahedral finite element mesh.

Figure 2

 2 Figure 2.1: Finite element meshes generated from a MRI image of a bovine trabecular (or cancellous) bone.

Figure 2

 2 Figure 2.3: High displacement value observed at a degenerated element. Left: u z displacement map. Middle: Two critical elements. Right: Deformed mesh.

  (a) Unsmooth surface. (b) Surface with coefficients c = G σ * I). (c) Facet size equal to 1.1 voxels. (d) Facet size equal to 0.9 voxels. (e) Facet size equal to 0.8 voxels. (f) Facet size equal to 0.4 voxels.

Figure 2

 2 Figure 2.4: Mesh generation of a polyurethane open cell foam specimen.

Figure 2

 2 Figure 2.5: Ω p is embedded in a rectangular domain. Image inspired from [Schillinger & Ruess 2015].

  (a) Complete integrand on an element. (b) Penalized integrand on an element.

Figure 2 . 7 :

 27 Figure 2.7: Trimming of a basis function.

Figure 2 . 8 :

 28 Figure 2.8: Cut elements are recursively decomposed until the maximum level is reached (l max = 5) (Image inspired from [Schillinger & Ruess 2015]).

Figure 2 .

 2 Figure 2.10: Triangulation of the last cut cells of the quadtree scheme.

  (a) Evaluation of the level-set function at the cell vertices. (b) Linear estimation of the cut interfaces.(c) Tetrahedralization of the set of retained vertices.

Figure 2 .

 2 Figure 2.11: Tetrahedralization of a cut-cell.

  their tensor product in order to get the integrated block matrix over Ω e = Ω ex × Ω ey × Ω ez (see the illustration of Fig. 2.14).

Figure 2 .

 2 Figure 2.13: Integration of two-dimensional homogeneous cells in the quadtree scheme.

Figure 2 .

 2 Figure 2.14: Integration of three-dimensional cell cells in the octree scheme.

  An example of the embedding B-spline elements with the corresponding boundary conditions.

Figure 2 .

 2 Figure 2.16: Mechanical problem definition: elastic plate with a quarter hole.

  Triangulation of the last cut cells.

Figure 2 .

 2 Figure 2.18: Effect of the integration of the last cell on the convergence under h refinement.

  Figure 2.19: h refinement for different degrees of B-splines (along with the linear triangular finite element method).

  (a) Finite element discretization. (b) Finite cell method based on a cubic B-spline discretization.

Figure 2 .

 2 Figure 2.21: Boundary conditions for the curved beam problem.

Figure 2 .

 2 Figure 2.22: Curved beam under shear displacement: geometry and boundary conditions.

  (a) Finite cell method solution.(b) Boundary conforming solution.

Figure 2 .

 2 Figure 2.23: Von Mises stress field defined by (2.23).

Figure 2 .

 2 Figure 2.24: Summary of the different steps of the construction of a mechanical image-based model.

  2.25, left and middle). We construct an image by assigning a maximum grey-level value denoted v max to all the pixels that are completely in the domain [0, c[ and a minimum value denoted v min to those in the domain ]c, l]. The pixel that cuts the discontinuity of the step function is assigned with a mean grey-level value weighted by the intersection lengths. Fig.2.25 (middle and right) shows the procedure that consists in constructing the smooth representation of the image.

Figure 2 .

 2 Figure 2.25: Steps for constructing the synthetic image and the use of a smooth representation of the image.

Figure 2

 2 Figure 2.26: The different interfaces of the geometry approximation.

Figure 2

 2 Figure 2.27: Evolution of the intrinsic geometry error in pixel size unit when translating the true geometric interface all over one pixel.

  12 pixels in the image.

Figure 2 .

 2 Figure 2.28: Intrinsic geometry error E displayed along with the evolution of the errors Ē and Ẽ

  30 × 30 pixels in the image. 60 × 60 pixels in the image.

Figure 2 .

 2 Figure 2.29: Intrinsic geometry error E displayed along with the evolution of the errors Ē and Ẽ with respect to the size of the smallest sub-cell in pixel size units (2D open-hole plate test case).

  Mechanical properties based on gray levels Marching squares algorithm Q4 Voxel based model (a) 30 × 30 pixels in the image. Mechanical properties based on gray levels Marching squares algorithm Q4 Voxel based model (b) 60 × 60 pixels in the image.

Figure 2 .

 2 Figure 2.30: Evolution of the error in energy norm under mesh refinement for the elastic plate with a quarter hole.

Figure 2 .

 2 Figure 2.31: Geometry of a metallic plate with multiples circular holes. An elementary part is framed with the red color.

Figure 2 .

 2 Figure 2.32: Mechanical problem: micro-structure with a hexagonal structure

  (a) 50 × 43 pixels. (b) 100 × 86 pixels. (c) Refined finite element solution.

Figure 2 .

 2 Figure 2.33: Finite cell method using the level-set geometric representation. Cubic B-spline mesh corresponding to Fig. 2.32b.

Figure 2 .

 2 Figure 2.34: Marching squares algorithm solution.

Figure 2 .

 2 Figure 2.35: Q 4 voxel based model

Figure 2 .

 2 Figure 2.36: Mechanical properties based on grey-levels. Post-processed stress field on the complete mesh.

  (a) Volume of interest (Compression of the top red part and clamping of the bottom black part). (b) Finite element mesh. (c) Coarse FCM cubic discretization.

Figure 2 .

 2 Figure 2.37: Compression test on a polyurethane foam geometry. The green points represent the Dirichlet control points and the red points represents the removed basis functions.

  (a) FEM solution. (b) FCM solution.(c) Color scale.

Figure 2 .

 2 Figure 2.38: Transverse displacement u x .

Figure 2 .

 2 Figure 2.39: Axial displacement u y .

Figure 2 .

 2 Figure 2.40: Transverse displacement u z .

Figure 2 .

 2 Figure 2.41: Von Mises stress field response for the compression test.

  Figure2.42: Synoptic view of the proposed approach to Digital Image Correlation (DIC) to perform displacement measurements in cellular materials below the cell scale. The images do not exhibit any texture under the cell scale and the DIC problem is therefore poorly conditioned. Inspired from[Réthoré et al. 2009, Leclerc et al. 2011], we propose here to use mechanical regularization to obtain sub-cellular displacement fields. In this work, a stiffness matrix that accounts for the underlying architecture is used. An automatic image-based mechanical model, based either on the finite-element method or the finite-cell method, is considered to construct the latter.

  approximated by ∇I r [Neggers et al. 2016, Passieux & Bouclier 2019]: H S (u) = ∇I r ⊗ ∇I r (3.21) and ∇ S (u) = -(I r -I d • (id + u))∇I r (3.22)

  .26) Box 3.1: Resolution of the correlation problem.

  x) T ∇I r (x) ⊗ ∇I r (x)N(x)dx. (3.30) Similarly, b k is obtained from the assembly of block vectors b k e = Ωe (I r (x) -I d (x + N(x)u e )) N(x) T ∇I r (x)dx.(3.31)

  (a) Uniform subdivision of a structured element (Trilinear finite element or arbitrary order B-spline element). (b) Uniform subdivision of the reference tetrahron in [0, 1] × [0, 1] × [0, 1].

  (a) Boundaries of a domain Ω. (b) Dirichlet nodes for the finite element method.(c) Dirichlet nodes for the finite-cell method.

Figure 3

 3 Figure 3.3: Degrees of freedom removal for the equilibrium gap regularization (two-dimensional case).

  Dirichlet nodes for the finite-cell method.

Figure 3 . 4 :

 34 Figure 3.4: Degrees of freedom removal for the equilibrium gap regularization (two-dimensional case).

Figure 3

 3 Figure 3.5: Illustration of the multi-level refinement strategy based on image resolution and mesh variations. The texture represents a Perlin noise.

Figure 3 . 6 :

 36 Figure 3.6: Different possibilities for displacement exchange between two meshes.

Figure 3 . 7 :

 37 Figure 3.7: Axis aligned bounding boxes for a three-dimensional mesh. Image taken from the documentation website of the CGAL library [The CGAL Project 2021].

  Figure 4.1: Synthetic image generation and procedure to assess the performance of the DIC measurements.

Figure 4

 4 Figure 4.2: Definition of the virtual experiment: FE mesh of the exact geometric object displayed with the boundary conditions. The sample corners are defined by x min = 0 mm, x max = 110 mm, y min = 0 mm, y max = 50 mm.

  (a) Elastic simulation (tensile test u 0 > 0). (b) Elasto-plastic simulation (tensile test u 0 > 0).(c) Elasto-plastic simulation with non-linear geometric analysis (compression test u 0 < 0).

Figure 4 . 3 :

 43 Figure 4.3: Evolution of the resultant of reaction forces at the right end of the specimen with respect to the prescribed displacement u 0 in x direction: (a) linear elasticity test (i), (b) elasto-plastic tension test (ii) and (c) geometric non-linear elasto-plastic compression test (iii). The red dots represent the mechanical states used to generate the deformed images.

  (a) Image of the reference configuration I r (load-free). (b) Image of the deformed configuration for the elastic model subjected to tension corresponding to Fig. 4.3a. (c) Image of the deformed configuration for the elasto-plastic model subjected to tension corresponding to Fig. 4.3b. (d) Image of the deformed configuration for the geometrically non-linear elasto-plastic model subjected to compression corresponding to Fig. 4.3c.

Figure 4 . 4 :

 44 Figure 4.4: Example of pairs of DIC test images based on the same sample but with different mechanical models. Image dynamic range is equal to 255 in the whole image area and equal to 127 in the cell area only.

Figure 4 . 5 :

 45 Figure 4.5: Zoom on an image area. The finite element mesh is superimposed on the image. Displayed points are the Gauss integration points of the reference triangular FE mesh used for the computation of the error.

Figure 4 . 7 :

 47 Figure 4.7: Main steps to build the specimen-specific, immersed B-spline image-based model.

Figure 4

 4 Figure 4.8: Cubic B-spline grid taken to discretize the measured displacement field for the considered two-dimensional cellular-like specimen. The white region represents the studied cellular specimen. The black region is the image background.

Figure 4

 4 Figure 4.9: Retained B-spline control points to describe the mechanically regularized DIC solution for the considered two-dimensional cellular-like specimen.

Figure 4 .

 4 Figure 4.10: DIC inputs resulting from a virtual compression test.

  (a) Euclidean norm of the displacement field (10 -3 mm).(b) Von Mises strain norm.

Figure 4 .

 4 Figure 4.11: Exact fields corresponding to the deformation configuration of Fig. 4.10.

  4.12a), large disparities are corrected gradually in order to capture more local displacement fields. The solution of the multi-level correlation is displayed in Fig. 4.13. (a) 2 0 × 2 0 elements. (b) 2 1 × 2 1 elements. (c) 2 4 × 2 4 elements.

Figure 4 .

 4 Figure 4.12: Multi-level refinement scheme using cubic B-splines. At each step, the orange grid represents the B-spline grid defined over the region of interest.

Figure 4 .

 4 Figure 4.13: Solution of the global method without any penalization using a multi-level refinement strategy.

  4.13 in comparison to Fig. 4.11).

  (a) Euclidean norm of the displacement field for λ T = 10 1 . (b) Von Mises strain norm for λ T = 10 1 . (c) Von Mises strain norm for λ T = 10 3 . (d) Von Mises strain norm for λ T = 10 5 .

Figure 4 .

 4 Figure 4.14: Solution using the multi-level Tikhonov regularization.

  (a) Euclidean norm of the displacement field. (b) Von Mises strain norm.

Figure 4 .

 4 Figure 4.15: Solution using the equilibrium gap regularization.

  Reference image I r .

Figure 4 .

 4 Figure 4.16: Zoom on a region in the ROI.

  Evolution of the Von Mises strain uncertainty.

Figure 4 .

 4 Figure 4.17: Noise sensitivity analysis for fixed regularization parameters.

Figure 4 .

 4 Figure 4.18: Periodic rectangular signal.

Figure 4 .

 4 Figure 4.19: Auto-correlation of the periodic rectangular signal.

   displays the considered image and its auto-correlation function. In order to estimate the characteristic distances of the image, we consider the radially averaged normalized auto-correlation function displayed in Fig.4.20c . The first minimum is approximately equal to 12 pixels which is the diameter ∆ 1 . The first maximum of the auto-correlation function has an abscissa equal to 23 pixels which approximates the distance ∆ 2 . (a) Image of a grid of disks. (b) Auto-correlation function. (c) Radially averaged normalized auto-correlation function.

Figure 4 .

 4 Figure 4.20: A two-dimensional application of the radially averaged normalized auto-correlation function.

Figure 4 .

 4 Figure 4.21: Radially averaged normalized auto-correlation function.

Figure 4 .

 4 Figure 4.22: Necessary discretization for the standard subset DIC. The measurement points are marked by the red dots. A large part of the boundary subsets are automatically removed in order to avoid high uncertainty measurements in theses zones. The orange square depicts the subset size.

Figure 4 .

 4 Figure 4.23: Horizontal component u x of the displacement field in the ROI of the subset method (in pixel units).

Figure 4 .

 4 Figure 4.24: Distribution of the root mean square of the relative longitudinal and transverse displacement error.

Figure 4 .

 4 Figure 4.25: Plot of the equivalent Von Mises strain field ε vm = ε 2 xx + ε 2 yy + 2ε 2 xy .

4. 3 .

 3 Numerical investigation of the influence of the model and parameters used for the regularization 106 (a) Variation of l L . (b) Variation of l K .

Figure 4 .

 4 Figure 4.26: Elastic regularization versus data fidelity for ADDICT on an elastic problem.

  (a) Variation of l K for the proposed ADDICT. l L is fixed to its optimal value. (b) Variation of λ = l T 4 (see Eq. (3.41) with T = L T L) for the Laplacian-based regularization.

Figure 4 .

 4 Figure 4.27: Influence of the regularization lengths for the three input models. M1: Elastic model (i), M2: Elasto-plastic model (ii), and M3: Geometrically non-linear elasto-plastic model (iii). There are missing data for the Laplacian-based regularization over M3 because the convergence of the minimization algorithm was not achieved on these points with our residual tolerance and maximum number of iterations allowed.

Figure 4 .

 4 Figure 4.28: Influence of the regularization parameter on the root mean square displacement error (U (u x ) 2 + U (u y )) 2 /2.

Figure 4 .

 4 Figure 4.29: Superposition of the deformed point clouds using the reference finite element field (red point cloud) and the measured field using the equilibrium gap method (green point cloud). Figures corresponding to the elasto-plastic problem (ii).

4. 3 .

 3 Figure 4.30: Deformed mesh displayed on the deformed configuration I d when considering the Levenberg-Marquardt method.

Figure 4 .

 4 Figure 4.31: Superposition of the deformed point clouds using the reference finite element field (red point cloud) and the measured field using the equilibrium gap method (green point cloud). Figures are corresponding to the geometrically non-linear elasto-plastic problem (iii). (The point clouds are amplified with amplification factor of 2).

  Measure using the elastic regularization. (c) ε xx range values.

Figure 4 .

 4 Figure 4.32: ε xx strain.

  (a) Texture free side for ADDICT. (b) Textured side for reference FE-DIC.

Figure 4 .

 4 Figure 4.34: Experimental setup for the speckle deposition.

Figure 4 .

 4 Figure 4.35: Printed speckle with a 300 dpi resolution.

Figure 4 .

 4 Figure 4.36: The components of the air-brush: an efficient tool for speckle deposit.

  (a) Basic optical setup. (b) Overview of the experimental setup.

Figure 4 .

 4 Figure 4.37: Experimental setup.

Figure 4 .

 4 Figure 4.38: Sample state at the last loading increment. The speckle was deposited on the matt white paint after washing the sample with water and of course drying it.

  (a) Image of the un-textured face provided to the ADDICT. Image resolution: 4.5 pixels per mm. Definition of the sub-image presented: 88 × 73. (b) Image of the textured face provided to the FE-DIC. Image resolution: 78 pixels per mm. Definition of the sub-image presented: 1218 × 1558.

Figure 4 .

 4 Figure 4.39: Zoom on a specific region of the sample.

Figure 4 .

 4 Figure 4.40: Experimental force ( F )-displacement ( Ū ) curve. Unloads were performed to highlight the non-linear nature of the deformation. Each point corresponds to the acquisition of images. The red one indicates the state that is analyzed in the remainder of this paper.

  4.45).

  (a) Grid and level-set used to perform ADDICT on the speckle-free face.(b) FE-DIC mesh used to measure the displacement field on the textured face.

Figure 4 .

 4 Figure 4.41: ADDICT (speckle-free face) and FE-DIC (textured face) discretizations.

Figure 4 .

 4 Figure 4.42: Influence of the regularization lengths for the experimental test-case. Variation of l K .

  4.45).

  (a) ADDICT measurements with low resolution speckle-free images. (b) FE-DIC measurement with high resolution textured images. (c) FE displacement field computed with Abaqus.

Figure 4 .

 4 Figure 4.43: Comparison of the longitudinal displacement fields u x (mm) measured with ADDICT (u 1 ), FE-DIC (u 2 ) and computed with Abaqus (section 4.1.1) for an imposed displacement Ū = 1.05 (Fig. 4.2).

Figure 4 .

 4 Figure 4.44: Comparison of the transverse displacement fields u y (mm) measured by ADDICT (u 1 ), FE-DIC (u 2 ) and computed with Abaqus (section 4.1.1) for an imposed displacement Ū = 1.05 (Fig. 4.2).

Figure 4 .

 4 Figure 4.45: Point cloud belonging to the intersection of the level-set geometry and the FE geometry.

Figure 4 .

 4 Figure 4.46: Relative displacement distance map between ADDICT (u 1 ) and FE-DIC measurements (u 2 ). The difference is scaled by the displacement Ū imposed to the grips.

Figure 4 .

 4 Figure 4.47: Measured equivalent Von Mises strain ε vm .

Figure 5

 5 Figure 5.1: Experimental setup of the in situ compression test using X-ray micro-tomography.Image taken from[START_REF] Pétureau | Stratégie de couplage expérimentation/modélisation dans les matériaux hétérogènes. Identification de propriétés mécaniques locales[END_REF]].

Figure 5 . 2 :

 52 Figure 5.2: Rendering of the reconstructed volumes at the different loading states. The left state is the reference configuration (at rest).

  (a) Evolution of the porosity of the segmented volume with respect to the threshold value. (b) Rendered surface for γ = 80. (c) Rendered surface for γ = 56.

Figure 5 . 3 :

 53 Figure 5.3: Effect of the threshold value on the porosity of the segmented volume.

  (a) Auto-correlation function of the image in the reference configuration. (b) Slice of the reconstructed volume stack.

Figure 5 . 4 :

 54 Figure 5.4: Determination of the approximate cell size using the radially averaged normalized auto-correlation function.

Figure 5 . 5 :

 55 Figure 5.5: Three finite element discretizations for treating the correlation problem.

Figure 5 . 6 :

 56 Figure 5.6: Coarse grid used for initialization.

Figure 5 . 7 :

 57 Figure 5.7: Deformed macro meshes for the different loading states.

Figure 5 . 8 :

 58 Figure 5.8: Displacement fields (in voxels unit) for 12 voxels macro-finite element DCV using the mesh defined in Fig. 5.5c.

Figure 5

 5 Figure 5.9: Position and definition of the region of interest considered in the ADDICT treatments.

Figure 5 .

 5 Figure 5.10: Extracted slice of the foam sample for the ADDICT measurements.

  (a) Structure of the Laplacian operator L. (b) Structured of the square of the Laplacian operator L 2 .

Figure 5 .

 5 Figure 5.11: Matrix structure of the second and fourth order regularization terms in three dimensions (obtained with finite-differences).

Figure 5 .

 5 Figure 5.12: Description of the boundary regularization applied on the extracted foam slice.

  (a) Correlation with macro elements. (b) Correlation with micro elements.

Figure 5 .

 5 Figure 5.13: Histogram of the residual field.

Figure

  Figure 5.14: Residual field

Figure 5 .

 5 Figure 5.15: Warped reference configuration with the micro displacement solution. Extreme irregular solutions at some finite element nodes are observed.

Figure 5 .

 5 Figure 5.18: Transverse displacement u z .

Figure 5 .

 5 Figure 5.19: Extraction of a region of analysis from the foam slice (Clipping in the direction (1,0,0) with the Paraview software).

  Figure 5.20: Transverse displacement u x .

Figure 5 .

 5 Figure 5.21: Axial displacement u y .

Figure 5 .

 5 Figure 5.22: Transverse displacement u z .

Figure 5 .

 5 Figure 5.23: Von Mises strain field ε vm = ||ε|| F .

Figure 5 .

 5 Figure 5.24: From top to bottom: 1. Segmented image I r . 2. Segmented image I d for the third loading state. 3. Warping the reference state with the solution coming from macro elements. 4.Warping the reference state with the solution coming from the micro solution.

  (a) Ω r (b) In grey color: Ω d . In red color: Ω r + u macro (Ω r ). (c) In grey color: Ω d . In red color: Ω r + u micro (Ω r ).

Figure 5 .

 5 Figure 5.25: A first region where local bending occurs.

Figure 5 .

 5 Figure 5.26: A second region where local bending occurs.

  (a) Ω r (b) In grey color: Ω d . In red color: Ω r + u macro (Ω r ). (c) In grey color: Ω d . In red color: Ω r + u micro (Ω r ).

Figure 5 .

 5 Figure 5.27: A third region where local bending occurs.

  Figure A.2: B-spline basis for ξ ξ ξ = (0, 1, 2, 3, 4) with p ∈ {0, 1, 2, 3}.

Fig

  Fig. A.2 and formula (A.4) show that the basis functions are equal to step functions when p = 0, to piece-wise linear functions for p = 1 and to piece-wise quadratic functions when p = 2 etc.Given a set of B-spline basis functions defined on a knot vector ξ ξ ξ, a B-spline curve C : R -→ R 3 is defined by the following linear combination that generalizes the successive linear interpolation subdivision shown previously in (A.3):

  Initial geometry with 1 × 1 elements. (b) Refined geometry with 2 × 2 elements. (c) 3 × 3 elements. (d) 10 × 10 elements.

Figure A. 11 :

 11 Figure A.11: Degenerated disk.

  Fig. A.12 shows the obtained geometries. In order to define a (a) Extrusion of the irregular disk. (b) Extrusion of the degnerated disk.

Figure A. 12 :

 12 Figure A.12: Solid Nurbs cylinder.

  (a) Initial beam geometry obtained from the extrusion of zero radius anulus. (b) Initial beam geometry obtained from the extrusion of a degenerated quadratic disk.

Figure A. 13 :

 13 Figure A.13: Minimal beam geometry.

  Fig. A.14 shows the mapping that exists between the parametric space and the physical space. If u h represents for example the displacement, then u h (x(ξ)) is the displacement at the physical point given by x(ξ) = i N i (ξ)p i where p i are the coordinates of the control points.

Figure A. 14 :

 14 Figure A.14: Mapping between the parametric space and the physical space.

  σ)I in Ω u = u 0 on Γ d σ.n = f on Γ n A.27), σ represents the Cauchy stress tensor, n is the normal to the domain Ω, ε is the strain tensor, ν is the Poisson coefficient, E is Young's elastic modulus, f represents the loading at the boundary of the domain Γ n , u 0 is an imposed displacement at the boundary Γ d and ∇ is the gradient differential operator. The relationship between the strain ε and stress σ on any point of Ω can be rewritten in the form (A.28) known as Hooke's law. σ = D(ε) with D(ε) = Eν (1 + ν)(1 -2ν)

  a(u, v) = l(v) u = u 0 on Γ d with a(u, v) = Ω σ(u), ε(v) F dx and l(v) = Γn f.vds (A.29)

  x u y + ∂ y u x ) 1 2 (∂ x u y + ∂ y u x )

p

  wi N i (ξ, η)

  η) ... ∂ y N w m (ξ, η) ∂ y N w 1 (ξ, η) ... ∂ y N w m (ξ, η) ∂ x N 1 (ξ, η) ... ∂ x N w m (ξ, η)

  where B is the matrix of derivatives of B-spline basis functions defined in (A.40), v a discrete vector of virtual displacements and N is the matrix of basis functions defined in(A.39). By injecting (A.44) in (A.43) we obtainΩ v T B T HBu dxdy = Γn v T .N T fx is to change the integration domain to the parametric space because the function B T HB is defined on the parametric space Ω p . We denote S the mapping between Ω p and Ω and γ the mapping between a parametric boundary denoted Γ p and one physical Neumann boundary:S : Ω p -→ Ω (ξ, η) -→ S(ξ,η) the control points that define one Neumann boundary Γ n . Finally the stiffness system is written as it follows: K = Ωp B T HB(ξ, η)|det(J S (ξ, η))|dξdη, is the jacobian matrix of the transformation S. Its determinant is given by det(J S ) = ∂ ξ x∂ η y -∂ ξ y∂ η x. The norm of the tangent vector to the Neumann boundary is given by ||γ′ (ζ)|| = x ′ (ζ) 2 + y ′ (ζ) 2 .In order to evaluate the integrals in (A.49), a univariate Gauss integration is used.

BB

  T HB(ϕ(u, v))|det(J S (ϕ(u, v)))| × ξ i+1ξ i 2 × η j+1η j 2 dudv.(A.51)A first method for assembling the stiffness matrix is to consider the elementary contributions by dividing the integral:T HB(ϕ(u, v))|det(J S (ϕ(u, v)))| × mesB T HB(ϕ(g i ))|det(J S (ϕ(g ig )))| × mes(e)4 dudv(A.52)where w g i is the Gauss quadrature weight associated to the Gauss point g i and mes(e) is the measure of an element. We denote B e the elementary differential matrix as it follows: η) ...∂ y N w me (ξ, η) ∂ y N w 1 (ξ, η) ... ∂ y N w me (ξ, η) ∂ x N 1 (ξ, η) ... ∂ x N w me (ξ, η)   ∈ R 3×2me (A.53)where m e = (p + 1) × (q + 1) is the number of nonzero basis functions on the element e and C e the elementary matrix defined byC e = nge i=1w g i B T e HB e (ϕ(g ij ))|det(J S (ϕ(g ig )))|.

  Properties of the plate deformation problem. Table A.1 shows that ε xx = σ xx /E, σ xx = f x and ν = -ε yy /ε xx are verified. Fig. A.17 shows the displacement field distribution on the plate.

  (a) Displacement field in x direction. (b) Displacement field in y direction.

Figure A. 17 :

 17 Figure A.17: IGA solution for the plate deformation problem with symmetric boundary conditions.

Figure A. 18 :

 18 Figure A.18: Solution fields of the problem of Fig. A.19a.

  Figure A.21: Solution fields of the problem of Fig. A.19b.

Figure A. 22 :

 22 Figure A.22: Description of the circular beam problem.

Figure A. 23 :

 23 Figure A.23: Deformed geometries.

  Number of degrees of freedomRelative error of the elastic energy

Figure A. 24 :

 24 Figure A.24: Evolution of the relative error on strain energy in log scale.

  .13) and the block

	The same previous two-dimensional procedure is performed and all the two-dimensional operators are
	pre-computed. Again as an example, suppose we want to compute
				Ωe	∂φ φ φ ∂x	⊗	∂φ φ φ ∂z	dxdydz,	(2.15)
	then we need to look up for the computed terms
				Ωe	∂φ φ φ ∂x	⊗	∂φ φ φ ∂x	dxdy	(2.13)
	is needed to be computed. With this factorization approach, we just look for the pre-computed
	integrated one-dimensional operators:
	Ωe x	dφ φ φ x dx	⊗	dφ φ φ x dx	dx and	Ωe y	φ φ φ y ⊗ φ φ φ y dy	(2.14)
	(for the x and y directions) corresponding to the adequat quad-tree level and perform their tensor
	product (see again Fig. 2.13). The idea is quite similar in three-dimensions as illustrated in Fig. 2.14.

Table 2 .

 2

		Marching squares	Q 4 voxel based	grey-level proper-
		algorithm	(Fig.	model (Fig. 2.35)	ties (Fig. 2.36)
	2.33)	2.34)
	50 × 43 pixels 798 100 × 86 pixels 798	3276 16432		1728 5916	4488 17574

1: Comparison of the number of degrees of freedom for each image-based model (hexagonal micro-structure problem).

  Table 2.2: Comparison of the simulation parameters of the three image-based discretizations.

	Number of	Mean vol-	Number of	Number of	LU decom-	Elastic
	elements	ume	of	DOFs (3×	non	zero	position		potential
		elements		the number	values in K	time	(4	energy
				of nodes)			CPUs		1 2 u T Ku
							used)
	FEM, p = 1 502 107 FCM, p = 3 6144, (16×24×16)	(1/2) 3 5 3		478 776 28 815	16 230 798 22 980 735	2 minutes 5.3 minutes 1.72 × 10 10 1.73 × 10 10

can mention for example iterative solvers: multigrid solvers [de Prenter et al. 2020, Jomo et al. 2019] and domain decomposition approaches based on inexact local solvers [Bosy et al. 2020].

Table 4 .

 4 1.

	Plastic strain 0%	0.2%	1%	10%
	Yield stress	230 MPa	295 MPa	340 MPa	425 MPa

Table 4

 4

.1: Elasto-plastic law used for the reference FE simulation.

Table 5 .

 5 1: Global residual indicator for the different loading states and mesh discretizations.

	Element size	Deformation state	State 1	State 2	State 3	State 4
	33 voxels		2.61%	3.62%	6.96%	No convergence
	23 voxels		2.56%	3.36%	6.32%	8.48%
	12 voxels		2.41%	3.18%	5.43%	7.63%

  H 0,1 H 0,2 H 1,0 H 1,1 H 1,2 H 2,0 H 2,1 H 2,2

								T  	H 0,0  	  	∂ux ∂x ∂uy ∂y ∂y + ∂ux	∂uy ∂x	   dxdy	(A.65)
	=	Ω	∂v x ∂x	H 0,0	∂u x ∂x	+ H 0,1	∂u y ∂y	+ H 0,2	∂u x ∂y	+	∂u y ∂x
	+	∂v y ∂y	H 1,0	∂u x ∂x	+ H 1,1	∂u y ∂y	+ H 1,2	∂u x ∂y	+	∂u y ∂x
	+	∂v x ∂y	+	∂v y ∂x		H 2,0	∂u x ∂x	+ H 2,1	∂u y ∂y	+ H 2,2	∂u x ∂y	+	∂u y ∂x	dxdy.

available in http://www.swig.org/.

(a) Reference configuration: image I r . (b) Deformed configuration: image I d .

A special thanks is adressed to Abdallah Bouzid from Institut Clément Ader for manufacturing the samples.

We warmly thank Laurent Crouzeix for his help and his supevision during this experiment.

We thank him for his precious help.

We would like to thank him warmly for his kind help and the fruitful exchanges.

https://anr.fr/Project-ANR-19-CE05-0005.

https://pprime.fr/en/research/mechanical-engineering-and-complex-systems/.

Remerciements

Caroline Chevillote pour ton aide à l'IUT. Merci Léo Lebrat pour m'avoir initié à l'interfaçage Python/C++

Part II Kinematic measurements using image-based mechanical modeling C that links the basis functions of the mesh M 1 to the basis functions of the mesh M 2 (N 1 = CN 2 ). In this case the two coarse and fine displacement vectors are related by the iso-geometric constraint:

(3.60)

The operator C can be efficiently constructed using the kronecker product of univariate knot insertion operators as it was defined in [Borden et al. 2011] and can be even generalized for multi-level local refinement [D' Angella et al. 2018].

3.7.1.3 L 2 projection (Case 3)

The multi-grid initialization step is usually faster with linear basis functions that is why the different refinement steps are performed using linear elements. When the final correlation (which corresponds to the best match) is performed with quadratic or cubic B-splines the displacement exchange between the linear mesh and the high order mesh is performed using an L 2 projection. In this case, we have

(3.61)

The optimal condition of problem (3.61) is given by:

(3.62)

The form (3.62) is discretized with the mesh M 2 . It finally leads to solving the following linear system (3.64)

The mass matrix M is integrated with a Gauss integration rule based on 2⌈(p 2 + 1)/2⌉ points for each univariate direction. The right hand side is integrated with ⌈(p 1 + 1)/2⌉ + ⌈(p 2 + 1)/2⌉ Gauss integration points for each direction. p 1 and p 2 represent respectively the degrees of the B-spline basis functions of M 1 and M 2 .

Point location in tetrahedral finite elements (Case 4)

When the displacement field exchange is performed between two C 0 linear tetrahedral finite element meshes, one need to determine in which elements of M 1 belong the nodes of M 2 in order to use the finite element interpolation. This geometric problem is very common and used in finite elements. In the case of tetrahedral finite element meshes we consider the following strategy:

• Location of the nearest face of the finite-element mesh to the point:

Efficient data structures (very common in collision detection) can be efficiently used to speed up point queries with respect to complex geometric objects represented by faces. In this work, we mainly use a rootine from the CGAL library [The CGAL Project 2021] and PyMesh [PyMesh 2020]. Efficient point queries such as intersections, distance computation, ray shooting

General conclusion and prospects

Synthesis

An efficient characterization of cellular solids must rely on accurate kinematic measurements below the cellular scale. This work attempts to improve the measurement of displacement and strain fields at the sub-cellular scale and therefore allows a non-negligible reduction of the measurement uncertainy. The originality of the present contribution lies in a rigorous link between the regularized optical flow [Horn & Schunck 1981] (a.k.a voxel scale DVC [Réthoré et al. 2009, Leclerc et al. 2011]) and automatic image-based mechanical modeling.

An image-processing strategy consisting in pre-filtering the considered images along with a spline representation is considered. It allows to remove noise effects and oscillations which is beneficial for both image-based mechanical modeling (for getting a regular boundary) and image-correlation (for getting regular image gradient fields).

For the image-based mechanical modeling aspect, two numerical methods are implemented and investigated: the finite element method and the finite cell method. For the finite element method, image-based mesh generation was mainly based on the use of the efficient CGAL library [The CGAL Project 2021]. For the finite cell method, we performed an efficient space-tree assembly strategy based on tensor-product of pre-computed integrated univariate or bivariate differential operators. Furthermore, when considering image-based reconstructed level-sets, non-negligible geometric errors are introduced by both the nature of the level-set and the threshold setting. A space-tree level integration rule based on the image resolution is considered in order to take into account the geometric error coming from the pixellation process.

We show that the finite cell method is largely more accurate than lower order finite element approaches given the same amout of degrees of freedom [START_REF] Rouwane | [END_REF]. This is obtained thanks to the regularity of high order B-spline basis functions. However, the suggested numerical implementation that consists in assembling the complete linear systems and solving them with standard direct solvers (Cholesky) seems to be quite expensive especially from a memory point of vue as the obtained linear systems are smaller but denser. This is the main point that must be ameliorated in order to solve large scale problems.

The developed image-based mechanical models are then used as a regularization for DIC/DVC measurements in cellular architectures (ADDICT approach). A first study based on artificial cellular geometries and reference finite element displacement fields allowed to assess the ability of the AD-DICT method to measure displacements and strains under the cell scale. The L-curve experimental study which consists in finding the best trade-off between data fidelity (grey-level dispersion) and regularity (elastic equilibrium gap) allowed to detect the optimal regularization length which corresponds approximately to the cell length. These results were published in [Rouwane et al. 2022].

The obtained regularization criteria was afterwards used in the three-dimensional generalization of the ADDICT approach. Similar results are obtained for volumetric measurements where complex localized kinematics such as buckling of multiple individual struts are identified.

We show that it is possible to estimate displacements and their derivatives within the strut thickness subdivision points are given by the equations in (A.1), (A.2) and (A.3).

3) shows that the Bézier curve at a parameter ξ is the linear combination of the control points by the coefficients given by the Bernstein polynomials

. For a degree n, the Bernstein polynomials form a basis for the vector space of polynomials of degree at most n and all sum to 1. B-splines are piece wise polynomial functions that generalize Bernstein polynomials. The idea of B-splines is to subdivide the parametric space into parametric elements called knot spans. The subdivision information is given by a non-decreasing set of parametric values called a knot vector and that we will denote ξ ξ ξ = (ξ 0 , ..., ξ m ). The sub-interval [ξ i , ξ i+1 [ is the i-th knot span. A B-spline basis function of degree p can be defined by the Cox-deBoor recursive formula [Piegl & Tiller 1996]

For a degree p, each basis function N i,p is a linear combination of two p -1 degree basis functions.

As an example, we compute the basis functions for the uniform knot vector ξ ξ ξ = (0, 1, 2, 3, 4).

.

with

, where X, Y, Z ∈ R n are the control points coordinates that define the control polygon. Usually, when using knot vectors, the first and last knots are repeated p + 1 times in order to have a boundary interpolation condition (so that we can impose boundary conditions). This type of knot vectors are called open knot vectors. In the opposite, if the first and last knots are not repeated, we call the knot vector a closed knot vector (see the example of suitable for analysis and refer the reader to the reference [Piegl & Tiller 1996] for additial examples and detailed mathematical proofs.

• Positivity:

• Partition of unity:

• Number of non zero basis functions on a knot span:

Given a knot span [ξ j , ξ j+1 [, only the p + 1 basis functions N j-p,p , .., N j,p are non zero.

• Regularity :

If ξ is a knot of multiplicity k (which means that it appears k times in the knot vector), then

A.1. B-splines: Definition and properties 150

• Linear independence :

• Basis :

N i,p form a basis for the vector space of piecewise polynomials of degree p. This vector space is of dimension mp where m + 1 is the number of knots in the knot vector ξ ξ ξ.

In addition to B-spline geometries, NURBS (Non uniform rational B-splines) allow to represent exact conics. The general formula of a NURBS curve is defined by Eq. (A.10). Once the initial B-spline curve is created, it is practically interesting to refine it for both geometry and analysis purposes. The refinement could be done in the parametric space by adding knots without altering the shape of the initial curve. This type of refinement is called knot insertion or h refinement with the iso-surface constraint. Another type of refinement could be obtained by increasing the degree of the B-spline basis functions and it is called order elevation or p refinement. Fig . A.5 illustrates the two refinement mechanisms. Formulas and algorithms for knot insertion and degree elevation are detailed in [Piegl & Tiller 1996].

The transition into B-spline surfaces and solids is done by performing the tensor product of the univariate B-spline basis functions as it follows:

• B-spline surface of degrees (p, q) defined on (ξ ξ ξ, η η η) knot vectors (S : R 2 -→ R 3 ): .11) with P i,j = (X i,j , Y i,j , Z i,j ) ∈ R 3 , where X, Y, Z ∈ R n×m are the control grid coordinates.

adding n x -1 and n y -1 knots) in the x and y directions. The initial control grid for one patch is defined as follows:

In this thesis, it is the only B-spline geometry (in addition to its volume generalization which is a cube) that is used for both the finite cell method and for image transformations. The approach of using this B-spline grid over images was first suggested in [Rueckert et al. 1999] but using the non-interpolatory formalism as defined in chapter 1.

A.2.2 A random surface

Let x and y the vectors representing respectively the uniform subdivision of a plate defined on [x xmin , x max ] with p elements and [y min , y max ] with q elements. We define 1 n the vector of dimension n containing only ones. The B-spline plate is defined by:

For p = 3, we modify the height of the control points by elevating the height coordinates as follows

We obtain the surfaces displayed in

A.2.3 Cuboid

In order to generate a B-spline grid of arbitrary order, we start from one trilinear C 8 element defined by the following tensor product control points X, Y, Z in R 2×2×2 defined by: A.2. Generation of some simple surfaces and volumes 154

A.2.4 Anulus

We define an anulus by its center, its minimal radius denoted r min and its maximal radius denoted r max . Let us consider a center equal to (0, 0). The minimal B-spline discretization for this geometry is defined by Eq. (A.19)-(A.20)-(A.21).

ξ ξ ξ = (0, 0, 0, 1/4, 1/4, 1/2, 1/2, 3/4, 3/4, 1, 1, 1), η η η = (0, 0, 1, 1), (A.19) X = 0 r max r max r max 0 -r max -r max -r max 0 0 r min r min r min 0 -r min -r min -r min 0 T , (A.20) Y = -r max -r max 0 r max r max r max 0 -r max -r max -r min -r min 0 r min r min r min 0 -r min -r min T , (A.21)

The control grids before and after knot insertion are shown in

A.2.5 Disk

A disk can be defined as an anulus with a zero radius value (see Fig. (K e )) I,J = (C e ) i,j

(K e ) I+nbf,J = (C e ) i+me,j

(K e ) I,J+nbf = (C e ) i,j+me

With the remaining terms of the matrix equal to 0, (A.55) where m is the total number of control points (m = ndof /dimension). We can finally assemble the global stiffness matrix by summing the elementary matrices:

The right hand side f is assembled similarly by summing the contributions of elementary forces defined on the boundary Γ n . We note that the methodology presented previously is well known and efficiently established and is exactly the same as the classical FE assembly rootines (generaly written in low level languages such as C/C++ or Fortran). The advantage of finite element assembly routines is that they usually define a reference element on which the basis functions are defined and evaluated once for all only one time at the reference element's integration points which is not the case for IGA. Multiple contributions have emerged in order to come closer to the finite element assembly using either Bézier extraction [Borden et al. 2011] or Lagrange extraction [Schillinger et al. 2016, Tirvaudey et al. 2019] or by considering user-defined elements in industrial softwares [Lai et al. 2017].

A second method for assembling the stiffness matrix is to construct global operators for the evaluated basis functions and their derivatives on all the integration points [Passieux 2018]. We set ng the total number of two dimensional integration points. We define W g the diagonal matrix of size ng × ng defined by:

mes(e(g i )) 4 where w(g i ) is the Gauss weight associated to the Gauss point g i and mes(e(g i )) is the measure of the parametric element e containing the integration point g i . We define N ∈ R ng×m a sparse matrix which row i contains the evaluation of the bi-variate basis functions at the Gauss integration point g i (see Fig . A.15 for a visualization of the matrix).

We define similarly the matrices ∂ x N ∈ R ng×nbf and ∂ y N ∈ R ng×nbf of derivatives of the basis functions at the Gauss points. Let N x , N y , ∂ x N x , ∂ x N y , ∂ y N x , ∂ y N y ∈ R ng×ndof the matrices defined by:

( A.57) This trick of stacking a zero matrix is used here because we are treating a vector field in our PDE.

If we only had a scalar field, we would have only needed the matrices N, ∂ x N and ∂ y N. Now that the differential operators and the integration weight operators are defined, we can write integration operations using matrix/vector product. Here are some computation of some integrals: .59) where 1 ∈ R ng is a vector containing only ones. As an other example, we have: Using the vectorized discretization method, the weak form is written:

• Boundary elimination:

Let b be the index sets of the boundary nodes for which a displacement is known (i.e. a dirichlet boundary condition is imposed) and let i its complementary set. The equilibrium stiffness system can be written as follows

where f i is the vector of internal and external forces at the nodes of index i. Imposing Dirichlet boundary conditions consists simply in solving

A reduced system is solved at the internal and Neumann nodes, then the final degrees of freedom vector u is rearranged in its original order. The advantage of this approach is that we can keep using the Cholesky solver and benefit from its numerical efficiency.

A.5 Examples of Isogeometric analysis applications in linear elasticity

A.5.2 Curved thick beam subject to end shear deformation

This second example treats the problem of a curved beam subjected to end shear. The geometry and the boundary conditions are illustrated in Fig. A.22. In order to represent exactly the circular beam, a first order parametrization is sufficient for the radial direction and a minimal second order parameterization is needed for the tangential direction. A constant displacement u x = u 0 is applied along the line y = 0. A symmetric boundary condition u y = 0 is applied for x = 0 in addition to fixing the point (x = 0, y = r min ) by setting u x = u y = 0. The geometrical parameters of the problem are defined by: r min = 5, r max = 10, u 0 = 0.01, E = 10 5 , ν = 0.25.