G Legrain

M Chevreuil

N Takano Gem

N Takano

Prediction of apparent properties with uncertain material parameters using high-order fictitious domain methods and PGD model reduction

Keywords: fictitious domain method, homogenization, high-order, stochastic, proper generalized decomposition

This contribution presents a numerical strategy to evaluate the effective properties of imagebased microstructures in the case of random material properties. The method relies on three points: (i) a high-order fictitious domain method; (ii) an accurate spectral stochastic model and (iii) an efficient model reduction method based on the Proper Generalized Decomposition in order to decrease the computational cost introduced by the stochastic model. A feedback procedure is proposed for an automatic estimation of the random effective properties with a given confidence. Numerical verifications highlight the convergence properties of the method for both deterministic and stochastic models. The method is finally applied to a real 3D bone microstructure where the empirical probability density function of the effective behaviour could be obtained.

Introduction

During the past few years, image-based analysis has been the object of a large interest for a wide range of applications. In particular, a large number of studies have been conducted in the domain of biomechanics since the 90's [START_REF] Keyak | Automated threedimensional finite element modelling of bone: a new method[END_REF][START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: Use of different densitometric measures[END_REF] for building patient specific models of bones for osteoporosis treatments. The finite element predictions are used to decide whether surgical interventions are needed, but also to understand the relationship between bone morphology and load transfer. Besides osteoporosis applications, patient specific finite element models have also been considered for the optimization and the design of orthopaedic implants (see e.g. [START_REF] Bongini | Hippin: A semiautomatic computer program for selecting hip prosthesis femoral components[END_REF][START_REF] Dick | Stress tensor field visualization for implant planning in orthopedics[END_REF]). These techniques have also been widely applied in material science [START_REF] Terada | Digital image-based model-ing applied to the homogenization analysis of composite materials[END_REF][START_REF] Coffman | OOF3D: An image-based finite element solver for materials science[END_REF] to understand the relationship between the morphology and the mechanical properties of the materials. Two main families of methods have been designed for building finite element meshes from geometrical data (CT scans, microscopy). Images being represented as a multidimensional array of pixels (voxels in 3D), voxel-based finite element approaches [START_REF] Keyak | Automated threedimensional finite element modelling of bone: a new method[END_REF][START_REF] Van Rietbergen | A new method to determine trabecular bone elastic properties and loading using micromechanical finite-element models[END_REF][START_REF] Hollister | Homogenization theory and digital imaging: a basis for studying the mechanics and design principles of bone tissue[END_REF][START_REF] Takano | Microstructure-based stress analysis and evaluation for porous ceramics by homogenization method with digital image-based modeling[END_REF][START_REF] Salleh Basaruddin | Stochastic multiscale prediction on the apparent elastic moduli of trabecular bone considering uncertainties of biological apatite (BAp) crystallite orientation and image-based modelling[END_REF] can directly translate geometrical acquisitions into numerical models by simply associating a quadrangular (resp. hexahedral) element to every pixel (resp. voxel). In this case, material properties can be assigned based on a segmentation, or directly based on the intensity level of the image (provided that a relationship between image intensity and material properties is given). This approach ie almost completely automatic, but leads to large numerical models (directly proportional to the image resolution). In addition, the stair steps approximation of the geometry is known to decrease stress accuracy near boundaries in addition to the convergence of the model [START_REF] Lian | Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxelbased approaches[END_REF]. The second family of methods is based on an image segmentation step, followed by surface triangulation and finally volume meshing [START_REF] Frey | Fully automatic mesh generation for 3-D domains based upon voxel sets[END_REF][START_REF] Ulrich | Finite element analysis of trabecular bone structure: a comparison of image-based meshing techniques[END_REF][START_REF] Maire | X-ray tomography applied to the characterization of cellular materials related finite element modeling problems[END_REF]. This method allows the use of non-homogeneous (user-defined) mesh density, and thus improve the efficiency of the numerical model. However, mesh generation of complex micro-structures is still a challenging task, especially if distorted elements have to be avoided. To answer the issues raised by these two approaches, alternative methods have emerged recently with the objective to simplify the interaction between geometry and numerical computations. One can cite for example the Finite Cell method from Rank, Düster and co-workers [START_REF] Düster | The finite cell method for three-dimensional problems of solid mechanics[END_REF][START_REF] Rank | Shell Finite Cell Method: A high order fictitious domain approach for thin-walled structures[END_REF][START_REF] Parvizian | Finite cell method -h and p extension for embedded domain problems in solid mechanics[END_REF][START_REF] Schillinger | Small and large deformation analysis with the p-and B-spline versions of the Finite Cell Method[END_REF], and the eXtended Finite Element Method (X-FEM in the following) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Sukumar | Extended Finite Element Method for three-dimensional crack modelling[END_REF][START_REF] Stolarska | Modelling Crack Growth by Level Sets and the Extended Finite Element Method[END_REF][START_REF] Sukumar | Modeling Holes and Inclusions by Level Sets in the Extended Finite Element Method[END_REF][START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF][START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model[END_REF][START_REF] Gravouil | Non-planar 3D crack growth by the extended finite element and level sets. Part II: level set update[END_REF][START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF][START_REF] Legrain | Numerical simulation of large strain fracture problems using X-FEM[END_REF][START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF][START_REF] Legrain | Stability of incompressible formulations enriched with X-FEM[END_REF][START_REF] Legrain | Treatment of nearly-singular problems with the X-FEM[END_REF]. The Finite Cell method is a high-order fictitious domain method which has proved to be very efficient thanks to the approximation properties of p-fem [START_REF] Babuška | On the rates of convergence of the finite element method[END_REF][START_REF] Szabó | The p-and h-p versions of the finite element method in solid mechanics[END_REF][START_REF] Szabó | Finite Element Analysis[END_REF]. The fictitious domain nature of the method makes possible to use meshes that do not conform to the geometry, avoiding the use of cumbersome meshing techniques. Instead, the effort is put on the integration of the weak formulation, as the geometrical error in the solution has to be mastered. The method has been applied in the context of design through analysis [START_REF] Ernst Rank | Geometric modeling, isogeometric analysis and the finite cell method[END_REF], but also to image-based analysis in both material science [START_REF] Düster | Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method[END_REF][START_REF] Heinze | Numerical homogenization of hybrid metal foams using the finite cell method[END_REF] and biomechanics [START_REF] Ruess | The finite cell method for bone simulations: Verification and validation[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF]. The use of the X-FEM in the context of CAD and fictitious domain traces back to the work of Belytschko [START_REF] Belytschko | Structured extended finite element methods for solids defined by implicit surfaces[END_REF]. The use of the X-FEM in the context of image-based analysis has been further developed in [START_REF] Lian | Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxelbased approaches[END_REF][START_REF] Legrain | An X-FEM and level set computational approach for image-based modelling: Application to homogenization[END_REF] where a complete strategy was proposed to build accurate numerical models from 2D or 3D images. The use of an octree data structure allowed to decrease the computational cost, but also required a special treatment for enriched hanging nodes [START_REF] Legrain | On the use of the extended finite element method with quadtree/octree meshes[END_REF]. The method has been extended to higher order approximations in [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF], and proved to be efficient for image-based analysis. Despite being very accurate, these methods assume the material properties to be deterministic. This is not the case in practice, as a large scattering of the properties was observed in the numerous experimental data that are taken as reference for the material properties [START_REF] Salleh Basaruddin | Stochastic multiscale prediction on the apparent elastic moduli of trabecular bone considering uncertainties of biological apatite (BAp) crystallite orientation and image-based modelling[END_REF].These uncertainties can have multiple sources: (i) geometrical uncertainties due to the geometrical acquisition process, the associated image treatments and the geometrical discretization; (ii) material uncertainties due to inter-individual variation in the material properties of the tissues or even intra-individual variations depending on the biological site where the bone is located and (iii) experimental uncertainties due to the processing of the samples, mis-alignments in the experimental devices, and sensors accuracy. In this contribution, we will restrict ourselves to material uncertainties alone. This source of uncertainty is of prime importance in the context of patient-specific models, as they allow to model both intra and inter-individual material variations. For patient-specific models, geometrical uncertainties are mainly linked to the resolution of the image which is considered as sufficiently large (see the example in section 4). Usually, Monte-Carlo strategies are considered to model the propagation of uncertainties in homogenized models. The X-FEM is appealing in this context, and was used in [START_REF] Hiriyur | Uncertainty quantification in homogenization of heterogeneous microstructures modeled by XFEM[END_REF] and [START_REF] Savvas | Homogenization of random heterogeneous media with inclusions of arbitrary shape modeled by XFEM[END_REF] to study the statistics of the effective properties against the distribution and shape of material inclusions. In the context of FEM, [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF] proposed statistical methods in order to estimate the size of the representative volume element of random composites and the associated homogenized properties. Mathematical bonds on the statistics of the effective properties were also obtained in [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF][START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF] when the RVE is sub-sampled and in [START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF] in the case of random micro-structures. So far, all these contribution assumed fixed material properties for the different phases. In Basaruddin et al. [START_REF] Salleh Basaruddin | Stochastic multiscale prediction on the apparent elastic moduli of trabecular bone considering uncertainties of biological apatite (BAp) crystallite orientation and image-based modelling[END_REF], both geometrical and material randomness were taken into account by means of perturbation techniques. Wille et al [START_REF] Wille | Prediction of the mechanical response of the femur with uncertain elastic properties[END_REF] also studied the impact of the variability of the material properties on the mechanical response of femurs. This approach was relying on a simple product of the deterministic response by an identified stochastic function.

The aim of this manuscript is to take into account state of the art stochastic methods for an accurate representation of the stochastic output of the model. The stochastic finite element method will be considered in the following [START_REF] Roger | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Matthies | Stochastic finite elements: Computational approaches to stochastic partial differential equations[END_REF]. However, this method is known to require high computational resources. This is why model reduction techniques will be considered so that large stochastic models can be solved accurately. The accuracy of the model will be ensured, not only at the deterministic level, like in [START_REF] Miranda Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] but for the whole stochastic model through feedback strategies [START_REF] Szabó | Finite Element Analysis[END_REF] that proved to be efficient in the p-fem context. As stated above, randomness is considered only for the material properties (to model inter-individual variations). In the case where the resolution of the input images is not sufficient, one could take into account the influence of the image treatment following the approach proposed in [START_REF] Salleh Basaruddin | Stochastic multiscale prediction on the apparent elastic moduli of trabecular bone considering uncertainties of biological apatite (BAp) crystallite orientation and image-based modelling[END_REF] which would lead to a random levelset.

The manuscript is organized as follows: first, the numerical strategy is introduced, both at the deterministic and stochastic levels. Emphasize is also put on the PGD (Proper Generalized Decomposition) model reduction technique [START_REF] Amine Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF] in the context of stochastic analysis. In a second part, the numerical strategy is validated on a 2D model problem, before being applied to a real 3D cortical bone micro-structure.

Numerical strategy

The numerical strategy which is proposed in this contribution is based on three ingredients: (i) an efficient and accurate numerical scheme for computing the overall properties of image-based structures, (ii) an accurate stochastic model to take into account the uncertainties in the model and (iii) an efficient model reduction method to decrease the computational cost of the stochastic model. These aspects are detailed in the following sections.

High order fictitious domain method

One of the main difficulties for simulating image-based micro-structures concerns the treatment of the amount of geometrical informations produced during the acquisition. Indeed, 1024 3 voxels images are now common: such high-resolution images require large computer resources for the simulations unless tedious geometrical and meshing treatments are used. In this contribution, an alternative path is followed: geometrically non-conforming approaches in the spirit of fictitious domain methods are considered in order to alleviate the meshing burden. The method is illustrated in figure 1 a: The approximation is defined on a regular background grid of simplices or hypercubes which does not conform to the geometry. Elements that are fully in the domaine are integrated with classical quadrature rules, whereas ad hoc strategies have been developed for elements lying partially in the domain. The objective is to integrate accurately the weak formulation on these elements. The Finite Cell method [START_REF] Düster | Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method[END_REF][START_REF] Parvizian | Finite cell method -h and p extension for embedded domain problems in solid mechanics[END_REF][START_REF] Rank | Shell Finite Cell Method: A high order fictitious domain approach for thin-walled structures[END_REF][START_REF] Schillinger | An unfitted hp-adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry[END_REF][START_REF] Schillinger | Small and large deformation analysis with the p-and B-spline versions of the Finite Cell Method[END_REF] initially used space-tree algorithms and integration points classification (see figure 1 b), whereas so-called sub-grid level-sets was advocated in [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF] (see figure 1 c) where the geometry is represented using the iso-zero of a level-set function [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF]. With respect to space-tree algorithms, this method employs a local tesselation near the interface. More recently, the use of suitable mappings was proposed in order to decrease the integration cost when the interfaces are represented by means of parametric surfaces [START_REF] Kwok | Higher-order XFEM for curved strong and weak discontinuities[END_REF][START_REF] Legrain | A NURBS enhanced extended finite element approach for unfitted CAD analysis[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary methods[END_REF][START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF][START_REF] Marco | Exact 3D boundary representation in finite element analysis based on Cartesian grids independent of the geometry[END_REF] (see figure 1 d). As we are dealing with discrete geometries, we will focus on sub-grid level-sets in the following. The computational strategy proposed in [START_REF] Legrain | An X-FEM and level set computational approach for image-based modelling: Application to homogenization[END_REF][START_REF] Lian | Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxelbased approaches[END_REF] is used to convert the images into numerical models using level-sets.

Thanks to this approach, it is possible to treat very complex geometries with an approximation defined on a regular grid. Accurate results with such coarse meshes are obtained thanks to the use of high-order approximations (typically up to p = 10 in 2D, and p = 8 in 3D). This method has been proved to be very effective in previous contributions. This stems from the fact that p methods can attain exponential convergence whereas classical h finite elements are limited to algebraic convergence [START_REF] Szabó | Finite Element Analysis[END_REF].

Numerical homogenization

We are interested in the computation of the overall properties of material micro-structures. The output of the procedure is the effective elasticity tensor C which relates the macroscopic stress and strain tensors, Σ and E:

Σ = C E (1) a) b) c) d)
Figure 1: Unfitted approximation of a body, and quadrature points for an element crossed by the interface: b) Finite Cell quad-tree algorithm, c) sub-grid levelset and d) transfinite mapping.

where Σ and E are defined by the following spatial averages:

Σ = σ , (2)
E = ε , (3)
< • > = 1 V Ω • dΩ (4)
where tensors σ and ε are the micro stress and strain tensors, Ω is the spatial domain, and V its volume. The homogenization procedure is based here on kinematic uniform boundary conditions (KUBC), but could also be based on stress uniform boundary conditions (SUBC) or periodic boundary conditions (PBC) [START_REF] Kanit | Determination of the size of the representative volume element for random composites: statistical and numerical approach[END_REF][START_REF] Kanit | Apparent and effective physical properties of heterogeneous materials: Representativity of samples of two materials from food industry[END_REF][START_REF] Ostoja-Starzewski | Material spatial randomness: From statistical to representative volume element[END_REF]:

u = E • x on ∂Ω ( 5
)
where u is the displacement. Using these boundary conditions, equation ( 3) is automatically verified, whereas the components of C are obtained by computing eq.( 2) for three macroscopic displacement modes associated to macroscopic strain fields E 1 , E 2 and E 3 (2D macroscopic tractions and shear):

E 1 = 1 0 0 0 , (6)
E 2 = 0 0 0 1 , (7)
E 3 = 0 0.5 0.5 0 . ( 8
)
If a so called energetic extraction is considered, then the components of C can also be obtained by means of the Hill-Mandel condition:

C ij = E i .C.E j = a(u i , u j ) (9)
where a(u, v) = Ω ε(u) : K : ε(v) dΩ is the bilinear form of the weak formulation of the mechanical problem with K the fourth order Hooke tensor, and u i is the microdisplacement associated to E i in [START_REF] Terada | Digital image-based model-ing applied to the homogenization analysis of composite materials[END_REF]. If u h,p i is the finite element solution of the weak formulation of the unit-cell problem, then the associated finite element effective tensor verifies:

C h,p ij = E i .C h,p .E j = a(u h,p i , u h,p j ) (10)
Following [START_REF] Miranda Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods[END_REF], we can write that:

C h,p ij -C ij = a(u h,p i , u h,p j ) -a(u i , u j ) (11) = a(u h,p i -u i , u h,p j -u j ) (12)
Using Schwarz' inequality, the expression can be bounded:

|C h,p ij -C ij | ≤ a(u h,p i -u i , u h,p i -u i ) 1/2 a(u h,p j -u j , u h,p j -u j ) 1/2 (13) Noting that a(u h,p i -u i , u h,p i -u i ) 1/2
is the energy norm of the error of micro-problem i, we can use classical finite elements error estimates:

e i E = a(u h,p i -u i , u h,p i -u i ) 1/2 ≤ Ch p CN -p/dim for h finite elements, ( 14
)
≤ k/ exp(γN θ ) for p finite elements, ( 15
)
where dim is the dimension of the problem, N the number of degrees of freedom and C, C, k, β and θ are problem dependent constants. It can be shown that θ > 1/2 in 2D (resp. 1/3 in 3D) if the solution is regular. Otherwise, the converge rate remains algebraic, see [START_REF] Szabó | Finite Element Analysis[END_REF] for the derivation of these estimates. This means that the convergence of the components of effective tensor behaves like the error of the micro-solution:

|C h,p ij -C ij | ≤ Ch 2p CN -2p/dim for h finite elements, (16)
≤ k/ exp(2γN θ ) for p finite elements.

Note however that C h,p ij may converge to C ij from above or below depending on the component. Thus, the absolute value in ( 16) can be dropped for the diagonal components, but not for the off-diagonal ones. Nevertheless it is still possible to apply the extrapolation estimate proposed in [START_REF] Szabó | Finite Element Analysis[END_REF] to estimate the error on these quantities of interest. Indeed, [START_REF] Szabó | Finite Element Analysis[END_REF] shows that C ij can be estimated from three computations with degrees p, p -1, p -2 by solving the following nonlinear equation:

C h,p ij -C ij C h,p-1 ij -C ij =   C h,p-1 ij -C ij C h,p-2 ij -C ij   Q (18) with Q = log N p-1 Np log N p-2 N p-1
. Although less powerful, this approach is easier to implement than goal-oriented error estimators which require 6 (21 in 3D) additional adjoint resolutions to estimate the error. The performances of this adaptive method will be studied in the following sections.

Stochastic Finite Elements

We are interested in the prediction of the impact of uncertainties on the overall behaviour of materials. This impact is traditionally obtained using sampling techniques, such as Monte Carlo simulations. However, such a strategy is intractable because of the large number of large simulations it requires (O(1/ N samples ) convergence). Perturbationbased methods [START_REF] Kleiber | The stochastic finite element method (basic perturbation technique and computer implementation)[END_REF][START_REF] Salleh Basaruddin | Stochastic multiscale prediction on the apparent elastic moduli of trabecular bone considering uncertainties of biological apatite (BAp) crystallite orientation and image-based modelling[END_REF] can be used with good results. However, the accuracy of such methods is satisfactory only for a small range of variation of the input parameters. If not, one has to rely on improved schemes like [START_REF] Falsone | A new approach for the stochastic analysis of finite element modelled structures with uncertain parameters[END_REF]. We will focus here on a functional point of view for stochastic analysis: random quantities are seen as functions of random parameters that represent basic uncertainties. This vision has led to the development of spectral stochastic methods [START_REF] Roger | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Xiu | Fast Numerical Methods for Stochastic Computations: A Review[END_REF][START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF][START_REF] Matthies | Stochastic finite elements: Computational approaches to stochastic partial differential equations[END_REF]. The main advantage of these methods is the explicit representation of the random response of the model as a function of the basic random parameters. They also allow for fast post-processing of the solution (envelope of the solution, sensitivity with respect to the input parameters,...).

In the following, uncertainties are represented by means of a finite number of random variables ξ = {ξ 1 , ξ 2 , • • • ξ r } which are random parameters of the mechanical problem (Young's modulus and Poisson's ratio in our case). Let (Θ, B, P ξ ) be the associated probability space, where Θ ∈ R r is the range of ξ, B the associated σ-algebra and P ξ the probability measure of ξ. The spatial weak formulation of the mechanical problem is:

Find u(x, ξ) ∈ V Ξ such that a(u, v; ξ) = b(v; ξ) ∀v ∈ V Ξ 0 (19)
where u(x, ξ) is the unknown random displacement field, Ω is the spatial domain of interest, and a and b represent respectively the bilinear and linear forms of the problem:

a(u, v; ξ) = Ω ε(u) : K(ξ) : ε(v) dΩ, (20) b
(v; ξ) = Ω v.b(ξ) dΩ + Γ N v.T(ξ) dΓ. (21)
In these expressions, K(ξ) is the random fourth order Hooke tensor, b(ξ) are the random body forces and T(ξ) are the random prescribed tractions on Γ N ∈ ∂Ω. Under suitable regularity assumptions on the random parameters, the solution of ( 19) is searched in L 2 (Θ, B, P ξ ; V) the space of real value square integrable functions defined on the probability space (Θ, B, P ξ ) with values in V = H 1 (Ω). Note that this space can be identified with the tensor product space V ⊗ S, with S = L 2 (Θ, B, P ξ ) the space of real value second order random variables. A finite element approximation of the solution can then be obtained by introducing a finite dimensional space V N ⊂ V such that:

u(x, ξ) ≈ u h,p (x, ξ) = N α=1 u α (ξ)N α (x) (22)
where {N α } N α=1 is the basis of functions of the approximation space V N and u α are the N random degrees of freedom (dofs). Spectral stochastic methods consist in looking for an approximation of the dofs u α (ξ) under the expansion:

u α (ξ) ≈ P k=1 u α k H k (ξ) (23)
where {H k (ξ)} P k=1 is a multidimensional basis of r parameters (ξ 1 , • • • ξ r ) of an approximation space S P ⊂ S (e.g. the polynomial chaos [START_REF] Roger | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Xiu | Fast Numerical Methods for Stochastic Computations: A Review[END_REF], piecewise polynoms, wavelets [START_REF] Manas | Solution of stochastic partial differential equations using Galerkin finite element techniques[END_REF][START_REF] Le Matîre | Uncertainty propagation using Wiener-Haar expansions[END_REF][START_REF] Wan | An adaptive multi-element generalized polynomial chaos method for stochastic differential equations[END_REF]). The coefficients of the expansion can be determined using the weak stochastic formulation of ( 19): Find u h,p (x, ξ) ∈ V N ⊗ S P such that:

A(u h,p , v) = B(v) ∀v ∈ V N ⊗ S P (24)
where:

A(u, v) =Exp (a(u(ξ), v(ξ); ξ)) , (25)
B(v) =Exp (b(v(ξ); ξ)) , (26)
where Exp(•) = Θ •dP ξ is the mathematical expectation. The corresponding linear system contains N × P equations, where N = dim(V N ) and P = dim(S P ) may both be large. This makes stochastic finite elements a computationally intensive method. For large deterministic models like those classically encountered in image-based analysis, such an approach would require large computer resources. This difficulty is also known as the "curse of the dimensionality", which motivates the use of a model reduction method for solving the stochastic equations.

PGD spectral model reduction

The Proper Generalized Decomposition [START_REF] Amine Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Chinesta | Recent Advances and New Challenges in the Use of the Proper Generalized Decomposition for Solving Multidimensional Models[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF] is presented in this section for the resolution of problem (24)2 . This method was initially called Generalized Spectral Decomposition (GSD) [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: Invariant subspace problem and dedicated algorithms[END_REF] in the stochastic framework. The idea is to exploit the tensor structure of the solution, and look for an optimal separated representation of u under the form:

u m (x, ξ) = m i=1 w i (x)λ i (ξ) (27)
where w i (x) ∈ V N are deterministic modes, λ i ∈ S P are stochastic functions and m is the rank which is supposed to be small. This separated representation is constructed without knowing the solution of [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model[END_REF]. More recently, it was extended to the full tensorisation of S P for independent random variables (ξ

1 , • • • , ξ r ) with Θ = Θ 1 × • • • × Θ r and dP ξ = dP ξ 1 • • • dP ξr : S P = S P 1 ⊗ S P 2 • • • ⊗ S Pr and: u m (x, ξ) = m i=1 w i (x)φ 1 i (ξ 1 ) • • • φ r i (ξ r ). ( 28
)
This full tensorisation will make possible the use of very fine stochastic discretization for each random variable ξ j at low cost. This decomposition can be built a priori by recasting problem [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model[END_REF] as the minimization of the functional

J(v) = 1 2 A(v, v) -B(v): u(x, ξ) = arg min v∈V N ⊗S P 1 ⊗•••⊗S Pr J(v). ( 29
)
The decomposition is built progressively by successive addition of rank-one elements w(x)φ 1 (ξ 1 ) • • • φ r (ξ r ) to a previously computed decomposition u m-1 as:

u m = u m-1 + w(x)φ 1 (ξ 1 ) • • • φ r (ξ r ) (30)
The additional term is obtained by finding the best additional term defined as:

w(x)φ 1 (ξ 1 ) • • • φ r (ξ r ) = arg min y(x)χ 1 (ξ 1 )•••χ r (ξr)∈ V N ⊗S P 1 ⊗•••⊗S Pr J(u m-1 + y(x)χ 1 (ξ 1 ) • • • χ r (ξ r )) (31)
This nonlinear minimization problem is solved using a simple fixed-point algorithm along each components of w(x)φ 1 (ξ 1 ) • • • φ r (ξ r ). This fixed point algorithm requires the resolution of r 1D stochastic problems and one 3D deterministic problem at each iteration. Usually, 4 to 5 iterations are sufficient to reach stagnation, which makes this algorithm attractive. In this paper, r = 2, as the two random variables parametrize respectively Young's modulus and Poisson's ratio. The fixed point algorithm consists in solving alternatively three problems:

1. Finding w(x) with φ 1 (ξ 1 ) and φ 2 (ξ 2 ) considered as known: this problem has size dim(V N ).

2. Finding φ 1 (ξ 1 ) with w(x) and φ 2 (ξ 2 ) considered as known: this problem has size dim(S P 1 ).

3. Finding φ 2 (ξ 2 ) with w(x) and φ 1 (ξ 1 ) considered as known: this problem has size dim(S P 2 ).

Thus, a collection of moderate size problems has to be solved (4 × m problems of size dim(V N ), 4 × m problems of size P 1 = dim(S P 1 ) and 4 × m problems of size P 2 = dim(S P 2 )), rather than one problem [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets-Part I: Mechanical model[END_REF] of size dim(V N ⊗ S P ) = N × P . The derivation of these problems is presented in appendix A. Note however that the efficiency of the method relies on the low rank separated representations of the operator A and righthand-side B of the weak formulation, these separated representations are precised in the appendix. It was shown that such an algorithm was very robust and converging very fast when efficient update strategies are considered [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: Invariant subspace problem and dedicated algorithms[END_REF]. In order to improve the efficiency of the method, an adaptive stochastic model could be defined at this level. It is not considered here as the model reduction method makes the cost of the stochastic dimensions negligible compared to the deterministic problem (set of 1D problems).

2D numerical verifications

First, 2D numerical verifications are conducted in order to assess the accuracy of the method and highlight its features for both stochastic and deterministic problems. Consider the porous unit-cell of size 1 mm × 1 mm presented in figure 2. It is composed of a linear elastic material whose Young's modulus is 1 Mpa and Poisson's ratio 0.3. The geometry of the hole's interface is represented by the iso-zero of a linear level-set function interpolated on an adapted mesh whose level of recursion is 12 (i.e. elements of size 2.4 10 -4 mm, see figure 2). Note that the geometrical accuracy is fixed for all the following computations. The approximation of the mechanical fields is defined on regular meshes composed of n × n linear elements along each side for h computations (n = 2, 4, • • • , 256), or alternatively 4 × 4 elements along each side for p computations.

Let introduce the relative error on the 2-norm:

(C h,p , C) = C h,p 2 -C 2 C 2 ( 32
)
where

•

h convergence

h-convergence is first considered in order to study the evolution of the different components of the homogenized tensor C. As stated above, the geometrical accuracy is fixed during the convergence, which contrasts with usual low-order finite elements where geometrical accuracy and mechanical accuracy are linked together. The results of the study are presented in figure 3(a) for the reference error on the tensor (C [START_REF] Schillinger | Small and large deformation analysis with the p-and B-spline versions of the Finite Cell Method[END_REF]. Figure 3(a) shows that all the components of the homogenized tensor converge at the expected algebraic rate (O(N ), see section 2.1). The estimated error (C h,p ij , C est ij ) also converges at the expected rate, but the effectivity of the estimator depends on the component of interest (see figures 3(b) and 4). It is very good for the diagonal components of C (and also C 12 ), even for a moderate number of dofs. On the contrary, the estimated error exhibits large oscillations for off-diagonal components. This can be explained by the behaviour of the reference error in figure 3(a): a large number of dofs is required in order to enter the asymptotic range of the convergence. These oscillations prevent the estimator to give reliable results until they stabilize for at least three consecutive computations. This behaviour can be explained as the Galerkin method gives the solution with the lowest error in the energy norm, which is the lowest error in the diagonal terms of C: Off-diagonal terms are not of primary interest. It is worth noting that in this example the components which exhibit these oscillations have small values (400 times smaller than the other ones). Still, they are not zero (because of the geometry of the problem of interest) and even the absolute error exhibits such oscillations. If rather than working component per component we focus on (C h,p , C) , the error on the 2-norm of C, then the influence of these components is discarded and the estimator gives very reliable results, as depicted in figures 3 and 4. Relative error (exact)

C11 C12 C13 C22 C23 C33 C 2

p convergence

The problem is solved by means of p refinement. The computational mesh is highlighted in blue in figure 2), and the spatial polynomial order p ranges from 1 to 7. The convergence of the components of C are shown in figure 5(a) and (b) for the reference and Effectivity of the error estimator

C11 C12 C13 C22 C23 C33 C 2
(C h,p ij ,C est ij ) (C h,p ij ,C ref ij )
and

(C h,p ,C est ) (C h,p ,C ref ) (h convergence).
estimated errors. As expected, the convergence is exponential for all the components. This exponential convergence is smoother for the diagonal coefficients, and is oscillatory for C 13 and C 23 . For the same reasons as in the last section, the estimated error is not able to give reliable results for these two components. Still, working with the 2-norm of C allows to obtain an accurate estimation of the error of this quantity (see figures 5 and 6). The worse performance of the estimator for p-convergence with respect to h convergence which is highlighted in figure 6 may also be related to the small number of dofs for this computation (1400 dofs compared to 100 000 dofs for h convergence). If the range of figure 4 is restricted to the same range of dofs for both h and p approximations, then the behaviour is similar but a little better for h convergence. Besides, the error level for p-convergence is three orders of magnitude smaller than h computations so that the effectivity becomes very sensitive for p computations, even for this low number of dofs. These numerical experiments also allow us to draw some partial conclusions:

• All the components of the homogenized tensor converge exponentially to their limit value.

• For low number of dofs, the convergence is monotonous only for the diagonal terms of C. This is why the extrapolation error estimator is not robust for the other components. More robust error estimators could be considered, see [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF] for a review, or [START_REF] Miranda Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] for homogenization applications.

• If the norm of the effective tensor is of interest, then the estimator was shown to be robust and can be used as a convergence criterion for the feedback procedure (note that there is no guaranty that this observation is true for other materials as it relies on the influence of the off-diagonal terms on the norm). In the following, we will consider the most unfavourable case as we will assume that we are interested in the value of all the components of C individually.

• As the estimator is not reliable in this case, a simple stagnation criterion will be used in the subsequent numerical experiments, as it is known to be a robust indicator (p-p-1 error indicator). The convergence of C will be monitored through the convergence of its individual components (i.e. C is considered as converged if the stagnation criterion is met for all of them). Two stagnation criteria are considered: absolute and relative stagnations. Absolute stagnation is met if the variation of C ij between the two last computations is small with respect to the max-norm of C : Relative error (estimated)

C h,p ij -C h,p-1 ij C h,p-1 max < ε a . Relative stagnation is met if (C h,p ij , C h,p-1 ij ) = C h,p ij -C h,p-1 ij C h,p-1 ij < ε r .
C11 C12 C13 C22 C23 C33 C 2 (

Deterministic Feedback

Based on the conclusions of the last section, the adaptive determination of the homogenized tensor C is now considered. A non-conforming p extension is used to make the overall properties converge. The convergence is monitored trough the stagnation of the components of C. In this numerical example, the computational mesh contains 4 × 4 elements per side and the approximation can range from p = 1 to p = 7. The geometrical mesh is unchanged with respect to the previous section. Finally, the material properties are still considered as deterministic. The process is considered as converged if the relative stagnation of the components of C is lower than ε r = 10 -2 or 5 10 -2 , and

(C h,p ij ,C est ij ) (C h,p ij ,C ref ij )
and

(C h,p ,C est ) (C h,p ,C ref ) (p convergence).
the absolute one lower than ε a = 10 -5 . The evolution of the components of C during the procedure are depicted in figure 7 for these two different stagnation thresholds. It can be seen that the converged quantities can be obtained respectively for p = 5 and p = 3.

Stochastic convergence

The case of stochastic material properties is now considered. We first focus on the convergence of the stochastic homogenized quantities during p refinement. The same porous problem is still considered, but now Young's modulus and Poisson's ratio are random. E ranges from 1/2 to 1.5, and ν from 0.15 to 0.45 both following a uniform probability distribution. Note that more complex probability distributions could be considered without any loss of generality, except unbounded distributions (like the Gaussian distribution). Indeed, negative values (for E and ν) and infinite (for E) or greater than 1/2 values (for ν) can occur with these probability distributions. A reference solution C ref is obtained by considering an overkill solution, both in the deterministic and stochastic spaces. The deterministic overkill discretization is obtained by selecting 256 elements per side with quadratic elements and the stochastic one by means of two degree 10 polynomial chaos (Legendre polynomials)3 for S P 1 and S P 2 . The final reference solution was obtained using 40 PGD modes. The p convergence was conducted on a 4 × 4 elements mesh with p from 1 to 7, and two quartic Legendre polynomial chaos for E and ν. At each step of the procedure, the reduced model approximates a spectral stochastic solution which should converge exponentially provided that the solution is smooth [START_REF] Babuška | Galerkin finite element approximations of stochastic elliptic partial differential equations[END_REF] in both space and stochastic spaces. This is consistent with the convergence results for the expectation and standard deviation of C shown in figure 8. All the components of Exp(C) and Std(C) converge exponentially, which is also the case for the 2-norm of these tensors. As shown in the last section, the extrapolation error estimate cannot be used to drive the adaptation process if the accuracy of the individual components is mandatory. However, figure 9 shows that the robust deterministic behaviour of the estimator for the two-norm of C extends to the stochastic framework. It is important to note here that although Exp(E) and Exp(ν) are equal to their deterministic counterparts, the expectation of the stochastic effective tensor is not equal to the deterministic one when mean values are used for both E and ν. This stems from the fact that the solution is nonlinear in ν. Relative error (Expectation) Relative error (Std) Relative error (exact and estimated)

Exp(C11) Exp(C12) Exp(C13) Exp(C22) Exp(C23) Exp(C33) Exp(C) 2
Std(C11) Std(C12) Std(C13) Std(C22) Std(C23) Std(C33) Std(C) 2 Exp(C) Std(C) Figure 8: 2D example: Evolution of (Exp(C h,p ij ), Exp(C ref ij )) and (Std(C h,p ij ), Std(C ref ij ))
Exp(C) 2 Exact Exp(C) 2 Estimated Std(C) 2 Exact Std(C) 2 Estimated

Stochastic Feedback

Based on the verification of the convergence the method, it is now possible to evaluate adaptively the expectation and standard deviation of C. In this section, a ε r = 10 -2 threshold is considered for the relative stagnation, and ε a = 10 -5 for the absolute stagnation of both Exp(C) and Std(C). The evolution of the 2-norm of the expectation and standard deviation of the C are represented in figure 10. It can be seen that a p = 4 spatial discretization and 13 PGD modes are sufficient to fall into this range of accuracy. Individual components are not presented here, but their evolution is similar to the one of the norm, even for the off-diagonal components. Next, the accuracy of the stochastic model is assessed. Indeed, this model allows to obtain the homogenized properties of the porous material for any value of E and ν in their range of evolution. Thus, 10 000 random evaluations of E and ν are considered. From these random evaluations, C can be obtained either by using the reduced order model (C P GD with m = 13 modes) or by a direct computation of the overall properties (C M C for the same h and p). The former approach is almost costless once the reduced order model has been computed while the second (Monte Carlo) is computationally demanding. The relative 2-norm error between these two tensors (C P GD , C M C ) is evaluated for each sample and studied. The distribution of the error levels is then presented in figure 11. It can be seen that the error is smaller than 6 10 -3 , which is very accurate. Note also that the majority of the error levels are lower than 2 10 -3 . Finally, the reduced order model can be used to build the empirical probability density functions (pdf) of the components of C. The pdf of selected components of C are compared to those obtained from 10 000 direct Monte-Carlo simulations in figure 12. It can be seen that there is a perfect accordance between these curves 4 , which highlights the accuracy if the proposed approach. Now that the method has been validated for simple problems, we consider its application to biomechanics.

Application to a real 3D bone micro-structure

The method is now applied to a real bone micro-structure (see figure 13) whose volumic fraction is 0.1382. The geometry was obtained by segmenting a bone CT scan following the procedure proposed in [START_REF] Legrain | An X-FEM and level set computational approach for image-based modelling: Application to homogenization[END_REF]. This raw CT scan is available on the Internet through the BEL Repository managed by the Istituti Ortopedici Rizzoli, Bologna, Italy. The objective is to apply the proposed approach in order to compute the effective properties of this micro-structure. The material is considered as linear-elastic with typical material properties: 10 000 MPa Young's modulus and 0.4 Poisson's ratio. The resolution of the volume is 129 3 voxels, which corresponds to a level 7 octree. Kinematic uniform boundary conditions are still used here in order to extract the homogenized behaviour of the bone. Note that this choice is known to lead to over-estimated values of the effective properties [START_REF] Huet | Application of variational concepts to size effects in elastic heterogeneous bodies[END_REF].

Deterministic feedback

The feedback method is first considered in order to evaluate the deterministic homogenized tensor C. The computational mesh is made from a regular grid of 16 elements along each side of the cube, and the spatial approximation order ranges from p = 1 to p = 8. The following criteria are considered for stopping the refinement process: ε r = 5% for relative stagnation and ε a = 5 10 -4 for absolute stagnation. It can be seen in figure 16 that these conditions are met for p = 5, and that the method is able to evaluate efficiently the effective properties of the material. Some corresponding micro-modes are presented in figures 14 and 15 for illustration. The results are consistent with the literature, as we found C 33 = 727 Mpa, which is in agreement with the experimental results given in [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: Use of different densitometric measures[END_REF] (scattering between 200 and 900 Mpa for this bone density).

Stochastic feedback

Finally, the evaluation of the stochastic effective properties of the bone is presented. Similarly to the 2D numerical examples, uniform probability distributions are considered. A large scattering is chosen for the material properties in order to highlight the accuracy of the method. Young's modulus thus ranges from 5000 MPa to 15 000 Mpa, and Poisson's ratio from 0.32 to 0.48. Note that such a high Poisson's ratio would lead to some volumetric locking with low-order finite elements, which is not the case with high-order finite elements (see [START_REF] Babuška | On the rates of convergence of the finite element method[END_REF]). The convergence is attained when relative and absolute stagnation respectively falls below ε r = 5% and ε a = 5 10 -4 for both expectation and standard deviation. This target accuracy is attained for p = 6 and m = 13 PGD modes. The evolution of the 2-norm of C and some of its components during the p refinement process are depicted in figure 18. It can be observed that the largest dispersion is obtained for the diagonal part of the effective tensor. The expectation of C 33 (729 Mpa) is very similar to the one obtained in the deterministic case. The corresponding standard deviation (211 Mpa) is smaller than the one highlighted in [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: Use of different densitometric measures[END_REF], as the range of variation of the parameters was not fitted with any experimental data. Note also that any comparison with [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: Use of different densitometric measures[END_REF] is difficult since the results are based on multiple donors whose age can have a large influence [START_REF] Ladd | Finite-element modeling of trabecular bone: comparison with mechanical testing and determination of tissue modulus[END_REF], and also because this study also includes experimental bias [START_REF] Helgason | Mathematical relationships between bone density and mechanical properties: a literature review[END_REF] that are out of the scope of this study. The accuracy of the reduced model is now assessed through 100 samplings of the material properties. The reduced order model is particularized for these 100 samplings, and the results are compared with the corresponding direct calculations. The error distribution (C P GD , C M C ) is plotted in figure 19: it can be seen that it is smaller than 1.5 10 -3 , which highlights the accuracy of the reduced model. For a finer comparison between direct and reduced solutions, we now consider the mechanical fields associated to a given sample. The local relative error on the von-mises stress between direct and reduced solutions is presented in figure 17 where it can be seen that this local error is very small.

Finally, the reduced model is used to build an empirical probability density function: 10 000 samples are generated at no additional cost. The empirical pdf of the different components of the tensor are presented in figure 20.

Conclusion

An innovative numerical strategy has been proposed in order to compute adaptively the effective properties of complex micro-structures with uncertain material properties. The method relies on two main ingredients: (i) a high-order fictitious domain method approximation which enables to avoid the meshing burden with such complex geometries, and (ii) a model reduction technique based on the Proper Generalized Decomposition of the spectral stochastic representation of uncertainties (Linear Elasticity was considered here, with random Young's modulus and Poisson's ratio). The accuracy of the output was ensured even with large uncertainties on the material properties thanks to a feedback method which is used to monitor the convergence of the effective properties. The method has been verified on a simple 2D example which has shown that error estimates based on extensions were not sufficiently reliable for this type of application. However, it was seen that the use of a stagnation criterion allowed to obtain very accurate results as both expectation and standard deviation errors converge exponentially upon p refinement. The resulting stochastic effective tensor was seen to match Monte-Carlo's results with only a fraction of the computational cost. The method was finally applied to a real bone micro-structure whose stochastic effective properties could be obtained very efficiently. In order to improve the proposed strategy, goal oriented error estimation could be considered, as in [START_REF] Miranda Guedes | Preprocessing and postprocessing for materials based on the homogenization method with adaptive finite element methods[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] for the deterministic side, but also on the stochastic side. Note however that the influence of the stochastic discretization on the computational cost is very low thanks to the model reduction approach. Finally, multi-material and/or nonlinear materials could be considered, as well as the influence of a possible geometrical randomness.

Acknowledgements: G.Legrain and N. Takano would like to acknowledge partial support of the Erasmus Mundus BEAM program.

A. Derivation of the reduced model

Recall the minimization problem [START_REF] Legrain | Stability of incompressible formulations enriched with X-FEM[END_REF]. Two independent random variables are considered: ξ 1 for Young's modulus and ξ 2 for Poisson's ratio. Using a progressive definition of the PGD, we look for a solution in the form: u m = u m-1 + w(x)φ 1 (ξ 1 )φ 2 (ξ 2 ) where an element w(x)φ 1 (ξ 1 )φ 2 (ξ 2 ) is solution of problem [START_REF] Babuška | On the rates of convergence of the finite element method[END_REF]. Functional J(u m-1 + w φ 1 φ 2 )

Figure 2 :

 2 Figure 2: Definition of the unit cell of interest for the 2D verifications.

Figure 3 :

 3 Figure 3: 2D example: (a) reference errors (C h,p , C ref ) and (C h,p ij , C ref ij ) and (b) estimated errors (C h,p , C est ) and (C h,p ij , C est ij ) of the homogenized tensor C and its components (h convergence).

Figure 4 :

 4 Figure 4: 2D example: Evolution of the effectivity of the estimators

Figure 5 :

 5 Figure 5: 2D example: (a) reference errors (C h,p , C ref ) and (C h,p ij , C ref ij ) and (b) estimated errors (C h,p , C est ) and (C h,p ij , C est ij ) of the homogenized tensor C and its components (p convergence).

2 Figure 6 :

 26 Figure 6: 2D example: Evolution of the effectivity of the estimators

2 Figure 7 :

 27 Figure 7: 2D example: Feedback history of the homogenized tensor C. The thin vertical line indicates the fulfilment of the convergence criterion in the case of a 5 10 -2 relative stagnation criterion.

Figure 9 :

 9 Figure 9: 2D example: (a) comparison between reference and estimated errors on Exp(C) and Std(C) ( (Exp(C h,p ), Exp(C ref )) and (Std(C h,p ), Std(C ref ))) ; (b) effectivity index.

2 Figure 10 :

 210 Figure 10: 2D example: Evolution of Exp(C) and Std(C) during the feedback procedure.

Figure 11 :

 11 Figure 11: 2D example: Distribution of the error (C P GD , C M C ) for 10000 random evaluations.

Figure 12 :

 12 Figure 12: 2D example: Empirical pdf of the components of C: left, reduced model; right, Monte Carlo simulations (x-axis in MPa)

Figure 13 :

 13 Figure 13: Bone geometry.

Figure 14 :

 14 Figure 14: Micro-modes 1 (left) and 4 (right): displacement field

Figure 15 : 2 Figure 16 :

 15216 Figure 15: Micro-modes 1 (left) and 4 (right): Von-Mises stress field

Figure 17 :

 17 Figure 17: 3D example: Local relative error between direct and reduced solution for a given sample.

Figure 18 :

 18 Figure 18: 3D example: Evolution of some of the components of Exp(C) and associated dispersion during the feedback procedure.

Figure 19 :Figure 20 :

 1920 Figure 19: 3D example: Distribution of the error (C P GD , C M C ) for 100 random evaluations.

Correspondence to: G. Legrain. E-Mail: gregory.legrain@ec-nantes.fr

Note that the method itself has already been applied to the resolution of image-based problems, see for instance[START_REF] Lamari | Routes for Efficient Computational Homogenization of Nonlinear Materials Using the Proper Generalized Decompositions[END_REF] and[START_REF] Giraldi | Tensor-based methods for numerical homogenization from high-resolution images[END_REF]

Note that high order approximations are used both in the deterministic and stochastic space.

It is also the case for the 3 other components that are not presented here.

can be expanded into:

where the last terms in parenthesis are denoted as

In addition, we can write:

In our case, the Hooke tensor can be separated as:

with:

where ξ 1 = E and ξ 2 = ν. Thus, we have:

All the integrals can be separated, leading to the following functional:

The expression of R is not given here, but can be obtained following the same rationale. Necessary conditions of optimality for a triplet (w, φ 1 , φ 2 ) are the three Galerkin criteria:

w = arg min

Thus, the nonlinear equation is solved by means of a fixed point algorithm which reduces in the following three linear problems: P 1 : φ 2 and w are known. φ 1 is solution of problem [START_REF] Legrain | On the use of the extended finite element method with quadtree/octree meshes[END_REF]. Stochastic problem of dimension dim(S N 1 ); P 2 : φ 1 and w are known. φ 2 is solution of problem [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF]. Stochastic problem of dimension dim(S N 2 ); P 3 : φ 1 and φ 2 are known. w is solution of problem [START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF]. Stochastic problem of dimension dim(V N );

Additional update strategies can be considered in order to improve the efficiency of the algorithm, the interested reader can refer to [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: Invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF][START_REF] Nouy | A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations[END_REF].